text:
//| Copyright Inria May 2015
//| This project has received funding from the European Research Council (ERC) under
//| the European Union's Horizon 2020 research and innovation programme (grant
//| agreement No 637972) - see http://www.resibots.eu
//|
//| Contributor(s):
//|   - Jean-Baptiste Mouret (jean-baptiste.mouret@inria.fr)
//|   - Antoine Cully (antoinecully@gmail.com)
//|   - Konstantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
//|   - Federico Allocati (fede.allocati@gmail.com)
//|   - Vaios Papaspyros (b.papaspyros@gmail.com)
//|   - Roberto Rama (bertoski@gmail.com)
//|
//| This software is a computer library whose purpose is to optimize continuous,
//| black-box functions. It mainly implements Gaussian processes and Bayesian
//| optimization.
//| Main repository: http://github.com/resibots/limbo
//| Documentation: http://www.resibots.eu/limbo
//|
//| This software is governed by the CeCILL-C license under French law and
//| abiding by the rules of distribution of free software.  You can  use,
//| modify and/ or redistribute the software under the terms of the CeCILL-C
//| license as circulated by CEA, CNRS and INRIA at the following URL
//| "http://www.cecill.info".
//|
//| As a counterpart to the access to the source code and  rights to copy,
//| modify and redistribute granted by the license, users are provided only
//| with a limited warranty  and the software's author,  the holder of the
//| economic rights,  and the successive licensors  have only  limited
//| liability.
//|
//| In this respect, the user's attention is drawn to the risks associated
//| with loading,  using,  modifying and/or developing or reproducing the
//| software by the user in light of its specific status of free software,
//| that may mean  that it is complicated to manipulate,  and  that  also
//| therefore means  that it is reserved for developers  and  experienced
//| professionals having in-depth computer knowledge. Users are therefore
//| encouraged to load and test the software's suitability as regards their
//| requirements in conditions enabling the security of their systems and/or
//| data to be ensured and,  more generally, to use and operate it in the
//| same conditions as regards security.
//|
//| The fact that you are presently reading this means that you have had
//| knowledge of the CeCILL-C license and that you accept its terms.
//|
#ifndef LIMBO_BAYES_OPT_IMGPO_HPP
#define LIMBO_BAYES_OPT_IMGPO_HPP

#include <algorithm>
#include <cmath>
#include <iostream>
#include <iterator>

#include <boost/parameter/aux_/void.hpp>

#include <Eigen/Core>

#include <limbo/bayes_opt/bo_base.hpp>
#include <limbo/tools/macros.hpp>
#include <limbo/tools/random_generator.hpp>

namespace limbo {
    namespace defaults {
        struct bayes_opt_imgpo {
            BO_PARAM(bool, hp_opt, false);
        };
    }

    namespace bayes_opt {
        namespace experimental {
            struct TreeNode {
                std::vector<Eigen::VectorXd> x_max, x_min, x, f;
                std::vector<bool> leaf, samp;
            };
            // clang-format off
        template <class Params,
          class A1 = boost::parameter::void_,
          class A2 = boost::parameter::void_,
          class A3 = boost::parameter::void_,
          class A4 = boost::parameter::void_,
          class A5 = boost::parameter::void_>
            // clang-format on
            // Bayesian Optimization with Exponential Convergence (NIPS 2015 paper)
            // Kenji Kawaguchi, Leslie Pack Kaelbling, Tomás Lozano-Pérez
            // http://papers.nips.cc/paper/5715-bayesian-optimization-with-exponential-convergence
            class IMGPO : public limbo::bayes_opt::BoBase<Params, A1, A2, A3, A4, A5> {
            public:
                using base_t = BoBase<Params, A1, A2, A3, A4, A5>;
                using model_t = typename base_t::model_t;
                using acquisition_function_t = typename base_t::acquisition_function_t;

                template <typename StateFunction, typename AggregatorFunction = FirstElem>
                void optimize(const StateFunction& sfun, const AggregatorFunction& afun = AggregatorFunction(), bool reset = true)
                {
                    this->_init(sfun, afun, reset);
                    size_t h_upper = 1000;
                    // Init tree
                    if (this->_total_iterations == 0)
                        _init_tree(h_upper);

                    // Init model
                    _model = model_t(StateFunction::dim_in(), StateFunction::dim_out());

                    // Init root
                    _tree[0].x_max.push_back(Eigen::VectorXd::Ones(StateFunction::dim_in()));
                    _tree[0].x_min.push_back(Eigen::VectorXd::Zero(StateFunction::dim_in()));
                    _tree[0].x.push_back(Eigen::VectorXd::Ones(StateFunction::dim_in()) * 0.5);
                    _tree[0].f.push_back(sfun(_tree[0].x[0]));
                    _tree[0].leaf.push_back(true);
                    _tree[0].samp.push_back(true);

                    this->add_new_sample(_tree[0].x[0], _tree[0].f[0]);
                    double LB = afun(_tree[0].f[0]);
                    double inf = std::numeric_limits<double>::infinity();

                    size_t depth_T = 0, M = 1;
                    double rho_avg = 0, rho_bar = 0;
                    size_t t = 0, XI_max = 4, split_n = 1;
                    int xi_max = 0;
                    double XI = 1;
                    double LB_old = LB;
                    // N = N+1
                    this->_current_iteration++;
                    this->_total_iterations++;

                    while (this->_samples.size() == 0 || !this->_stop(*this, afun)) {
                        std::vector<int> i_max(depth_T + 1, -1);
                        std::vector<double> b_max(depth_T + 1, -inf);
                        size_t h_max = depth_T + 1;
                        double b_hi_max = -inf;
                        t = t + 1;

                        // Steps (i)-(ii)
                        for (size_t h = 0; h <= depth_T; h++) {
                            if (h >= h_max)
                                break;

                            while (true) {
                                for (size_t i = 0; i < _tree[h].x.size(); i++) {
                                    if (_tree[h].leaf[i] == true) {
                                        double b_hi = afun(_tree[h].f[i]);
                                        if (b_hi > b_hi_max) {
                                            b_hi_max = b_hi;
                                            i_max[h] = i;
                                            b_max[h] = b_hi;
                                        }
                                    }
                                }
                                if (i_max[h] == -1)
                                    break;

                                if (_tree[h].samp[i_max[h]] == true) {
                                    break;
                                }
                                else {
                                    Eigen::VectorXd xxx = _tree[h].x[i_max[h]];

                                    auto tmp_sample = sfun(xxx);
                                    _tree[h].samp[i_max[h]] = true;
                                    this->add_new_sample(xxx, tmp_sample);

                                    // N = N+1
                                    this->_current_iteration++;
                                    this->_total_iterations++;
                                }
                            }
                        }

                        // Steps (iii)
                        for (size_t h = 0; h <= depth_T; h++) {
                            if (h >= h_max)
                                break;
                            if (i_max[h] != -1) {
                                int xi = -1;
                                for (size_t h2 = h + 1; h2 <= std::min(depth_T, h + std::min((size_t)std::ceil(XI), XI_max)); h2++) {
                                    if (i_max[h2] != -1) {
                                        xi = h2 - h;
                                        break;
                                    }
                                }

                                double z_max = -inf;
                                size_t M2 = M;
                                if (xi != -1) {
                                    std::vector<TreeNode> tmp_tree = std::vector<TreeNode>(h + xi + 1);
                                    tmp_tree[h].x_max.push_back(_tree[h].x_max[i_max[h]]);
                                    tmp_tree[h].x_min.push_back(_tree[h].x_min[i_max[h]]);
                                    tmp_tree[h].x.push_back(_tree[h].x[i_max[h]]);

                                    M2 = M;
                                    for (size_t h2 = h; h2 <= h + xi - 1; h2++) {
                                        for (size_t ii = 0; ii < std::pow(3, h2 - h); ii++) {
                                            if (ii >= tmp_tree[h].x.size())
                                                break;
                                            Eigen::VectorXd xx = tmp_tree[h].x[ii];
                                            Eigen::VectorXd to_split = tmp_tree[h2].x_max[ii].array() - tmp_tree[h2].x_min[ii].array();
                                            size_t tmp, splitd;
                                            to_split.maxCoeff(&splitd, &tmp);
                                            auto x_g = xx, x_d = xx;
                                            x_g(splitd) = (5 * tmp_tree[h2].x_min[ii](splitd) + tmp_tree[h2].x_max[ii](splitd)) / 6.0;
                                            x_d(splitd) = (tmp_tree[h2].x_min[ii](splitd) + 5 * tmp_tree[h2].x_max[ii](splitd)) / 6.0;

                                            // TO-DO: Properly handle bl_samples etc
                                            _model.compute(this->_samples, this->_observations);
                                            acquisition_function_t acqui_g(_model, M2);
                                            z_max = std::max(z_max, acqui_g(x_g, afun));
                                            M2++;

                                            _model.compute(this->_samples, this->_observations);
                                            acquisition_function_t acqui_d(_model, M2);
                                            z_max = std::max(z_max, acqui_d(x_d, afun));
                                            M2++;

                                            if (z_max >= b_max[h + xi])
                                                break;

                                            tmp_tree[h2 + 1].x.push_back(x_g);
                                            Eigen::VectorXd newmin = tmp_tree[h2].x_min[ii];
                                            tmp_tree[h2 + 1].x_min.push_back(newmin);
                                            Eigen::VectorXd newmax = tmp_tree[h2].x_max[ii];
                                            newmax(splitd) = (2 * tmp_tree[h2].x_min[ii](splitd) + tmp_tree[h2].x_max[ii](splitd)) / 3.0;
                                            tmp_tree[h2 + 1].x_max.push_back(newmax);

                                            tmp_tree[h2 + 1].x.push_back(x_d);
                                            Eigen::VectorXd newmax2 = tmp_tree[h2].x_max[ii];
                                            tmp_tree[h2 + 1].x_max.push_back(newmax2);
                                            Eigen::VectorXd newmin2 = tmp_tree[h2].x_min[ii];
                                            newmin2(splitd) = (tmp_tree[h2].x_min[ii](splitd) + 2 * tmp_tree[h2].x_max[ii](splitd)) / 3.0;
                                            tmp_tree[h2 + 1].x_min.push_back(newmin2);

                                            tmp_tree[h2 + 1].x.push_back(xx);
                                            Eigen::VectorXd newmin3 = tmp_tree[h2].x_min[ii];
                                            newmin3(splitd) = (2 * tmp_tree[h2].x_min[ii](splitd) + tmp_tree[h2].x_max[ii](splitd)) / 3.0;
                                            tmp_tree[h2 + 1].x_min.push_back(newmin3);
                                            Eigen::VectorXd newmax3 = tmp_tree[h2].x_max[ii];
                                            newmax3(splitd) = (tmp_tree[h2].x_min[ii](splitd) + 2 * tmp_tree[h2].x_max[ii](splitd)) / 3.0;
                                            tmp_tree[h2 + 1].x_max.push_back(newmax3);
                                        }
                                        if (z_max >= b_max[h + xi])
                                            break;
                                    }
                                }

                                if (xi != -1 && z_max < b_max[h + xi]) {
                                    M = M2;
                                    i_max[h] = -1;
                                    xi_max = std::max(xi, xi_max);
                                }
                            }
                        }

                        // Steps (iv)-(v)
                        double b_hi_max_2 = -inf, rho_t = 0.0;
                        for (size_t h = 0; h <= depth_T; h++) {
                            if (h >= h_max)
                                break;
                            if (i_max[h] != -1 && b_max[h] > b_hi_max_2) {
                                rho_t += 1.0;
                                depth_T = std::max(depth_T, h + 1);
                                split_n++;
                                _tree[h].leaf[i_max[h]] = 0;

                                Eigen::VectorXd xx = _tree[h].x[i_max[h]];
                                Eigen::VectorXd to_split = _tree[h].x_max[i_max[h]].array() - _tree[h].x_min[i_max[h]].array();
                                size_t tmp, splitd;
                                to_split.maxCoeff(&splitd, &tmp);
                                auto x_g = xx, x_d = xx;
                                x_g(splitd) = (5 * _tree[h].x_min[i_max[h]](splitd) + _tree[h].x_max[i_max[h]](splitd)) / 6.0;
                                x_d(splitd) = (_tree[h].x_min[i_max[h]](splitd) + 5 * _tree[h].x_max[i_max[h]](splitd)) / 6.0;

                                // left node
                                _tree[h + 1].x.push_back(x_g);
                                // TO-DO: Properly handle bl_samples etc
                                _model.compute(this->_samples, this->_observations);
                                acquisition_function_t acqui_g(_model, M);
                                double UCB = acqui_g(x_g, afun);
                                Eigen::VectorXd fsample_g;
                                if ((UCB - LB) < 1e-6) {
                                    Eigen::VectorXd mu;
                                    double sigma;
                                    std::tie(mu, sigma) = _model.query(x_g);
                                    // TO-DO: we should fix this somehow for general acquis
                                    double var = (std::sqrt(2.0 * std::log(std::pow(M_PI, 2.0) * std::pow(M, 2.0) / (12.0 * 0.05))) + 0.2) * std::sqrt(sigma);
                                    M++;
                                    fsample_g = mu.array() + var;
                                    _tree[h + 1].samp.push_back(false);
                                }
                                else {
                                    fsample_g = sfun(x_g);
                                    _tree[h + 1].samp.push_back(true);

                                    this->add_new_sample(x_g, fsample_g);

                                    b_hi_max_2 = std::max(b_hi_max_2, fsample_g(0));

                                    // N = N+1
                                    this->_current_iteration++;
                                    this->_total_iterations++;
                                }

                                _tree[h + 1].f.push_back(fsample_g);

                                Eigen::VectorXd newmin = _tree[h].x_min[i_max[h]];
                                _tree[h + 1].x_min.push_back(newmin);
                                Eigen::VectorXd newmax = _tree[h].x_max[i_max[h]];
                                newmax(splitd) = (2 * _tree[h].x_min[i_max[h]](splitd) + _tree[h].x_max[i_max[h]](splitd)) / 3.0;
                                _tree[h + 1].x_max.push_back(newmax);
                                _tree[h + 1].leaf.push_back(true);

                                // right node
                                _tree[h + 1].x.push_back(x_d);
                                // TO-DO: Properly handle bl_samples etc
                                _model.compute(this->_samples, this->_observations);
                                acquisition_function_t acqui_d(_model, M);
                                double UCB2 = acqui_d(x_d, afun);
                                Eigen::VectorXd fsample_d;
                                if ((UCB2 - LB) < 1e-6) {
                                    Eigen::VectorXd mu;
                                    double sigma;
                                    std::tie(mu, sigma) = _model.query(x_d);
                                    // TO-DO: we should fix this somehow for general acquis
                                    double var = (std::sqrt(2.0 * std::log(std::pow(M_PI, 2.0) * std::pow(M, 2.0) / (12.0 * 0.05))) + 0.2) * std::sqrt(sigma);
                                    M++;
                                    fsample_d = mu.array() + var;
                                    _tree[h + 1].samp.push_back(false);
                                }
                                else {
                                    fsample_d = sfun(x_d);
                                    _tree[h + 1].samp.push_back(true);

                                    this->add_new_sample(x_d, fsample_d);

                                    b_hi_max_2 = std::max(b_hi_max_2, fsample_d(0));

                                    // N = N+1
                                    this->_current_iteration++;
                                    this->_total_iterations++;
                                }

                                _tree[h + 1].f.push_back(fsample_d);

                                Eigen::VectorXd newmax2 = _tree[h].x_max[i_max[h]];
                                _tree[h + 1].x_max.push_back(newmax2);
                                Eigen::VectorXd newmin2 = _tree[h].x_min[i_max[h]];
                                newmin2(splitd) = (_tree[h].x_min[i_max[h]](splitd) + 2 * _tree[h].x_max[i_max[h]](splitd)) / 3.0;
                                _tree[h + 1].x_min.push_back(newmin2);
                                _tree[h + 1].leaf.push_back(true);

                                // central node
                                _tree[h + 1].x.push_back(xx);
                                _tree[h + 1].f.push_back(_tree[h].f[i_max[h]]);
                                _tree[h + 1].samp.push_back(true);
                                Eigen::VectorXd newmin3 = _tree[h].x_min[i_max[h]];
                                Eigen::VectorXd newmax3 = _tree[h].x_max[i_max[h]];
                                newmin3(splitd) = (2 * _tree[h].x_min[i_max[h]](splitd) + _tree[h].x_max[i_max[h]](splitd)) / 3.0;
                                newmax3(splitd) = (_tree[h].x_min[i_max[h]](splitd) + 2 * _tree[h].x_max[i_max[h]](splitd)) / 3.0;
                                _tree[h + 1].x_min.push_back(newmin3);
                                _tree[h + 1].x_max.push_back(newmax3);
                                _tree[h + 1].leaf.push_back(true);

                                // update LB
                                LB = afun(this->best_observation());

                                // std::cout << t << " N= " << this->_current_iteration << " n= " << split_n << " fmax_hat= " << LB << " rho= " << rho_bar << " xi= " << xi_max << " h= " << depth_T << std::endl;
                            }
                        }

                        // Finalizing iteration
                        rho_avg = (rho_avg * (t - 1) + rho_t) / t;
                        rho_bar = std::max(rho_bar, rho_avg);
                        // update XI
                        if (std::abs(LB_old - LB) < 1e-6)
                            XI = std::max(XI - 0.5, 1.0);
                        else
                            XI = XI + 4.0;
                        LB_old = LB;
                    }

                    if (Params::bayes_opt_imgpo::hp_opt())
                        _model.optimize_hyperparams();
                }

                template <typename AggregatorFunction = FirstElem>
                const Eigen::VectorXd& best_observation(const AggregatorFunction& afun = AggregatorFunction()) const
                {
                    auto rewards = std::vector<double>(this->_observations.size());
                    std::transform(this->_observations.begin(), this->_observations.end(), rewards.begin(), afun);
                    auto max_e = std::max_element(rewards.begin(), rewards.end());
                    return this->_observations[std::distance(rewards.begin(), max_e)];
                }

                template <typename AggregatorFunction = FirstElem>
                const Eigen::VectorXd& best_sample(const AggregatorFunction& afun = AggregatorFunction()) const
                {
                    auto rewards = std::vector<double>(this->_observations.size());
                    std::transform(this->_observations.begin(), this->_observations.end(), rewards.begin(), afun);
                    auto max_e = std::max_element(rewards.begin(), rewards.end());
                    return this->_samples[std::distance(rewards.begin(), max_e)];
                }

                const model_t& model() const { return _model; }

            protected:
                void _init_tree(size_t h_max = 1000)
                {
                    _tree.clear();
                    _tree = std::vector<TreeNode>(h_max);
                }

                model_t _model;
                std::vector<TreeNode> _tree;
            };
        }
    }
}

#endif

ext: hpp


text:
/*
 * Copyright (C) 2011-2018 Daniel Scharrer
 *
 * This software is provided 'as-is', without any express or implied
 * warranty.  In no event will the author(s) be held liable for any damages
 * arising from the use of this software.
 *
 * Permission is granted to anyone to use this software for any purpose,
 * including commercial applications, and to alter it and redistribute it
 * freely, subject to the following restrictions:
 *
 * 1. The origin of this software must not be misrepresented; you must not
 *    claim that you wrote the original software. If you use this software
 *    in a product, an acknowledgment in the product documentation would be
 *    appreciated but is not required.
 * 2. Altered source versions must be plainly marked as such, and must not be
 *    misrepresented as being the original software.
 * 3. This notice may not be removed or altered from any source distribution.
 */

#include "crypto/checksum.hpp"

#include <ios>
#include <cstring>
#include <ostream>
#include <iomanip>

#include <boost/range/size.hpp>

namespace crypto {

bool checksum::operator==(const checksum & other) const {

	if(other.type != type) {
		return false;
	}

	switch(type) {
		case None: return true;
		case Adler32: return (adler32 == other.adler32);
		case CRC32: return (crc32 == other.crc32);
		case MD5: return !memcmp(md5, other.md5, sizeof(md5));
		case SHA1: return !memcmp(sha1, other.sha1, sizeof(sha1));
		default: return false;
	};
}

} // namespace crypto

NAMES(crypto::checksum_type, "Checksum Type",
	"None",
	"Adler32",
	"CRC32",
	"MD5",
	"SHA-1",
)

std::ostream & operator<<(std::ostream & os, const crypto::checksum & checksum) {

	std::ios_base::fmtflags old_fmtflags = os.flags();

	os << checksum.type << ' ';

	switch(checksum.type) {
		case crypto::None: {
			os << "(no checksum)";
			break;
		}
		case crypto::Adler32: {
			os << "0x" << std::hex << std::setw(8) << checksum.adler32;
			break;
		}
		case crypto::CRC32: {
			os << "0x" << std::hex << std::setw(8) << checksum.crc32;
			break;
		}
		case crypto::MD5: {
			for(size_t i = 0; i < size_t(boost::size(checksum.md5)); i++) {
				os << std::setfill('0') << std::hex << std::setw(2) << int(boost::uint8_t(checksum.md5[i]));
			}
			break;
		}
		case crypto::SHA1: {
			for(size_t i = 0; i < size_t(boost::size(checksum.sha1)); i++) {
				os << std::setfill('0') << std::hex << std::setw(2) << int(boost::uint8_t(checksum.sha1[i]));
			}
			break;
		}
	}

	os.setf(old_fmtflags, std::ios_base::basefield);

	return os;
}

ext: cpp


text: #include <boost/units/physical_dimensions/time.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include <boost/atomic.hpp>
#include <thread>
#include <iostream>

boost::atomic<int> a{0};
void thread()
{
  a.fetch_add(1, boost::memory_order_relaxed);
}

int main(int argc, char *argv[]) {
    std::cout << "CMAKE Boost -- Hello world" << std::endl;
    std::cout << " - App is built using " << BUILD_TYPE << std::endl;

    std::thread t1{thread};
    std::thread t2{thread};
    t1.join();
    t2.join();
    std::cout << a << '\n';

    (void)argc;
    (void)argv;

    return 0;
}

ext: cc


text:
/*
 *            Copyright 2009-2016 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#include <votca/ctp/apolarsite.h>
#include <boost/math/special_functions/round.hpp>
#include <boost/format.hpp>
#include <fstream>
#include <string>

namespace votca { namespace ctp {

APolarSite::APolarSite(APolarSite *templ, bool do_depolarize)
    : _id(templ->_id), _name(templ->_name), _isVirtual(templ->_isVirtual),
       _pos(templ->_pos),

      _locX(templ->_locX), _locY(templ->_locY), _locZ(templ->_locZ),

      _top(templ->_top), _seg(templ->_seg), _frag(templ->_frag),

      _Qs(templ->_Qs), _rank(templ->_rank), _Ps(templ->_Ps),
      Pxx(templ->Pxx), Pxy(templ->Pxy), Pxz(templ->Pxz), Pyy(templ->Pyy),
      Pyz(templ->Pyz), Pzz(templ->Pzz), pax(templ->pax), pay(templ->pay),
      paz(templ->paz), eigenpxx(templ->eigenpxx), eigenpyy(templ->eigenpyy),
      eigenpzz(templ->eigenpzz), eigendamp(templ->eigendamp),

      Q00(templ->Q00), Q1x(templ->Q1x), Q1y(templ->Q1y), Q1z(templ->Q1z),
      Q20(templ->Q20), Q21c(templ->Q21c), Q21s(templ->Q21s), Q22c(templ->Q22c),
      Q22s(templ->Q22s), Qxx(templ->Qxx), Qxy(templ->Qxy), Qxz(templ->Qxz),
      Qyy(templ->Qyy), Qyz(templ->Qyz), Qzz(templ->Qzz),

      U1x(templ->U1x), U1y(templ->U1y), U1z(templ->U1z),
      FPx(templ->FPx), FPy(templ->FPy), FPz(templ->FPz),
      FUx(templ->FUx), FUy(templ->FUy), FUz(templ->FUz),

     _resolution(templ->_resolution),PhiP(templ->PhiP), PhiU(templ->PhiU) {

    if (do_depolarize) this->Depolarize();
}

void APolarSite::ImportFrom(APolarSite *templ, string tag) {

    _pos = templ->getPos();
    _resolution = templ->getResolution();

    if (tag == "basic") {
        _Qs[0] = templ->getQs(-1);
        _Qs[1] = templ->getQs(0);
        _Qs[2] = templ->getQs(1);
        _rank  = templ->_rank;
        _Ps    = templ->_Ps;
    }

    if (tag == "full") {
        _locX = templ->_locX;
        _locY = templ->_locY;
        _locZ = templ->_locZ;
        _top  = templ->_top;
        _seg  = templ->_seg;
        _frag = templ->_frag;
        _rank = templ->_rank;
        _Qs   = templ->_Qs;
        _Ps   = templ->_Ps;
        _id   = templ->_id;
        _name = templ->_name;
    }
}

void APolarSite::Rotate(const matrix &rot, const vec &refPos) {

    vec dir = _pos - refPos;
    dir = rot * dir;
    _pos = refPos + dir;

    _locX = rot.getCol(0);
    _locY = rot.getCol(1);
    _locZ = rot.getCol(2);

    // Rotate multipoles into global frame >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    for (int state = -1; state < 2; state++) {

        //   0    1    2    3    4    5    6    7    8    9    10   ...
        //   Q00  Q10  Q1c  Q1s  Q20  Q21c Q21s Q22c Q22s Q30  Q31c ...

        matrix R = rot;
        matrix R_T = matrix(R.getRow(0),R.getRow(1),R.getRow(2));

        // Transform polarizability tensor into global frame
        matrix P_Global = R * _Ps[state+1] * R_T;
        _Ps[state+1] = P_Global;

        // Any multipoles for this charge state available?
        if (_Qs[state+1].size() < 1) { continue; }

        // Transform dipole moment into global frame
        if (_Qs[state+1].size() > 1) {

            double Qz = _Qs[state+1][1];
            double Qx = _Qs[state+1][2];
            double Qy = _Qs[state+1][3];

            vec d = vec(Qx, Qy, Qz);
            d = R * d;

            _Qs[state+1][1] = d.getZ();
            _Qs[state+1][2] = d.getX();
            _Qs[state+1][3] = d.getY();
        }

        // Transform quadrupole moment into global frame
        if (_Qs[state+1].size() > 4) {

            double Qzz =      _Qs[state+1][4];
            double Qxx = -0.5*_Qs[state+1][4] + 0.5*sqrt(3)*_Qs[state+1][7];
            double Qyy = -0.5*_Qs[state+1][4] - 0.5*sqrt(3)*_Qs[state+1][7];

            double Qxy =  0.5*sqrt(3)*_Qs[state+1][8];
            double Qxz =  0.5*sqrt(3)*_Qs[state+1][5];
            double Qyz =  0.5*sqrt(3)*_Qs[state+1][6];

            matrix Q = matrix(vec(Qxx,Qxy,Qxz),
                              vec(Qxy,Qyy,Qyz),
                              vec(Qxz,Qyz,Qzz));

            matrix Q_Global  = R * Q * R_T;

            /* if (this->getId() == 1) {
                cout << endl;
                cout << "  " << Q_Global.get(0,0);
                cout << "  " << Q_Global.get(0,1);
                cout << "  " << Q_Global.get(0,2);
                cout << endl;
                cout << "  " << Q_Global.get(1,0);
                cout << "  " << Q_Global.get(1,1);
                cout << "  " << Q_Global.get(1,2);
                cout << endl;
                cout << "  " << Q_Global.get(2,0);
                cout << "  " << Q_Global.get(2,1);
                cout << "  " << Q_Global.get(2,2);
                cout << endl;
            }                                   */

            _Qs[state+1][4] =               Q_Global.get(2,2);  // Q20
            _Qs[state+1][5] = 2 / sqrt(3) * Q_Global.get(0,2);  // Q21c
            _Qs[state+1][6] = 2 / sqrt(3) * Q_Global.get(1,2);  // Q21s
            _Qs[state+1][7] = 1 / sqrt(3) *(Q_Global.get(0,0)
                                          - Q_Global.get(1,1)); // Q22c
            _Qs[state+1][8] = 2 / sqrt(3) * Q_Global.get(0,1);  // Q22s
        }

    }
    // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

}

bool APolarSite::getIsActive(bool estatics_only) {
    bool is_active = false;
    bool is_charged = IsCharged();
    bool is_polarizable = IsPolarizable();

    if (estatics_only) is_active = is_charged;
    else is_active = is_charged || is_polarizable;

    return is_active;
}

bool APolarSite::IsCharged() {
    bool is_charged = false;
    // Tolerances
    double q_tol = 1e-9; // [e]
    double d_tol = 1e-9; // [enm]
    double Q_tol = 1e-9; // [enm^2]
    // Magnitudes
    double q_mag = sqrt(Q00*Q00);
    double d_mag = sqrt(Q1x*Q1x + Q1y*Q1y + Q1z*Q1z);
    double Q_mag = sqrt(Q20*Q20 + Q22c*Q22c + Q22s*Q22s + Q21c*Q21c + Q21s*Q21s);
    // Compare
    if (q_mag>q_tol) is_charged = true;
    if (_rank > 0 && d_mag>d_tol) is_charged = true;
    if (_rank > 1 && Q_mag>Q_tol) is_charged = true;
    return is_charged;
}

bool APolarSite::IsPolarizable() {
    bool is_polarizable = false;
    // Tolerances
    double p_tol = 1e-9; // [nm^3]
    // Compare
    if (getIsoP() > p_tol) is_polarizable = true;
    return is_polarizable;
}

void APolarSite::Translate(const vec &shift) {
    _pos += shift;
    return;
}

void APolarSite::Charge(int state) {
    int idx = state + 1;
    // Adjust polarizability to charge state
    Pxx = _Ps[idx].get(0,0);
    Pxy = _Ps[idx].get(0,1);
    Pxz = _Ps[idx].get(0,2);
    Pyy = _Ps[idx].get(1,1);
    Pyz = _Ps[idx].get(1,2);
    Pzz = _Ps[idx].get(2,2);
    // Calculate principal axes
    matrix polarity = matrix( vec(Pxx,Pxy,Pxz),
                              vec(Pxy,Pyy,Pyz),
                              vec(Pxz,Pyz,Pzz) );
    matrix::eigensystem_t eigensystem_polarity;
    polarity.SolveEigensystem(eigensystem_polarity);
    pax         = eigensystem_polarity.eigenvecs[0];
    pay         = eigensystem_polarity.eigenvecs[1];
    paz         = eigensystem_polarity.eigenvecs[2];
    eigenpxx    = eigensystem_polarity.eigenvalues[0];
    eigenpyy    = eigensystem_polarity.eigenvalues[1];
    eigenpzz    = eigensystem_polarity.eigenvalues[2];
    eigendamp =  (eigenpxx > eigenpyy) ?
                     ((eigenpxx > eigenpzz) ? eigenpxx : eigenpzz)
                   : ((eigenpyy > eigenpzz) ? eigenpyy : eigenpzz);
    //cout << endl << "eigensystem ...";
    //cout << endl << pax << " --- " << eigenpxx;
    //cout << endl << pay << " --- " << eigenpyy;
    //cout << endl << paz << " --- " << eigenpzz;

    // Adjust multipole moments to charge state
    Q00 = _Qs[idx][0];

    if (_rank > 0) {
        Q1z = _Qs[idx][1];   // |
        Q1x = _Qs[idx][2];   // |-> NOTE: order z - x - y
        Q1y = _Qs[idx][3];   // |
    }
    else {
        Q1z = Q1x = Q1y = 0.0;
    }
    if (_rank > 1) {
        // Spherical tensor
        Q20  = _Qs[idx][4];
        Q21c = _Qs[idx][5];
        Q21s = _Qs[idx][6];
        Q22c = _Qs[idx][7];
        Q22s = _Qs[idx][8];

        // Cartesian tensor * 1/3
        Qzz =      Q20;
        Qxx = -0.5*Q20 + 0.5*sqrt(3)*Q22c;
        Qyy = -0.5*Q20 - 0.5*sqrt(3)*Q22c;
        Qxy =          + 0.5*sqrt(3)*Q22s;
        Qxz =          + 0.5*sqrt(3)*Q21c;
        Qyz =          + 0.5*sqrt(3)*Q21s;

        Qzz *= 1./3.;
        Qxx *= 1./3.;
        Qyy *= 1./3.;
        Qxy *= 1./3.;
        Qxz *= 1./3.;
        Qyz *= 1./3.;
    }
    else {
        Q20 = Q21c = Q21s = Q22c = Q22s = 0.0;
        Qzz = Qxx = Qyy = Qxy = Qxz = Qyz = 0.0;
    }
    return;
}

void APolarSite::ChargeDelta(int state1, int state2) {
    int idx1 = state1 + 1;
    int idx2 = state2 + 1;
    // Adjust polarizability to charge state
    Pxx = 0.0;
    Pxy = 0.0;
    Pxz = 0.0;
    Pyy = 0.0;
    Pyz = 0.0;
    Pzz = 0.0;
    // Adjust multipole moments to charge state
    Q00 = _Qs[idx2][0] - _Qs[idx1][0];
    if (_rank > 0) {
        Q1z = _Qs[idx2][1] - _Qs[idx1][1];   // |
        Q1x = _Qs[idx2][2] - _Qs[idx1][2];   // |-> NOTE: order z - x - y
        Q1y = _Qs[idx2][3] - _Qs[idx1][3];   // |
    }
    if (_rank > 1) {
        Q20  = _Qs[idx2][4] - _Qs[idx1][4];
        Q21c = _Qs[idx2][5] - _Qs[idx1][5];
        Q21s = _Qs[idx2][6] - _Qs[idx1][6];
        Q22c = _Qs[idx2][7] - _Qs[idx1][7];
        Q22s = _Qs[idx2][8] - _Qs[idx1][8];
    }
    return;
}

double APolarSite::getProjP(vec &dir) {
    double alpha_proj = fabs(dir * pax) * fabs(dir * pax) * eigenpxx
         + fabs(dir * pay) * fabs(dir * pay) * eigenpyy
         + fabs(dir * paz) * fabs(dir * paz) * eigenpzz;

    // Mix ?
    //alpha_proj = 0.5* (1./3. * (Pxx+Pyy+Pzz)) + 0.5*alpha_proj;

    //assert("APS::getProjP THIS FEATURE SHOULD NOT BE USED" == "" && false);
    return alpha_proj;
}

void APolarSite::Induce(double wSOR) {
    // SUCCESSIVE OVERRELAXATION
    U1_Hist.push_back( vec(U1x,U1y,U1z) ); // Remember all previous moments
    //U1_Hist[0] = vec(U1x,U1y,U1z);           // Remember previous moment
    U1x = (1 - wSOR) * U1x + wSOR * ( - Pxx * (FPx + FUx) - Pxy * (FPy + FUy) - Pxz * (FPz + FUz) );
    U1y = (1 - wSOR) * U1y + wSOR * ( - Pxy * (FPx + FUx) - Pyy * (FPy + FUy) - Pyz * (FPz + FUz) );
    U1z = (1 - wSOR) * U1z + wSOR * ( - Pxz * (FPx + FUx) - Pyz * (FPy + FUy) - Pzz * (FPz + FUz) );
    return;

    /*
    // ANDERSON + SUCCESSIVE OVERRELAXATION
    U1_Hist.push_back( vec(U1x,U1y,U1z) );

    // Back up U(N,i)
    U1_i.push_back(vec(U1x, U1y, U1z));
    // Compute U(N,o)
    U1x = - Pxx * (FPx + FUx) - Pxy * (FPy + FUy) - Pxz * (FPz + FUz);
    U1y = - Pxy * (FPx + FUx) - Pyy * (FPy + FUy) - Pyz * (FPz + FUz);
    U1z = - Pxz * (FPx + FUx) - Pyz * (FPy + FUy) - Pzz * (FPz + FUz);
    // Back up U(N,o)
    U1_o.push_back(vec(U1x, U1y, U1z));

    // Compile Anderson in/out values from U1 record
    vec U1_i_mix = vec(0,0,0);
    vec U1_o_mix = vec(0,0,0);
    if (U1_i.size() > 2 && U1_o.size() > 2) {
        // Compute Anderson mixing factor from minimal rms deviation
        vec D_N   = U1_o[U1_o.size()-1] - U1_i[U1_i.size()-1];
        vec D_N_1 = U1_o[U1_o.size()-2] - U1_i[U1_i.size()-2];
        vec D_N_N_1 = D_N - D_N_1;
        double wAND = - D_N_1*D_N_N_1/(D_N_N_1*D_N_N_1);
        // Mingle new and old via wAND ...
        U1_i_mix = wAND*U1_i[U1_i.size()-1] + (1-wAND)*U1_i[U1_i.size()-2];
        U1_o_mix = wAND*U1_o[U1_o.size()-1] + (1-wAND)*U1_o[U1_o.size()-2];
    }
    else {
        U1_i_mix = U1_i.back();
        U1_o_mix = U1_o.back();
    }

    // Apply SOR
    vec U1 = (1-wSOR)*U1_i_mix + wSOR*U1_o_mix;
    U1x = U1.getX(); U1y = U1.getY(); U1z = U1.getZ();
    return;
    */
}

void APolarSite::InduceDirect() {
    U1_Hist.push_back( vec(0.,0.,0.) );
    U1x =  - Pxx * FPx - Pxy * FPy - Pxz * FPz;
    U1y =  - Pxy * FPx - Pyy * FPy - Pyz * FPz;
    U1z =  - Pxz * FPx - Pyz * FPy - Pzz * FPz;
    return;
}

double APolarSite::HistdU() {
    vec U0 = U1_Hist.back();
    vec U1 = vec(U1x, U1y, U1z);
    vec dU = U1 - U0;

    double abs_U0 = votca::tools::abs(U0);
    double abs_U1 = votca::tools::abs(U1);
    double abs_dU = votca::tools::abs(dU);
    double small = 1e-10; // in e*nm
    double abs_U = (abs_U1 > abs_U0 && small > abs_U0) ? abs_U1 : abs_U0;

    double dU_U = 1.;
    if (small > abs_U) dU_U = abs_U; // i.e.: U1, U0 << 1 enm
    else dU_U = abs_dU / abs_U;
    return dU_U;
}

double APolarSite::HistdU2() {
    vec U0 = U1_Hist.back();
    vec U1 = vec(U1x, U1y, U1z);
    double dU2 = (U1-U0)*(U1-U0);
    return dU2;
}

void APolarSite::Depolarize() {
    // Zero out induced moments
    U1x = U1y = U1z = 0.0;
    U1_Hist.clear();
    // Zero out fields
    FPx = FPy = FPz = 0.0;
    FUx = FUy = FUz = 0.0;
    return;
}

void APolarSite::PrintInfo(std::ostream &out) {

    printf("\n%2s %2d POS %+2.3f %+2.3f %+2.3f "
           "POL %+2.3f %+2.3f %+2.3f %+2.3f %+2.3f %+2.3f ",
            _name.c_str(), _id, _pos.getX(), _pos.getY(), _pos.getZ(),
            Pxx*1000, Pxy*1000, Pxz*1000, Pyy*1000, Pyz*1000, Pzz*1000);
}

void APolarSite::PrintTensorPDB(FILE *out, int state) {

    this->Charge(state);

    double scale = 30.;

    vec rxyz = _pos*scale;

    double rx = rxyz.getX();
    double ry = rxyz.getY();
    double rz = rxyz.getZ();

    vec end_pax = rxyz + 1.0*this->pax*this->eigenpxx/this->eigendamp;
    vec end_pay = rxyz + 1.0*this->pay*this->eigenpyy/this->eigendamp;
    vec end_paz = rxyz + 1.0*this->paz*this->eigenpzz/this->eigendamp;

    fprintf(out, "%-5s %+4.7f %+4.7f %+4.7f\n",
            this->_name.c_str(),
            rx,
            ry,
            rz);
    fprintf(out, "X     %+4.7f %+4.7f %+4.7f\n",
            end_pax.getX(),
            end_pax.getY(),
            end_pax.getZ());
    fprintf(out, "Y     %+4.7f %+4.7f %+4.7f\n",
            end_pay.getX(),
            end_pay.getY(),
            end_pay.getZ());
    fprintf(out, "Z     %+4.7f %+4.7f %+4.7f\n",
            end_paz.getX(),
            end_paz.getY(),
            end_paz.getZ());

}

void APolarSite::WritePdbLine(FILE *out, const string &tag) {

    fprintf(out, "ATOM  %5d %4s%1s%3s %1s%4d%1s   "
              "%8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s%4.7f\n",
         _id % 100000,          // Atom serial number           %5d
         _name.c_str(),         // Atom name                    %4s
         " ",                   // alternate location indicator.%1s
         tag.c_str(),           // Residue name.                %3s
         "A",                   // Chain identifier             %1s
         _id % 10000,           // Residue sequence number      %4d
         " ",                   // Insertion of residues.       %1s
         _pos.getX()*10,        // X in Angstroms               %8.3f
         _pos.getY()*10,        // Y in Angstroms               %8.3f
         _pos.getZ()*10,        // Z in Angstroms               %8.3f
         1.0,                   // Occupancy                    %6.2f
         0.0,                   // Temperature factor           %6.2f
         " ",                   // Segment identifier           %4s
         _name.c_str(),          // Element symbol               %2s
         " ",                   // Charge on the atom.          %2s
         Q00
         );
}

void APolarSite::WriteXyzLine(FILE *out, vec &shift, string format) {
/*
    double int2ext = 1.0;

    if (format == "gaussian") {
        int2ext = 10.;
    }
    else {
        int2ext = 10.;
    }
*/
    vec pos = _pos + shift;
    fprintf(out, "%-2s %+4.9f %+4.9f %+4.9f \n",
            _name.c_str(),
            pos.getX()*10, pos.getY()*10, pos.getZ()*10);
}

void APolarSite::WriteChkLine(FILE *out, vec &shift, bool split_dpl,
                             string format, double spacing) {

    vec pos = _pos + shift;

    string unit = "";

    if (format == "xyz") {
        unit = "angstrom";
    }
    else if (format == "gaussian") {
        unit = "angstrom";
    }

    // Take care of unit conversion
    double int2ext=0;

    if (unit == "nanometer") {
        int2ext = 1.;
    }
    else if (unit == "angstrom") {
        int2ext = 10.;
    }
    else if (unit == "bohr") {
        assert(false);
    }

    if (format == "xyz") {
        fprintf(out, "%2s ", _name.c_str());
    }

    // Print charge line
    fprintf(out, "%+4.9f %+4.9f %+4.9f %+4.7f \n",
            pos.getX()*int2ext,
            pos.getY()*int2ext,
            pos.getZ()*int2ext,
            Q00);

    // Split dipole moment onto charges (if desired)
    if (split_dpl) {

        vec tot_dpl = vec(U1x,U1y,U1z);

        if (_rank > 0) { tot_dpl += vec(Q1x,Q1y,Q1z); }

        matrix::eigensystem_t EIGEN;

        if (_rank == 2) {
            //cout << endl
            //     << "WARNING: Quadrupoles are not split onto point charges."
            //     << endl;

            int state = 0;

            double Qzz =      _Qs[state+1][4];
            double Qxx = -0.5*_Qs[state+1][4] + 0.5*sqrt(3)*_Qs[state+1][7];
            double Qyy = -0.5*_Qs[state+1][4] - 0.5*sqrt(3)*_Qs[state+1][7];

            double Qxy =  0.5*sqrt(3)*_Qs[state+1][8];
            double Qxz =  0.5*sqrt(3)*_Qs[state+1][5];
            double Qyz =  0.5*sqrt(3)*_Qs[state+1][6];

            matrix Q = matrix(vec(Qxx,Qxy,Qxz),
                              vec(Qxy,Qyy,Qyz),
                              vec(Qxz,Qyz,Qzz));

            Q.SolveEigensystem(EIGEN);
        }

        double a        = spacing;
        double mag_d    = abs(tot_dpl);
        vec    dir_d_0  = tot_dpl.normalize();
        vec    dir_d    = dir_d_0.normalize();
        vec    A        = pos + 0.5 * a * dir_d;
        vec    B        = pos - 0.5 * a * dir_d;
        double qA       = mag_d / a;
        double qB       = - qA;

        if (this->eigendamp == 0 || mag_d < 1e-9) {
            A = pos + 0.1*a*vec(1,0,0); // != pos since self-energy may diverge
            B = pos - 0.1*a*vec(1,0,0);
            qA = 0;
            qB = 0;
        }

        if (format == "xyz") {
            fprintf(out, " A ");
        }
        fprintf(out, "%+4.9f %+4.9f %+4.9f %+4.7f \n",
                A.getX()*int2ext,
                A.getY()*int2ext,
                A.getZ()*int2ext,
                qA);

        if (format == "xyz") {
            fprintf(out, " B ");
        }
        fprintf(out, "%+4.9f %+4.9f %+4.9f %+4.7f \n",
                B.getX()*int2ext,
                B.getY()*int2ext,
                B.getZ()*int2ext,
                qB);

        if (format == "xyz" && _rank == 2) {
            vec D1 = pos + 0.5 * a * EIGEN.eigenvecs[0];
            vec D2 = pos - 0.5 * a * EIGEN.eigenvecs[0];
            vec E1 = pos + 0.5 * a * EIGEN.eigenvecs[1];
            vec E2 = pos - 0.5 * a * EIGEN.eigenvecs[1];
            vec F1 = pos + 0.5 * a * EIGEN.eigenvecs[2];
            vec F2 = pos - 0.5 * a * EIGEN.eigenvecs[2];
            fprintf(out, " D %+4.9f %+4.9f %+4.9f \n",
                    D1.getX()*int2ext,
                    D1.getY()*int2ext,
                    D1.getZ()*int2ext);
            fprintf(out, " D %+4.9f %+4.9f %+4.9f \n",
                    D2.getX()*int2ext,
                    D2.getY()*int2ext,
                    D2.getZ()*int2ext);
            fprintf(out, " E %+4.9f %+4.9f %+4.9f \n",
                    E1.getX()*int2ext,
                    E1.getY()*int2ext,
                    E1.getZ()*int2ext);
            fprintf(out, " E %+4.9f %+4.9f %+4.9f \n",
                    E2.getX()*int2ext,
                    E2.getY()*int2ext,
                    E2.getZ()*int2ext);
            fprintf(out, " F %+4.9f %+4.9f %+4.9f \n",
                    F1.getX()*int2ext,
                    F1.getY()*int2ext,
                    F1.getZ()*int2ext);
            fprintf(out, " F %+4.9f %+4.9f %+4.9f \n",
                    F2.getX()*int2ext,
                    F2.getY()*int2ext,
                    F2.getZ()*int2ext);
        }

    }

}

void APolarSite::WriteMpsLine(std::ostream &out, string unit = "angstrom") {

    // Set conversion factor for higher-rank moments (e*nm**k to e*a0**k)
    double conv_dpl = 1./0.0529189379;
    double conv_qdr = conv_dpl*conv_dpl;
    // Set conversion factor for polarizabilities (nm**3 to A**3)
    double conv_pol = 1000;
    // Set conversion factor for positions (nm to ??)
    double conv_pos = 1.;
    if (unit == "angstrom") {
        conv_pos = 10.;
    }
    else if (unit == "nanometer") {
        conv_pos = 1.;
    }
    else assert(false); // Units error

    out << (boost::format(" %1$2s %2$+1.7f %3$+1.7f %4$+1.7f Rank %5$d\n")
            % _name % (_pos.getX()*conv_pos)
            % (_pos.getY()*conv_pos) % (_pos.getZ()*conv_pos)
            % _rank);
    // Charged
    out << (boost::format("    %1$+1.7f\n") % Q00);
    if (_rank > 0) {
        // Dipole z x y
        out << (boost::format("    %1$+1.7f %2$+1.7f %3$+1.7f\n")
            % (Q1z*conv_dpl) % (Q1x*conv_dpl) % (Q1y*conv_dpl));
        if (_rank > 1) {
            // Quadrupole 20 21c 21s 22c 22s
            out << (boost::format("    %1$+1.7f %2$+1.7f %3$+1.7f %4$+1.7f %5$+1.7f\n")
                % (Q20*conv_qdr) % (Q21c*conv_qdr) % (Q21s*conv_qdr)
                % (Q22c*conv_qdr) % (Q22s*conv_qdr));
        }
    }
    // Polarizability
    out << (boost::format("     P %1$+1.7f %2$+1.7f %3$+1.7f %4$+1.7f %5$+1.7f %6$+1.7f \n")
        % (Pxx*conv_pol) % (Pxy*conv_pol) % (Pxz*conv_pol)
        % (Pyy*conv_pol) % (Pyz*conv_pol) % (Pzz*conv_pol));

}

void APolarSite::WriteXmlLine(std::ostream &out) {
    out << "<aps>" << endl;
    out << _id << endl;
    out << _name << endl;
    out << ((_isVirtual) ? "virtual" : "notvirtual") << endl;
    out << "resolution" << (int)_resolution << endl;
    out << _pos << endl;
    out <<  _locX << " "; out <<  _locY << " "; out <<  _locZ << " "; out << endl;
    for (int state = -1; state < 2; ++state) {
        out << "<state>" << endl;
        out << state << endl;
        for (unsigned i = 0; i < _Qs[state+1].size(); ++i) {
            out << _Qs[state+1][i] << " ";
        }
        out << endl;
        out << "</state>" << endl;
    }
    out <<  _rank; out << endl;

    out << "<pol>" << endl;
    for (int state = -1; state < 2; ++state) {
        out << "<state>" << endl;
        out << _Ps[state+1].get(0,0) << " " << _Ps[state+1].get(0,1) << " " << _Ps[state+1].get(0,2) << endl;
        out << _Ps[state+1].get(1,0) << " " << _Ps[state+1].get(1,1) << " " << _Ps[state+1].get(1,2) << endl;
        out << _Ps[state+1].get(2,0) << " " << _Ps[state+1].get(2,1) << " " << _Ps[state+1].get(2,2) << endl;
        out << "</state>" << endl;
    }
    out << "</pol>" << endl;

    out << "<config>" << endl;
    out <<  Pxx << " "; out <<  Pxy << " "; out <<  Pxz << " "; out << endl;
    out <<  Pyy << " "; out <<  Pyz << " "; out << endl;
    out <<  Pzz << " "; out << endl;

    out <<  pax << " "; out <<  eigenpxx; out << endl;
    out <<  pay << " "; out <<  eigenpyy; out << endl;
    out <<  paz << " "; out <<  eigenpzz; out << endl;

    out <<  eigendamp; out << endl;

    out <<  Q00 << " "; out << endl;
    out <<  Q1x << " "; out <<  Q1y << " ";  out <<  Q1z << " "; out << endl;
    out <<  Q20 << " "; out <<  Q21c << " "; out <<  Q21s << " "; out <<  Q22c << " "; out <<  Q22s << " "; out << endl;
    out <<  Qxx << " "; out <<  Qxy << " ";  out <<  Qxz << " ";  out <<  Qyy << " ";  out <<  Qyz << " "; out <<  Qzz << " "; out << endl;

    out <<  U1x << " "; out <<  U1y << " "; out <<  U1z << " "; out << endl;
    out <<  FPx << " "; out <<  FPy << " "; out <<  FPz << " "; out << endl;
    out <<  FUx << " "; out <<  FUy << " "; out <<  FUz << " "; out << endl;
    out << "</config>" << endl;
    out << "</aps>" << endl;
    return;
}

vector<APolarSite*> APS_FROM_MPS(string filename, int state, QMThread *thread) {

    int poleCount = 1;
    double Q0_total = 0.0;
    string units = "";
    bool useDefaultPs = true;

    vector<APolarSite*> poles;
    APolarSite *thisPole = NULL;

    vector<double> Qs; // <- multipole moments
    matrix         P1; // <- dipole polarizability

    std::string line;
    std::ifstream intt;
    intt.open(filename.c_str());

    if (intt.is_open() ) {
        while ( intt.good() ) {

            std::getline(intt, line);
            vector<string> split;
            Tokenizer toker(line, " \t");
            toker.ToVector(split);

            if ( !split.size()      ||
                  split[0] == "!"   ||
                  split[0].substr(0,1) == "!" ) { continue; }

    // ! Interesting information here, e.g.
    // ! DCV2T opt
    // ! SP        RB3LYP          6-311+G(d,p)
    // Units bohr
    //
    // C          -4.2414603400   -3.8124751600    0.0017575736    Rank  2
    //  -0.3853409355
    //  -0.0002321905   0.2401559510   0.6602334308
    //  -0.7220625314   0.0004894995  -0.0003833545   0.4526409813  -0.50937399
    //  P 1.75

            // Units used
            if ( split[0] == "Units") {
                units = split[1];
                if (units != "bohr" && units != "angstrom") {
                    throw std::runtime_error( "Unit " + units + " in file "
                                            + filename + " not supported.");
                }
            }
            // element,  position,  rank limit
            else if ( split.size() == 6 ) {
                Qs.clear();
                int id = poleCount++;  // <- starts from 1
                string name = split[0];

                double BOHR2NM = 0.0529189379;
                double ANGSTROM2NM = 0.1;
                double x, y, z;

                if (units == "bohr") {
                    x = BOHR2NM * boost::lexical_cast<double>(split[1]);
                    y = BOHR2NM * boost::lexical_cast<double>(split[2]);
                    z = BOHR2NM * boost::lexical_cast<double>(split[3]);
                }
                else if (units == "angstrom") {
                    x = ANGSTROM2NM * boost::lexical_cast<double>(split[1]);
                    y = ANGSTROM2NM * boost::lexical_cast<double>(split[2]);
                    z = ANGSTROM2NM * boost::lexical_cast<double>(split[3]);
                }
                else {
                    throw std::runtime_error( "Unit " + units + " in file "
                                            + filename + " not supported.");
                }

                vec pos = vec(x,y,z);
                int rank = boost::lexical_cast<int>(split[5]);
                APolarSite *newPole = new APolarSite(id, name);
                newPole->setRank(rank);
                newPole->setPos(pos);
                poles.push_back(newPole);
                thisPole = newPole;
            }
            // 'P', dipole polarizability
            else if ( split[0] == "P") {
                double pxx, pxy, pxz;
                double      pyy, pyz;
                double           pzz;
                if (split.size() == 7) {
                    pxx = 1e-3 * boost::lexical_cast<double>(split[1]);
                    pxy = 1e-3 * boost::lexical_cast<double>(split[2]);
                    pxz = 1e-3 * boost::lexical_cast<double>(split[3]);
                    pyy = 1e-3 * boost::lexical_cast<double>(split[4]);
                    pyz = 1e-3 * boost::lexical_cast<double>(split[5]);
                    pzz = 1e-3 * boost::lexical_cast<double>(split[6]);
                    P1 = matrix(vec(pxx,pxy,pxz),
                                vec(pxy,pyy,pyz),
                                vec(pxz,pyz,pzz));
                }
                else if (split.size() == 2) {
                    pxx = 1e-3 * boost::lexical_cast<double>(split[1]);
                    pxy = 0.0;
                    pxz = 0.0;
                    pyy = pxx;
                    pyz = 0.0;
                    pzz = pxx;
                    P1 = matrix(vec(pxx,pxy,pxz),
                                vec(pxy,pyy,pyz),
                                vec(pxz,pyz,pzz));
                }
                else {
                    throw std::runtime_error("Invalid line in " + filename
                                             + ": " + line);
                }
                thisPole->setPs(P1, state);
                useDefaultPs = false;
            }
            // Multipole line
            else {
                int lineRank = int( sqrt(Qs.size()) + 0.5 );
                if (lineRank == 0) {
                    Q0_total += boost::lexical_cast<double>(split[0]);
                }
                for (unsigned i = 0; i < split.size(); i++) {
                    double qXYZ = boost::lexical_cast<double>(split[i]);
                    // Convert e*(a_0)^k to e*(nm)^k where k = rank
                    double BOHR2NM = 0.0529189379;
                    qXYZ *= pow(BOHR2NM, lineRank); // OVERRIDE
                    Qs.push_back(qXYZ);
                }
                thisPole->setQs(Qs, state);
            }
        } /* Exit loop over lines */
    }
    else { cout << endl << "ERROR: No such file " << filename << endl;
           throw runtime_error("Please supply input file.");           }

    if (thread == NULL)
    printf("\n... ... ... Reading %-25s -> N = %2d Q0(Sum) = %+1.7f ",
                          filename.c_str(), (int)poles.size(),  Q0_total);

    // Apply charge correction: Sum to closest integer
    int Q_integer = boost::math::iround(Q0_total);
    double dQ = (double(Q_integer) - Q0_total)/poles.size();

    double Q0_total_corr = 0.0;
    vector<APolarSite*>::iterator pit;
    for (pit = poles.begin(); pit < poles.end(); ++pit) {
        double Q_uncorr = (*pit)->getQs(state)[0];
        double Q_corr = Q_uncorr + dQ;
        (*pit)->setQ00(Q_corr, state);
        Q0_total_corr += (*pit)->getQs(state)[0];
    }

    if (thread == NULL)
    printf("=> dQ0 = %+1.1e, Q0(corr.) = %+1.0f",
            dQ, Q0_total_corr);

    if (useDefaultPs) {
        if (thread == NULL)
        cout << endl << "... ... ... NOTE Using default Thole polarizabilities "
             << "for charge state " << state << ". ";

        std::map<string,double> polar_table = POLAR_TABLE();

        vector< APolarSite* > ::iterator pol;
        for (pol = poles.begin(); pol < poles.end(); ++pol) {
            string elem = (*pol)->getName();
            double alpha = 0.0;

            try {
                alpha = polar_table.at(elem);
            }
            catch(out_of_range&){
                cout << endl << "WARNING No default polarizability given for "
                    << "polar site type '" << elem << "'. Defaulting to 1 A**3. "
                    << flush;
                alpha = 1e-3;
            }
//            // Original set of Thole polarizabilites
//            if      (elem == "C") { alpha = 1.75e-3;  } // <- conversion from
//            else if (elem == "H") { alpha = 0.696e-3; } //    A³ to nm³ = 10⁻³
//            else if (elem == "N") { alpha = 1.073e-3; }
//            else if (elem == "O") { alpha = 0.837e-3; }
//            else if (elem == "S") { alpha = 2.926e-3; }
//            // Different set of Thole polarizabilities
//            //if      (elem == "C") { alpha = 1.334e-3; } // <- conversion from
//            //else if (elem == "H") { alpha = 0.496e-3; } //    A³ to nm³ = 10⁻³
//            //else if (elem == "N") { alpha = 1.073e-3; }
//            //else if (elem == "O") { alpha = 0.837e-3; }
//            //else if (elem == "S") { alpha = 3.300e-3; }
//            else { throw runtime_error("No polarizability given "
//                                       "for polar site type " + elem + ". "); }

            P1 = matrix(vec(alpha,0,0),vec(0,alpha,0),vec(0,0,alpha));

            (*pol)->setPs(P1, state);
        }
    }

    vector< APolarSite* > ::iterator pol;
    for (pol = poles.begin(); pol < poles.end(); ++pol) {
        (*pol)->Charge(state);
    }

    return poles;
}

// All polarizabilities in nm**3
// Isotropic polarizability volume is evaluated from the tensor
// as (a_xx * a_yy * a_zz )^(1/3) for eigenvalues of the polarizability tensor
// GAUSSIAN: Bohr^3 to Angstrom^3 : k =  0.529177^3*10^(-3)
map<string,double> POLAR_TABLE() {
    map<string,double> polar_table;
    polar_table["H"] = 0.496e-3;
    polar_table["B"] = 2.797e-3;   // B3LYP/6-311+g(2d,2p)
    polar_table["C"] = 1.334e-3;
    polar_table["N"] = 1.073e-3;
    polar_table["O"] = 0.837e-3;
    polar_table["S"] = 2.926e-3;
    polar_table["F"] = 0.440e-3;
    polar_table["Si"] = 3.962e-3;   // B3LYP/6-311+g(2d,2p)
    polar_table["Cl"] = 1.504e-3;   // B3LYP/6-311+g(2d,2p)
    polar_table["Zn"] = 5.962e-3;   // B3LYP/6-311+g(2d,2p)
    polar_table["Cu"] = 6.295e-3;   // B3LYP/6-311+g(2d,2p)
    polar_table["Al"] = 5.80e-3;   //[1]P. Fuentealba, “The static dipole polarizability of aluminium atom: discrepancy between theory and experiment,” Chemical physics letters, vol. 397, no. 4, pp. 459–461, 2004.

    return polar_table;
}

}}

max_stars_repo_name: mbarbry/ctp

max_issues_repo_name: mbarbry/ctp

max_forks_repo_name: mbarbry/ctp


text:
#ifndef BNF_H
#define BNF_H

#include <boost/algorithm/string/trim.hpp>

#include "bnf_bootstrap_parser.hpp"
#include "../common/io.hpp"

namespace splicpp
{
	class bnf : public grammar
	{
	public:
		const stid L_RULE_NAME, L_ASSIGNMENT, L_EXPR_SEP, L_OPTIONAL, L_KLEENE, L_PLUS, L_NEWLINE, NL_SYNTAX;

		bnf()
		: grammar()

		, L_RULE_NAME(add_symbol(	new regex_literal("l_rule_name", "([a-zA-Z0-9\\-\\_]+)")))
		, L_ASSIGNMENT(add_symbol(	new static_literal("l_assignment", ":==")))
		, L_EXPR_SEP(add_symbol(	new static_literal("l_expr_sep", "|")))
		, L_OPTIONAL(add_symbol(	new static_literal("l_optional", "?")))
		, L_KLEENE(add_symbol(		new static_literal("l_kleene", "*")))
		, L_PLUS(add_symbol(		new static_literal("l_plus", "+")))
		, L_NEWLINE(add_symbol(		new static_literal("l_newline", "\n")))

		, NL_SYNTAX(add_symbol(		new non_literal("nl_syntax")))
		{
			ignore(' ');
			ignore('\t');

			splicpp::bnf_bootstrap_parser p;
			std::string lang = readfile("grammars/bnf.bnf");
			boost::algorithm::trim(lang);

			p.parse((*this), lang);

			add_rule(rule(NL_START) + NL_SYNTAX);
		}
	};
}

#endif

max_stars_repo_name: Wassasin/splicpp

max_issues_repo_name: Wassasin/splicpp

max_forks_repo_name: Wassasin/splicpp


text:
/*
 * Copyright (c) 2014 Codilime.
 */

#include "ifmap/ifmap_config_listener.h"

#include <boost/assign/list_of.hpp>

#include "ifmap/ifmap_dependency_tracker.h"

IFMapConfigListener::ConfigDelta::ConfigDelta() {
}

IFMapConfigListener::ConfigDelta::ConfigDelta(
  const IFMapConfigListener::ConfigDelta &rhs)
    : id_type(rhs.id_type), id_name(rhs.id_name),
      node(rhs.node), obj(rhs.obj) {
}

IFMapConfigListener::IFMapConfigListener(ConfigManager *manager,
                                          const char *concurrency)
    : manager_(manager), kConcurrency_(concurrency) {}

IFMapConfigListener::~IFMapConfigListener() {}

//
// Initialize the IFMapConfigListener.
//
// Create and initialize the DependencyTracker.
// We register one listener for the IFMapLinkTable and a listener for each of
// the relevant IFMapTables.
//
// Pure virtual DependencyTrackerInit() should initialize dependency tracker's
// policy_map, which is then used to register listeners for specified tables.
//
void IFMapConfigListener::Initialize() {
    DB *database = manager_->database();

    tracker_.reset(
        new IFMapDependencyTracker(
            database, manager_->graph(),
            boost::bind(&IFMapConfigListener::ChangeListAdd, this, _1)));
    DependencyTrackerInit();

    DBTable *link_table = static_cast<DBTable *>(
        database->FindTable("__ifmap_metadata__.0"));
    assert(link_table != NULL);

    DBTable::ListenerId id = link_table->Register(
        boost::bind(&IFMapConfigListener::LinkObserver, this, _1, _2));
    table_map_.insert(make_pair(link_table->name(), id));

    BOOST_FOREACH(
      const IFMapDependencyTracker::NodeEventPolicy::value_type &policy,
        *tracker_->policy_map()) {
        const char *schema_typename= policy.first.c_str();
        IFMapTable *table = IFMapTable::FindTable(database, schema_typename);
        assert(table);
        DBTable::ListenerId id = table->Register(
                boost::bind(&IFMapConfigListener::NodeObserver, this, _1, _2));
        table_map_.insert(make_pair(table->name(), id));
    }
}

//
// Unregister listeners for all the IFMapTables.
//
void IFMapConfigListener::Terminate() {
    DB *database = manager_->database();

    for (TableMap::iterator iter = table_map_.begin();
         iter != table_map_.end(); ++iter) {
        IFMapTable *table =
            static_cast<IFMapTable *>(database->FindTable(iter->first));
        assert(table);
        table->Unregister(iter->second);
    }
    table_map_.clear();
}

//
// Get the the DB in the ConfigManager.
//
DB *IFMapConfigListener::database() {
    return manager_->database();
}

//
// Ask the DependencyTracker to build up the ChangeList.
//
void IFMapConfigListener::GetChangeList(ChangeList *change_list) {
    CHECK_CONCURRENCY(kConcurrency_.c_str());

    tracker_->PropagateChanges();
    tracker_->Clear();
    change_list->swap(change_list_);
}

IFMapDependencyTracker *IFMapConfigListener::get_dependency_tracker() {
    return tracker_.get();
}

//
// Add an IFMapNode to the ChangeList by creating a ConfigDelta.
//
// We take references on the IFMapNode and the IFMapObject.
//
void IFMapConfigListener::ChangeListAdd(IFMapNode *node) {
    CHECK_CONCURRENCY(kConcurrency_.c_str(), "db::DBTable");

    IFMapTable *table = node->table();
    TableMap::const_iterator tid = table_map_.find(table->name());
    if (tid == table_map_.end()) {
        return;
    }

    ConfigDelta delta;
    delta.id_type = table->Typename();
    delta.id_name = node->name();
    if (!node->IsDeleted()) {
        DBState *current = node->GetState(table, tid->second);
        if (current == NULL) {
            delta.node = IFMapNodeRef(new IFMapNodeProxy(node, tid->second));
        }
        delta.obj = IFMapObjectRef(node->GetObject());
    }
    change_list_.push_back(delta);
}

//
// Callback to handle node events in the IFMapTables.
// Adds the node itself to the ChangeList and informs the DependencyTracker
// about the node event.
//
void IFMapConfigListener::NodeObserver(
    DBTablePartBase *root, DBEntryBase *db_entry) {
    CHECK_CONCURRENCY("db::DBTable");

    // Ignore deleted nodes for which the configuration code doesn't hold
    // state. This is the case with ie. BgpRouter objects other than the local
    // node.
    IFMapNode *node = static_cast<IFMapNode *>(db_entry);
    if (node->IsDeleted()) {
        IFMapTable *table = node->table();
        TableMap::const_iterator tid = table_map_.find(table->name());
        assert(tid != table_map_.end());
        if (node->GetState(table, tid->second) == NULL) {
            return;
        }
    }

    tracker_->NodeEvent(node);
    manager_->OnChange();
}

//
// Callback to handle link events in the IFMapLinkTable.
// Informs the DependencyTracker about the link event.
//
void IFMapConfigListener::LinkObserver(
    DBTablePartBase *root, DBEntryBase *db_entry) {
    CHECK_CONCURRENCY("db::DBTable");

    IFMapLink *link = static_cast<IFMapLink *>(db_entry);
    IFMapNode *left = link->LeftNode(database());
    IFMapNode *right = link->RightNode(database());
    if (tracker_->LinkEvent(link->metadata(), left, right)) {
        manager_->OnChange();
    }
}

max_stars_repo_name: Juniper/contrail-dev-controller

max_issues_repo_name: Juniper/contrail-dev-controller

max_forks_repo_name: Juniper/contrail-dev-controller


text:
//==============================================================================
//         Copyright 2003 - 2011   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2011   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_DISPATCH_FUNCTOR_META_CALL_HPP_INCLUDED
#define BOOST_DISPATCH_FUNCTOR_META_CALL_HPP_INCLUDED

/*!
 * \file
 * \brief Defines and implements the \ref boost::dispatch::meta::call \metafunction
 */

#include <boost/config.hpp>
#include <boost/dispatch/functor/details/dispatch.hpp>
#include <boost/dispatch/meta/result_of.hpp>
#include <boost/dispatch/meta/as_ref.hpp>

namespace boost { namespace dispatch { namespace meta
{
  //============================================================================
  /*!
   * \ref call is a \metafunction that, for a given \c Tag, \c EvalContext and
   * argument lists \c ...Args, computes the return type of
   * functor<Tag,EvalContext>::operator()(Args...).
   *
   * \param Sig         function type build from Tag(Args...)
   * \param EvalContext evaluation context to use for type deduction
   *
   * \par Models:
   * \metafunction
   *
   * \par Semantics:
   * For a given \c Tag, \c EvalContext and argument lists \c ...Args,
   *
   * \code
   * typedef call<Tag(Args...),EvalContext>::type type;
   * \endcode
   *
   * is equivalent to:
   * \code
   * typedef decltype(functor<Tag,EvalContext>(declval<Args...>())) type;
   * \endcode
   */
  //============================================================================
  template< class Sig
          , class EvalContext = void
          > struct call {};

  template<class Tag, class... Args, class Site>
  struct  call<Tag(Args...),Site>
        : meta::result_of<typename meta::dispatch_call<Tag(typename meta::as_ref<Args>::type...), Site>::type(Args...)>
  {
  };
} } }

#endif


text:
#include <iostream>
#include <string>
#include <vector>

#include <boost/program_options.hpp>

#include "query.h"
#include "queryplan.h"
#include "serde.h"
#include "options.h"

class CutOptions : public Options {
public:
    CutOptions() {
        namespace po = boost::program_options;

        description_.add_options()
            ("field,f", po::value(&fields_)->composing(), "Include this field in output.")
        ;
        add_positional_argument("field", 1, -1);
    }

    std::vector<std::string> get_fields() const { return fields_; }

private:
    std::vector<std::string> fields_;
};

int main(int argc, char **argv) {
    CutOptions options;
    if (!options.parse(argc, argv)) {
        return 1;
    }

    auto query_plan = load_query_plan(std::cin);
    if (!query_plan) {
        std::cerr << "Unable to parse query plan from standard input." << std::endl;
        return 1;
    }

    query_plan->select.emplace(
        query_plan->select->get_tablename(),
        options.get_fields()
    );

    dump_query_plan(*query_plan, std::cout);
    return 0;
}

max_stars_repo_name: StefansM/dataset-coreutils

max_issues_repo_name: StefansM/dataset-coreutils

max_forks_repo_name: StefansM/dataset-coreutils


text:
/*
 * MIT License
 *
 * Copyright (c) 2018 Qwant Research
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 *
 * Author:
 * Noel Martin (n.martin@qwantresearch.com)
 *
 */

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <thread>
#include <atomic>

#include <unistd.h>
#include <curl/curl.h>

#ifdef MMZ_PROFILE
#include <gperftools/heap-profiler.h>
#endif

#include <boost/program_options.hpp>
namespace po = boost::program_options;

#include "tsafe/thread_safe_queue.h"

#include "common/common.hpp"
#include "spider/spider.hpp"
#include "urlserver/urlserver.hpp"

using namespace mermoz;

int main (int argc, char** argv)
{
# ifdef MMZ_PROFILE
  HeapProfilerStart("mmz");
# endif

  /*
   * Defining options for code
   */
  po::options_description desc("Allowed options");
  desc.add_options()
  ("help", "displays this message")
  ("settings", po::value<std::string>(), "setting file for the whole Mermoz run")
  ("seeds", po::value<std::string>(), "list of seeds for the current run")
  ;

  po::variables_map vmap;
  po::store(po::parse_command_line(argc, argv, desc), vmap);
  po::notify(vmap);

  /*
   * If one asks '--help'
   */
  if (vmap.count("help")) {
    std::cout << desc << std::endl;
    return 1;
  }

  /*
   * Parsing the 'settings' file
   */
  size_t pos;
  std::string line;
  std::ifstream settingsfile(vmap["settings"].as<std::string>());

  std::string user_agent;
  unsigned int nfetchers {0};
  unsigned int nparsers {0};
  int max_ram {0};

  while(!settingsfile.eof()) {
    line.clear();
    std::getline(settingsfile, line);

    unsigned char c;
    while ((c = *(line.end()-1)) < 0x20 && !line.empty())
      line.pop_back();

    if ((pos = line.find("fetchers")) != std::string::npos)
      nfetchers = static_cast<unsigned int>(std::atoi(line.substr(pos + 9).c_str()));
    else if ((pos = line.find("parsers")) != std::string::npos)
      nparsers = static_cast<unsigned int>(std::atoi(line.substr(pos + 8).c_str()));
    else if ((pos = line.find("user-agent")) != std::string::npos)
      user_agent = line.substr(pos + 11);
    else if ((pos = line.find("max-ram")) != std::string::npos)
      max_ram = std::atoi(line.substr(pos + 8).c_str());
  }
  settingsfile.close();

  /*
   * One verifies that all the mandatory settings
   * where included
   */
  if (user_agent.empty() ||
      nfetchers == 0 ||
      nparsers == 0 ||
      max_ram == 0) {
    print_error("Wrong settings Mermoz cannot start");
  } else {
    print_strong_log("Staring of Mermoz with following settings:");

    std::ostringstream oss;

    oss << "Fetchers: " << nfetchers;
    print_strong_log(oss.str());

    oss.str("");
    oss << "Parsers: " << nparsers;
    print_strong_log(oss.str());

    oss.str("");
    oss << "User-agent: " << user_agent;
    print_strong_log(oss.str());

    oss.str("");
    oss << "Max-ram (GB): " << max_ram;
    print_strong_log(oss.str());
  }

  bool status = true;

  TSQueueVector url_queues(nfetchers);
  TSQueueVector content_queues(nparsers);
  MemSec mem_sec(max_ram * MemSec::GB);

  unsigned int queue_id {0};
  std::string link;
  std::ifstream seedfile(vmap["seeds"].as<std::string>());
  while (!seedfile.eof()) {
    link = {""};
    seedfile >> link;
    if (!link.empty()) {
      urlfactory::UrlParser up(link);
      std::string message;
      std::string host {up.get_host()};
      pack(message, {&host, &link});
      url_queues[queue_id].push(message);
      mem_sec += message.size();
      queue_id++;

      if (queue_id >= url_queues.size()) {
        queue_id = 0;
      }
    }
  }
  seedfile.close();

  /* Must initialize libcurl before any threads are started */
  curl_global_init(CURL_GLOBAL_ALL);

  /*
   * Settings for the UrlServer
   */
  UrlServerSettings uset = {
    user_agent,
    &mem_sec
  };

  std::thread userv(urlserver,
                    &status,
                    &uset,
                    &content_queues,
                    &url_queues);

  std::atomic<uint64_t> nfetched;
  nfetched = 0;

  std::atomic<uint64_t> nparsed;
  nparsed = 0;

  /*
   * Settings for the Spider
   */
  SpiderSettings sset = {
    nfetchers,
    nparsers,
    user_agent,
    &nfetched,
    &nparsed,
    &mem_sec,
  };

  std::thread spdr(spider,
                   &status,
                   &sset,
                   &url_queues,
                   &content_queues);

  std::ofstream ofp("log.out");

  ofp << "# time urls contents fetched parsed mem(MB)" << std::endl;

  while (status) {
    sleep(10);
    std::time_t t = std::time(nullptr);
    std::tm tm = *std::localtime(&t);

    unsigned long wait_url {0};
    for (auto& queue : url_queues) {
      wait_url += queue.size();
    }

    unsigned long wait_content {0};
    for (auto& queue : content_queues) {
      wait_content += queue.size();
    }

    ofp << tm.tm_hour*3600 + tm.tm_min*60 + tm.tm_sec << " ";
    ofp << wait_url << " ";
    ofp << wait_content << " ";
    ofp << nfetched << " ";
    ofp << nparsed << " ";
    ofp << mem_sec.get_mem()/(1UL << 20) << std::endl;
  }

# ifdef MMZ_PROFILE
  HeapProfilerStop();
# endif

  userv.join();
  spdr.join();

  return 0;
}

max_stars_repo_name: stephaneerard/mermoz

max_issues_repo_name: stephaneerard/mermoz

max_forks_repo_name: stephaneerard/mermoz


text:
#include <stdio.h>
#include <string.h>
#include <string>
#include <algorithm>
#include <iostream>
#include <functional>
#include <boost/any.hpp>
#include "traversal/P.h"
#include "structure/Vertex.h"
#include "traversal/GraphTraversal.h"
#include "structure/CPUGraph.h"

std::vector<std::function<void()>> tests;

int main(int argc, char* argv[]) {
	tests.push_back([]() {
		CPUGraph graph;
		auto explanation = graph.traversal()->addV()->has("a", "b")->explain();
		printf("%s\n", explanation.c_str());

		graph.traversal()->addV()->addV()->addV()->iterate();
		printf("Added 3 vertices to the graph.\n");
		graph.traversal()->V()->property("a", std::string("b"))->iterate();

		std::vector<Vertex*> vertices;
		std::for_each(graph.vertices().begin(), graph.vertices().end(), [&vertices](Vertex* v) {
			vertices.push_back(v);
		});

		VertexProperty<boost::any>* p = vertices.front()->property("a");
		std::cout << boost::any_cast<uint64_t>(vertices.front()->id()) << " Set property " << p->key() << " to " << boost::any_cast<std::string>(p->value()) << "\n";

		p = vertices[1]->property("a");
		std::cout << boost::any_cast<uint64_t>(vertices[1]->id()) << " Set property " << p->key() << " to " << boost::any_cast<std::string>(p->value()) << "\n";

		p = vertices[2]->property("a");
		std::cout << boost::any_cast<uint64_t>(vertices[2]->id()) << " Set property " << p->key() << " to " << boost::any_cast<std::string>(p->value()) << "\n";

		std::cout << "size: " << graph.vertices().size() << "\n";

		graph.traversal()->V()->has("a", std::string("b"))->forEachRemaining([](boost::any& v){
			Vertex* w = boost::any_cast<Vertex*>(v);
			std::cout << "The Vertex with id " << boost::any_cast<uint64_t>(w->id()) << " has property a = b\n";
		});

		/*
		graph.traversal()->addE()
			->from(__->V()->hasId(1))
			->to(__->V()->hasId(1))
			->V().hasId(1)
			->addE()
			->to(__->V()->hasId(1))
			->iterate();
		*/
		graph.traversal()->addE("basic_edge")->from(vertices[0])->to(vertices[1])->iterate();
		cout << "Edge successfully added!\n";

		std::cout << graph.vertices().size() << "\n";
		for_each(vertices.begin(), vertices.end(), [](Vertex* v){ std::cout << boost::any_cast<uint64_t>(v->id()) << " "; });
		std::cout << "\n---------\n";

		std::vector<BitEdge*> edges = ((BitVertex*)graph.vertices().front())->edges(OUT);
		cout << "# edges: " << edges.size() << "\n";
		for_each(edges.begin(), edges.end(), [](BitEdge* edg) {
			uint64_t edge_id = boost::any_cast<uint64_t>(edg->id());
			uint64_t v1_id = boost::any_cast<uint64_t>(edg->outV()->id());
			uint64_t v2_id = boost::any_cast<uint64_t>(edg->inV()->id());
			std::cout << edge_id << ": " << v1_id << "->" << v2_id << "\n";
		});

	});

	for_each(tests.begin(), tests.end(), [](auto n){
		try {
			n();
		} catch(const std::exception & err) { cout << err.what() << "\n"; return -1; }
	});
}

max_stars_repo_name: dylanht/BitGraph

max_issues_repo_name: dylanht/BitGraph

max_forks_repo_name: dylanht/BitGraph


text:
/**
 *  @file fc/cmt/asio.hpp
 *  @brief defines wrappers for boost::asio functions
 */
#pragma once
#include <boost/asio.hpp>
#include <boost/bind.hpp>
#include <fc/thread/future.hpp>
#include <fc/io/iostream.hpp>

namespace fc {
/**
 *  @brief defines fc wrappers for boost::asio functions.
 */
namespace asio {
    /**
     *  @brief internal implementation types/methods for fc::asio
     */
    namespace detail {
        using namespace fc;

        class read_write_handler
        {
        public:
          read_write_handler(const promise<size_t>::ptr& p);
          void operator()(const boost::syCC::error_code& ec, size_t bytes_transferred);
        private:
          promise<size_t>::ptr _completion_promise;
        };

        class read_write_handler_with_buffer
        {
        public:
          read_write_handler_with_buffer(const promise<size_t>::ptr& p,
                                         const std::shared_ptr<const char>& buffer);
          void operator()(const boost::syCC::error_code& ec, size_t bytes_transferred);
        private:
          promise<size_t>::ptr _completion_promise;
          std::shared_ptr<const char> _buffer;
        };

        //void read_write_handler( const promise<size_t>::ptr& p,
        //                         const boost::syCC::error_code& ec,
        //                         size_t bytes_transferred );
        void read_write_handler_ec( promise<size_t>* p,
                                    boost::syCC::error_code* oec,
                                    const boost::syCC::error_code& ec,
                                    size_t bytes_transferred );
        void error_handler( const promise<void>::ptr& p,
                              const boost::syCC::error_code& ec );
        void error_handler_ec( promise<boost::syCC::error_code>* p,
                              const boost::syCC::error_code& ec );

        template<typename C>
        struct non_blocking {
          bool operator()( C& c ) { return c.non_blocking(); }
          bool operator()( C& c, bool s ) { c.non_blocking(s); return true; }
        };

        #if WIN32  // windows stream handles do not support non blocking!
	       template<>
         struct non_blocking<boost::asio::windows::stream_handle> {
	          typedef boost::asio::windows::stream_handle C;
            bool operator()( C& ) { return false; }
            bool operator()( C&, bool ) { return false; }
        };
        #endif
    }
    /**
     * @return the default boost::asio::io_service for use with fc::asio
     *
     * This IO service is automatically running in its own thread to service asynchronous
     * requests without blocking any other threads.
     */
    boost::asio::io_service& default_io_service(bool cleanup = false);

    /**
     *  @brief wraps boost::asio::async_read
     *  @pre s.non_blocking() == true
     *  @return the number of bytes read.
     */
    template<typename AsyncReadStream, typename MutableBufferSequence>
    size_t read( AsyncReadStream& s, const MutableBufferSequence& buf ) {
        promise<size_t>::ptr p(new promise<size_t>("fc::asio::read"));
        boost::asio::async_read( s, buf, detail::read_write_handler(p) );
        return p->wait();
    }
    /**
     *  This method will read at least 1 byte from the stream and will
     *  cooperatively block until that byte is available or an error occurs.
     *
     *  If the stream is not in 'non-blocking' mode it will be put in 'non-blocking'
     *  mode it the stream supports s.non_blocking() and s.non_blocking(bool).
     *
     *  If in non blocking mode, the call will be synchronous avoiding heap allocs
     *  and context switching. If the sync call returns 'would block' then an
     *  promise is created and an async read is generated.
     *
     *  @return the number of bytes read.
     */
    template<typename AsyncReadStream, typename MutableBufferSequence>
    future<size_t> read_some(AsyncReadStream& s, const MutableBufferSequence& buf)
    {
      promise<size_t>::ptr completion_promise(new promise<size_t>("fc::asio::async_read_some"));
      s.async_read_some(buf, detail::read_write_handler(completion_promise));
      return completion_promise;//->wait();
    }

    template<typename AsyncReadStream>
    future<size_t> read_some(AsyncReadStream& s, char* buffer, size_t length, size_t offset = 0)
    {
      promise<size_t>::ptr completion_promise(new promise<size_t>("fc::asio::async_read_some"));
      s.async_read_some(boost::asio::buffer(buffer + offset, length),
                        detail::read_write_handler(completion_promise));
      return completion_promise;//->wait();
    }

    template<typename AsyncReadStream>
    future<size_t> read_some(AsyncReadStream& s, const std::shared_ptr<char>& buffer, size_t length, size_t offset)
    {
      promise<size_t>::ptr completion_promise(new promise<size_t>("fc::asio::async_read_some"));
      s.async_read_some(boost::asio::buffer(buffer.get() + offset, length),
                        detail::read_write_handler_with_buffer(completion_promise, buffer));
      return completion_promise;//->wait();
    }

    template<typename AsyncReadStream, typename MutableBufferSequence>
    void async_read_some(AsyncReadStream& s, const MutableBufferSequence& buf, promise<size_t>::ptr completion_promise)
    {
      s.async_read_some(buf, detail::read_write_handler(completion_promise));
    }

    template<typename AsyncReadStream>
    void async_read_some(AsyncReadStream& s, char* buffer,
                         size_t length, promise<size_t>::ptr completion_promise)
    {
      s.async_read_some(boost::asio::buffer(buffer, length), detail::read_write_handler(completion_promise));
    }

    template<typename AsyncReadStream>
    void async_read_some(AsyncReadStream& s, const std::shared_ptr<char>& buffer,
                         size_t length, size_t offset, promise<size_t>::ptr completion_promise)
    {
      s.async_read_some(boost::asio::buffer(buffer.get() + offset, length), detail::read_write_handler_with_buffer(completion_promise, buffer));
    }

    template<typename AsyncReadStream>
    size_t read_some( AsyncReadStream& s, boost::asio::streambuf& buf )
    {
        char buffer[1024];
        size_t bytes_read = read_some( s, boost::asio::buffer( buffer, sizeof(buffer) ) );
        buf.sputn( buffer, bytes_read );
        return bytes_read;
    }

    /** @brief wraps boost::asio::async_write
     *  @return the number of bytes written
     */
    template<typename AsyncWriteStream, typename ConstBufferSequence>
    size_t write( AsyncWriteStream& s, const ConstBufferSequence& buf ) {
        promise<size_t>::ptr p(new promise<size_t>("fc::asio::write"));
        boost::asio::async_write(s, buf, detail::read_write_handler(p));
        return p->wait();
    }

    /**
     *  @pre s.non_blocking() == true
     *  @brief wraps boost::asio::async_write_some
     *  @return the number of bytes written
     */
    template<typename AsyncWriteStream, typename ConstBufferSequence>
    future<size_t> write_some( AsyncWriteStream& s, const ConstBufferSequence& buf ) {
        promise<size_t>::ptr p(new promise<size_t>("fc::asio::write_some"));
        s.async_write_some( buf, detail::read_write_handler(p));
        return p; //->wait();
    }

    template<typename AsyncWriteStream>
    future<size_t> write_some( AsyncWriteStream& s, const char* buffer,
                               size_t length, size_t offset = 0) {
        promise<size_t>::ptr p(new promise<size_t>("fc::asio::write_some"));
        s.async_write_some( boost::asio::buffer(buffer + offset, length), detail::read_write_handler(p));
        return p; //->wait();
    }

    template<typename AsyncWriteStream>
    future<size_t> write_some( AsyncWriteStream& s, const std::shared_ptr<const char>& buffer,
                               size_t length, size_t offset ) {
        promise<size_t>::ptr p(new promise<size_t>("fc::asio::write_some"));
        s.async_write_some( boost::asio::buffer(buffer.get() + offset, length), detail::read_write_handler_with_buffer(p, buffer));
        return p; //->wait();
    }

    /**
    *  @pre s.non_blocking() == true
    *  @brief wraps boost::asio::async_write_some
    *  @return the number of bytes written
    */
    template<typename AsyncWriteStream, typename ConstBufferSequence>
    void async_write_some(AsyncWriteStream& s, const ConstBufferSequence& buf, promise<size_t>::ptr completion_promise) {
      s.async_write_some(buf, detail::read_write_handler(completion_promise));
    }

    template<typename AsyncWriteStream>
    void async_write_some(AsyncWriteStream& s, const char* buffer,
                          size_t length, promise<size_t>::ptr completion_promise) {
      s.async_write_some(boost::asio::buffer(buffer, length),
                         detail::read_write_handler(completion_promise));
    }

    template<typename AsyncWriteStream>
    void async_write_some(AsyncWriteStream& s, const std::shared_ptr<const char>& buffer,
                          size_t length, size_t offset, promise<size_t>::ptr completion_promise) {
      s.async_write_some(boost::asio::buffer(buffer.get() + offset, length),
                         detail::read_write_handler_with_buffer(completion_promise, buffer));
    }

    namespace tcp {
        typedef boost::asio::ip::tcp::endpoint endpoint;
        typedef boost::asio::ip::tcp::resolver::iterator resolver_iterator;
        typedef boost::asio::ip::tcp::resolver resolver;
        std::vector<endpoint> resolve( const std::string& hostname, const std::string& port );

        /** @brief wraps boost::asio::async_accept
          * @post sock is connected
          * @post sock.non_blocking() == true
          * @throw on error.
          */
        template<typename SocketType, typename AcceptorType>
        void accept( AcceptorType& acc, SocketType& sock ) {
            //promise<boost::syCC::error_code>::ptr p( new promise<boost::syCC::error_code>("fc::asio::tcp::accept") );
            promise<void>::ptr p( new promise<void>("fc::asio::tcp::accept") );
            acc.async_accept( sock, boost::bind( fc::asio::detail::error_handler, p, _1 ) );
            p->wait();
            //if( ec ) BOOST_THROW_EXCEPTION( boost::syCC::syCC_error(ec) );
        }

        /** @brief wraps boost::asio::socket::async_connect
          * @post sock.non_blocking() == true
          * @throw on error
          */
        template<typename AsyncSocket, typename EndpointType>
        void connect( AsyncSocket& sock, const EndpointType& ep ) {
            promise<void>::ptr p(new promise<void>("fc::asio::tcp::connect"));
            sock.async_connect( ep, boost::bind( fc::asio::detail::error_handler, p, _1 ) );
            p->wait();
            //if( ec ) BOOST_THROW_EXCEPTION( boost::syCC::syCC_error(ec) );
        }
    }
    namespace udp {
        typedef boost::asio::ip::udp::endpoint endpoint;
        typedef boost::asio::ip::udp::resolver::iterator resolver_iterator;
        typedef boost::asio::ip::udp::resolver resolver;
        /// @brief resolve all udp::endpoints for hostname:port
        std::vector<endpoint> resolve( resolver& r, const std::string& hostname,
                                                         const std::string& port );
    }

    template<typename AsyncReadStream>
    class istream : public virtual fc::istream
    {
       public:
          istream( std::shared_ptr<AsyncReadStream> str )
          :_stream( fc::move(str) ){}

          virtual size_t readsome( char* buf, size_t len )
          {
             return fc::asio::read_some(*_stream, buf, len).wait();
          }
          virtual size_t readsome( const std::shared_ptr<char>& buf, size_t len, size_t offset )
          {
             return fc::asio::read_some(*_stream, buf, len, offset).wait();
          }

       private:
          std::shared_ptr<AsyncReadStream> _stream;
    };

    template<typename AsyncWriteStream>
    class ostream : public virtual fc::ostream
    {
       public:
          ostream( std::shared_ptr<AsyncWriteStream> str )
          :_stream( fc::move(str) ){}

          virtual size_t writesome( const char* buf, size_t len )
          {
             return fc::asio::write_some(*_stream, buf, len).wait();
          }

          virtual size_t     writesome( const std::shared_ptr<const char>& buf, size_t len, size_t offset )
          {
             return fc::asio::write_some(*_stream, buf, len, offset).wait();
          }

          virtual void       close(){ _stream->close(); }
          virtual void       flush() {}
       private:
          std::shared_ptr<AsyncWriteStream> _stream;
    };

} } // namespace fc::asio

max_stars_repo_name: DunnCreativeSS/ccd

max_issues_repo_name: DunnCreativeSS/ccd

max_forks_repo_name: DunnCreativeSS/ccd