text:
#pragma once
#ifndef INC_SERVER_H
#define INC_SERVER_H

#include <string>
#include <list>
#include <vector>
#include <queue>
#include <memory>
#include <boost/unordered_map.hpp>
#include <boost/regex.hpp>
#include <boost/atomic.hpp>
#include <boost/thread/mutex.hpp>
#include <boost/thread/lock_guard.hpp>
#include <boost/lockfree/spsc_queue.hpp>

#include "CSingleton.hpp"

using std::string;
using std::list;
using std::vector;
using std::queue;
using std::shared_ptr;
using boost::atomic;
using boost::unordered_map;
using boost::mutex;

class CCallback;
typedef shared_ptr<CCallback> Callback_t;

struct Channel
{
	typedef unsigned int Id_t;
	static const Id_t Invalid = -1;

	enum class Types
	{
		INVALID,
		TEMPORARY,
		SEMI_PERMANENT,
		PERMANENT
	};

	enum class QueryData
	{
		INVALID,
		CHANNEL_TOPIC, //channel_topic [string]
		CHANNEL_DESCRIPTION, //channel_description [string]
		CHANNEL_CODEC, //channel_codec [int]
		CHANNEL_CODEC_QUALITY, //channel_codec_quality [int]
		CHANNEL_FORCED_SILENCE, //channel_forced_silence [int]
		CHANNEL_ICON_ID, //channel_icon_id [int]
		CHANNEL_CODEC_IS_UNENCRYPTED, //channel_codec_is_unencrypted [int]
		CHANNEL_SECONDS_EMPTY //seconds_empty [int]
	};

	Id_t
		ParentId = Invalid,
		OrderId = Invalid;

	string Name;
	Types Type = Types::INVALID;
	bool
		HasPassword = false,
		WasPasswordToggled = false;
	int RequiredTalkPower = 0;

	list<unsigned int> Clients;
	int MaxClients = -1;
};
typedef shared_ptr<Channel> Channel_t;

struct Client
{
	typedef unsigned int Id_t;
	static const Id_t Invalid = 0;

	enum class KickTypes
	{
		INVALID,
		CHANNEL,
		SERVER
	};

	enum class QueryData
	{
		INVALID,
		CLIENT_NICKNAME, //client_nickname [string]
		CLIENT_VERSION, //client_version [string]
		CLIENT_PLATFORM, //client_platform [string]
		CLIENT_INPUT_MUTED, //client_input_muted [int]
		CLIENT_OUTPUT_MUTED, //client_output_muted [int]
		CLIENT_INPUT_HARDWARE, //client_input_hardware [int]
		CLIENT_OUTPUT_HARDWARE, //client_output_hardware [int]
		CLIENT_CHANNEL_GROUP_ID, //client_channel_group_id [int]
		CLIENT_SERVER_GROUPS, //client_servergroups [string (list of numbers, separated by a comma)]
		CLIENT_FIRSTCONNECTED, //client_created [int (UTC timestamp of first connection)]
		CLIENT_LASTCONNECTED, //client_lastconnected [int (UTC timestamp of last connection)]
		CLIENT_TOTALCONNECTIONS, //client_totalconnections [int]
		CLIENT_AWAY, //client_away [int]
		CLIENT_AWAY_MESSAGE, //client_away_message [string]
		CLIENT_AVATAR, //client_flag_avatar [string] (CRC checksum?)
		CLIENT_TALK_POWER, //client_talk_power [int]
		CLIENT_TALK_REQUEST, //client_talk_request [int]
		CLIENT_TALK_REQUEST_MSG, //client_talk_request_msg [string]
		CLIENT_IS_TALKER, //client_is_talker [int]
		CLIENT_IS_PRIORITY_SPEAKER, //client_is_priority_speaker [int]
		CLIENT_DESCRIPTION, //client_description [string]
		CLIENT_IS_CHANNEL_COMMANDER, //client_is_channel_commander [int]
		CLIENT_ICON_ID, //client_icon_id [int]
		CLIENT_COUNTRY, //client_country [string]
		CLIENT_IDLE_TIME, //client_idle_time [int] (seconds)
		CLIENT_IS_RECORDING, //client_is_recording [int]
	};

	Id_t DatabaseId = Invalid;
	string
		Uid,
		IpAddress;

	Channel::Id_t CurrentChannel = Channel::Invalid;
};
typedef shared_ptr<Client> Client_t;

class CServer : public CSingleton <CServer>
{
	friend class CSingleton <CServer>;
private: //variables
	unordered_map<Channel::Id_t, Channel_t> m_Channels;
	Channel::Id_t m_DefaultChannel = Channel::Invalid;
	mutex m_ChannelMtx;

	unordered_map<Client::Id_t, Client_t> m_Clients;
	mutex m_ClientMtx;

	atomic<bool> m_IsLoggedIn;

	unsigned int m_ServerId = 0;

	boost::lockfree::spsc_queue<
			string,
			boost::lockfree::fixed_sized<true>,
			boost::lockfree::capacity<32678>
		> m_QueriedData;
	string m_ActiveQueriedData;

private: //constructor / deconstructor
	CServer() :
		m_IsLoggedIn(false)
	{}
	~CServer() = default;

private: //functions (internal)
	void Initialize();

public: //server functions
	bool Login(string login, string pass);
	bool ChangeNickname(string nickname);
	inline bool IsLoggedIn() const
	{
		return m_IsLoggedIn;
	}

	bool SendServerMessage(string msg);

public: //data query functions
	bool QueryChannelData(Channel::Id_t cid, Channel::QueryData data, Callback_t callback);
	bool QueryClientData(Client::Id_t clid, Client::QueryData data, Callback_t callback);
	bool GetQueriedData(string &dest);
	bool GetQueriedData(int &dest);

public: //channel functions
	bool CreateChannel(string name, Channel::Types type = Channel::Types::TEMPORARY,
		int maxusers = -1, Channel::Id_t pcid = Channel::Invalid, Channel::Id_t ocid = Channel::Invalid,
		int talkpower = 0);
	bool DeleteChannel(Channel::Id_t cid);
	Channel::Id_t GetChannelIdByName(string name);
	inline bool IsValidChannel(Channel::Id_t cid)
	{
		boost::lock_guard<mutex> channel_mtx_guard(m_ChannelMtx);
		return (m_Channels.find(cid) != m_Channels.end());
	}
	bool SetChannelName(Channel::Id_t cid, string name);
	string GetChannelName(Channel::Id_t cid);
	bool SetChannelDescription(Channel::Id_t cid, string desc);
	bool SetChannelType(Channel::Id_t cid, Channel::Types type);
	Channel::Types GetChannelType(Channel::Id_t cid);
	bool SetChannelPassword(Channel::Id_t cid, string password);
	bool HasChannelPassword(Channel::Id_t cid);
	bool SetChannelRequiredTalkPower(Channel::Id_t cid, int talkpower);
	int GetChannelRequiredTalkPower(Channel::Id_t cid);
	bool SetChannelUserLimit(Channel::Id_t cid, int maxusers);
	int GetChannelUserLimit(Channel::Id_t cid);
	bool SetChannelParentId(Channel::Id_t cid, Channel::Id_t pcid);
	Channel::Id_t GetChannelParentId(Channel::Id_t cid);
	bool SetChannelOrderId(Channel::Id_t cid, Channel::Id_t ocid);
	Channel::Id_t GetChannelOrderId(Channel::Id_t cid);
	inline Channel::Id_t GetDefaultChannelId()
	{
		boost::lock_guard<mutex> channel_mtx_guard(m_ChannelMtx);
		return m_DefaultChannel;
	}

public: //client functions
	Client::Id_t GetClientIdByUid(string uid);
	Client::Id_t GetClientIdByIpAddress(string ip);
	inline bool IsValidClient(Client::Id_t clid)
	{
		boost::lock_guard<mutex> client_mtx_guard(m_ClientMtx);
		return (m_Clients.find(clid) != m_Clients.end());
	}
	string GetClientUid(Client::Id_t clid);
	Client::Id_t GetClientDatabaseId(Client::Id_t clid);
	Channel::Id_t GetClientChannelId(Client::Id_t clid);
	string GetClientIpAddress(Client::Id_t clid);

	bool KickClient(Client::Id_t clid, Client::KickTypes type, string reasonmsg);
	bool BanClient(string uid, int seconds, string reasonmsg);
	bool MoveClient(Client::Id_t clid, Channel::Id_t cid);

	bool SetClientChannelGroup(Client::Id_t clid, int groupid, Channel::Id_t cid);
	bool AddClientToServerGroup(Client::Id_t clid, int groupid);
	bool RemoveClientFromServerGroup(Client::Id_t clid, int groupid);
	bool SetClientTalkerStatus(Client::Id_t clid, bool status);
	bool SetClientDescription(Client::Id_t clid, string desc);

	bool PokeClient(Client::Id_t clid, string msg);
	bool SendClientMessage(Client::Id_t clid, string msg);

public: //network callbacks
	void OnLogin(vector<string> &res);
	void OnChannelList(vector<string> &res);
	void OnClientList(vector<string> &res);

public: //event callbacks
	void OnChannelCreated(boost::smatch &result);
	void OnChannelDeleted(boost::smatch &result);
	void OnChannelReorder(boost::smatch &result);
	void OnChannelMoved(boost::smatch &result);
	void OnChannelRenamed(boost::smatch &result);
	void OnChannelPasswordToggled(boost::smatch &result);
	void OnChannelPasswordChanged(boost::smatch &result);
	void OnChannelTypeChanged(boost::smatch &result);
	void OnChannelSetDefault(boost::smatch &result);
	void OnChannelMaxClientsChanged(boost::smatch &result);
	void OnChannelRequiredTalkPowerChanged(boost::smatch &result);

	void OnClientConnect(boost::smatch &result);
	void OnClientDisconnect(boost::smatch &result);
	void OnClientMoved(boost::smatch &result);
	void OnClientServerText(boost::smatch &result);
	void OnClientPrivateText(boost::smatch &result);
};

#endif // INC_SERVER_H


text:
#include "Avatar.h"

#include <fstream>
#include <boost/filesystem.hpp>

#include "Version.h"
#include "Util.h"

namespace ark {
Eigen::VectorXd AvatarPoseSequence::getFrame(size_t frame_id) const {
    if (preloaded) return data.col(frame_id);
    std::ifstream ifs(sequencePath, std::ios::in | std::ios::binary);
    ifs.seekg(frame_id * frameSize * sizeof(double), std::ios_base::beg);
    Eigen::VectorXd result(frameSize);
    ifs.read(reinterpret_cast<char*>(result.data()),
             frameSize * sizeof(double));
    return result;
}

AvatarPoseSequence::AvatarPoseSequence(const std::string& pose_sequence_path) {
    using namespace boost::filesystem;
    path seqPath =
        pose_sequence_path.empty()
            ? util::resolveRootPath("data/avatar-mocap/cmu-mocap.dat")
            : pose_sequence_path;
    path metaPath(std::string(seqPath.string()).append(".txt"));

    if (!exists(seqPath) || !exists(metaPath)) {
        numFrames = 0;
        return;
    }
    sequencePath = seqPath.string();

    std::ifstream metaIfs(metaPath.string());
    size_t nSubseq, frameSizeBytes, subseqStart;
    metaIfs >> nSubseq >> numFrames >> frameSizeBytes;
    std::string subseqName;
    for (int sid = 0; sid < nSubseq; ++sid) {
        metaIfs >> subseqStart >> subseqName;
        subsequences[subseqName] = subseqStart / frameSizeBytes;
    }
    metaIfs.close();

    frameSize = frameSizeBytes / sizeof(double);
}

void AvatarPoseSequence::poseAvatar(Avatar& ava, size_t frame_id) const {
    if (preloaded) {
        auto frameData = data.col(frame_id);
        ava.p.noalias() = frameData.head<3>();
        Eigen::Quaterniond q;
        for (int i = 0; i < ava.r.size(); ++i) {
            q.coeffs().noalias() = frameData.segment<4>(i * 4 + 3);
            ava.r[i].noalias() = q.toRotationMatrix();
        }
    } else {
        Eigen::VectorXd frameData = getFrame(frame_id);
        ava.p = frameData.head<3>();
        Eigen::Quaterniond q;
        for (int i = 0; i < ava.r.size(); ++i) {
            q.coeffs().noalias() = frameData.segment<4>(i * 4 + 3);
            ava.r[i].noalias() = q.toRotationMatrix();
        }
    }
}

void AvatarPoseSequence::preload() {
    data.resize(frameSize, numFrames);
    std::ifstream ifs(sequencePath, std::ios::in | std::ios::binary);
    ifs.read(reinterpret_cast<char*>(data.data()),
             numFrames * frameSize * sizeof(double));
    preloaded = true;
}
}  // namespace ark

ext: cpp

max_stars_repo_name: jyuatsfl/avatar


text:
//  Copyright (c) 2014 Anuj R. Sharma
//  Copyright (c) 2014-2015 Hartmut Kaiser
//
//  Distributed under the Boost Software License, Version 1.0. (See accompanying
//  file LICENSE_1_0.txt or copy at http:// ww.boost.org/LICENSE_1_0.txt)

#ifndef HPX_VECTOR_SEGMENTED_ITERATOR_HPP
#define HPX_VECTOR_SEGMENTED_ITERATOR_HPP

/// \file hpx/components/vector/segmented_iterator.hpp
/// \brief This file contains the implementation of iterators for hpx::vector.

 // The idea for these iterators is taken from
 // http://afstern.org/matt/segmented.pdf.

#include <hpx/include/naming.hpp>
#include <hpx/include/util.hpp>
#include <hpx/include/traits.hpp>
#include <hpx/include/serialization.hpp>

#include <hpx/components/containers/vector/partition_vector_component.hpp>

#include <cstdint>
#include <iterator>
#include <limits>
#include <type_traits>

#include <boost/integer.hpp>
#include <boost/iterator/iterator_facade.hpp>
#include <boost/iterator/iterator_adaptor.hpp>
#include <boost/iterator/filter_iterator.hpp>

namespace hpx
{
    ///////////////////////////////////////////////////////////////////////////
    template <typename T> class vector;       // forward declaration

    template <typename T> class local_vector_iterator;
    template <typename T> class const_local_vector_iterator;

    template <typename T, typename BaseIter>
    class local_raw_vector_iterator;

    template <typename T> class vector_iterator;
    template <typename T> class const_vector_iterator;

    template <typename T, typename BaseIter>
    class segment_vector_iterator;
    template <typename T, typename BaseIter>
    class const_segment_vector_iterator;

    template <typename T, typename BaseIter>
    class local_segment_vector_iterator;

    namespace server
    {
        template <typename T> class partition_vector;
    }

    ///////////////////////////////////////////////////////////////////////////
    // This class wraps plain a vector<>::iterator or vector<>::const_iterator
    template <typename T, typename BaseIter>
    class local_raw_vector_iterator
      : public boost::iterator_adaptor<
            local_raw_vector_iterator<T, BaseIter>, BaseIter
        >
    {
    private:
        typedef boost::iterator_adaptor<
                local_raw_vector_iterator<T, BaseIter>, BaseIter
            > base_type;
        typedef BaseIter base_iterator;

    public:
        typedef local_vector_iterator<T> local_iterator;
        typedef const_local_vector_iterator<T> local_const_iterator;

        local_raw_vector_iterator(base_iterator const& it,
                boost::shared_ptr<server::partition_vector<T> > const& data)
          : base_type(it), data_(data)
        {}

        local_iterator remote()
        {
            HPX_ASSERT(data_);
            std::size_t local_index = std::distance(data_->begin(), this->base());
            return local_iterator(partition_vector<T>(data_->get_gid()),
                local_index, data_);
        }
        local_const_iterator remote() const
        {
            HPX_ASSERT(data_);
            std::size_t local_index = std::distance(data_->begin(), this->base());
            return local_const_iterator(partition_vector<T>(data_->get_gid()),
                local_index, data_);
        }

    private:
        boost::shared_ptr<server::partition_vector<T> > data_;
    };

    ///////////////////////////////////////////////////////////////////////////
    namespace detail
    {
        ///////////////////////////////////////////////////////////////////////
        template <typename T>
        struct local_vector_value_proxy
        {
            local_vector_value_proxy(local_vector_iterator<T> const& it)
              : it_(it)
            {}

            operator T() const
            {
                if (!it_.get_data())
                {
                    return it_.get_partition().get_value_sync(
                        it_.get_local_index());
                }
                return *(it_.get_data()->begin() + it_.get_local_index());
            }

            template <typename T_>
            local_vector_value_proxy& operator=(T_ && value)
            {
                if (!it_.get_data())
                {
                    it_.get_partition().set_value_sync(
                        it_.get_local_index(), std::forward<T_>(value));
                }
                else
                {
                    *(it_.get_data()->begin() + it_.get_local_index()) =
                        std::forward<T_>(value);
                }
                return *this;
            }

            local_vector_iterator<T> const& it_;
        };

        template <typename T>
        struct const_local_vector_value_proxy
        {
            const_local_vector_value_proxy(
                    const_local_vector_iterator<T> const& it)
              : it_(it)
            {}

            operator T() const
            {
                if (!it_.get_data())
                {
                    return it_.get_partition().get_value_sync(
                        it_.get_local_index());
                }
                return *it_.local();
            }

            const_local_vector_iterator<T> const& it_;
        };

        ///////////////////////////////////////////////////////////////////////
        template <typename T>
        struct vector_value_proxy
        {
            vector_value_proxy(hpx::vector<T>& v, std::size_t index)
              : v_(v), index_(index)
            {}

            operator T() const
            {
                return v_.get_value_sync(index_);
            }

            template <typename T_>
            vector_value_proxy& operator=(T_ && value)
            {
                v_.set_value_sync(index_, std::forward<T_>(value));
                return *this;
            }

            vector<T>& v_;
            std::size_t index_;
        };
    }

    ///////////////////////////////////////////////////////////////////////////
    /// This class implements the local iterator functionality for the
    /// partitioned backend of a hpx::vector.
    template <typename T>
    class local_vector_iterator
      : public boost::iterator_facade<
            local_vector_iterator<T>, T, std::random_access_iterator_tag,
            detail::local_vector_value_proxy<T>
        >
    {
    private:
        typedef boost::iterator_facade<
                local_vector_iterator<T>, T, std::random_access_iterator_tag,
                detail::local_vector_value_proxy<T>
            > base_type;

    public:
        typedef std::size_t size_type;

        // constructors
        local_vector_iterator()
          : partition_(), local_index_(size_type(-1))
        {}

        local_vector_iterator(partition_vector<T> partition,
                size_type local_index,
                boost::shared_ptr<server::partition_vector<T> > const& data)
          : partition_(partition),
            local_index_(local_index),
            data_(data)
        {}

        typedef local_raw_vector_iterator<
                T, typename std::vector<T>::iterator
            > local_raw_iterator;
        typedef local_raw_vector_iterator<
                T, typename std::vector<T>::const_iterator
            > local_raw_const_iterator;

        ///////////////////////////////////////////////////////////////////////
        local_raw_iterator local()
        {
            HPX_ASSERT(data_);
            return local_raw_iterator(data_->begin() + local_index_, data_);
        }
        local_raw_const_iterator local() const
        {
            HPX_ASSERT(data_);
            return local_raw_iterator(data_->cbegin() + local_index_, data_);
        }

    private:
        friend class hpx::serialization::access;

        template <typename Archive>
        void load(Archive& ar, unsigned version)
        {
            ar & partition_ & local_index_;
            if (partition_)
                data_ = partition_.get_ptr();
        }
        template <typename Archive>
        void save(Archive& ar, unsigned version) const
        {
            ar & partition_ & local_index_;
        }

        HPX_SERIALIZATION_SPLIT_MEMBER()

    protected:
        friend class boost::iterator_core_access;

        bool equal(local_vector_iterator const& other) const
        {
            return partition_ == other.partition_ &&
                local_index_ == other.local_index_;
        }

        typename base_type::reference dereference() const
        {
            return detail::local_vector_value_proxy<T>(*this);
        }

        void increment()
        {
            ++local_index_;
        }

        void decrement()
        {
            --local_index_;
        }

        void advance(std::ptrdiff_t n)
        {
            local_index_ += n;
        }

        std::ptrdiff_t distance_to(local_vector_iterator const& other) const
        {
            HPX_ASSERT(partition_ == other.partition_);
            return other.local_index_ - local_index_;
        }

    public:
        partition_vector<T>& get_partition() { return partition_; }
        partition_vector<T> get_partition() const { return partition_; }

        size_type get_local_index() const { return local_index_; }

        boost::shared_ptr<server::partition_vector<T> >& get_data()
        {
            return data_;
        }
        boost::shared_ptr<server::partition_vector<T> > const& get_data() const
        {
            return data_;
        }

    protected:
        // refer to a partition of the vector
        partition_vector<T> partition_;

        // local position in the referenced partition
        size_type local_index_;

        // caching address of component
        boost::shared_ptr<server::partition_vector<T> > data_;
    };

    template <typename T>
    class const_local_vector_iterator
      : public boost::iterator_facade<
            const_local_vector_iterator<T>, T const,
            std::random_access_iterator_tag,
            detail::const_local_vector_value_proxy<T>
        >
    {
    private:
        typedef boost::iterator_facade<
                const_local_vector_iterator<T>, T const,
                std::random_access_iterator_tag,
                detail::const_local_vector_value_proxy<T>
            > base_type;

    public:
        typedef std::size_t size_type;

        // constructors
        const_local_vector_iterator()
          : partition_(), local_index_(size_type(-1))
        {}

        const_local_vector_iterator(partition_vector<T> partition,
                size_type local_index,
                boost::shared_ptr<server::partition_vector<T> > const& data)
          : partition_(partition),
            local_index_(local_index),
            data_(data)
        {}

        typedef local_raw_vector_iterator<
                T, typename std::vector<T>::const_iterator
            > local_raw_iterator;
        typedef local_raw_iterator local_raw_const_iterator;

        ///////////////////////////////////////////////////////////////////////
        local_raw_iterator local()
        {
            HPX_ASSERT(data_);
            return local_raw_iterator(data_->cbegin() + local_index_, data_);
        }
        local_raw_const_iterator local() const
        {
            HPX_ASSERT(data_);
            return local_raw_const_iterator(data_->cbegin() + local_index_, data_);
        }

    private:
        friend class hpx::serialization::access;

        template <typename Archive>
        void load(Archive& ar, unsigned version)
        {
            ar & partition_ & local_index_;
            if (partition_)
                data_ = partition_.get_ptr();
        }
        template <typename Archive>
        void save(Archive& ar, unsigned version) const
        {
            ar & partition_ & local_index_;
        }

        HPX_SERIALIZATION_SPLIT_MEMBER()

    protected:
        friend class boost::iterator_core_access;

        bool equal(const_local_vector_iterator const& other) const
        {
            return partition_ == other.partition_ &&
                local_index_ == other.local_index_;
        }

        typename base_type::reference dereference() const
        {
            return detail::const_local_vector_value_proxy<T>(*this);
        }

        void increment()
        {
            ++local_index_;
        }

        void decrement()
        {
            --local_index_;
        }

        void advance(std::ptrdiff_t n)
        {
            local_index_ += n;
        }

        std::ptrdiff_t distance_to(const_local_vector_iterator const& other) const
        {
            HPX_ASSERT(partition_ == other.partition_);
            return other.local_index_ - local_index_;
        }

    public:
        partition_vector<T> const& get_partition() const { return partition_; }
        size_type get_local_index() const { return local_index_; }

        boost::shared_ptr<server::partition_vector<T> >& get_data()
        {
            return data_;
        }
        boost::shared_ptr<server::partition_vector<T> > const& get_data() const
        {
            return data_;
        }

    protected:
        // refer to a partition of the vector
        partition_vector<T> partition_;

        // local position in the referenced partition
        size_type local_index_;

        // caching address of component
        boost::shared_ptr<server::partition_vector<T> > data_;
    };

    ///////////////////////////////////////////////////////////////////////////
    /// This class implement the segmented iterator for the hpx::vector
    template <typename T, typename BaseIter>
    class segment_vector_iterator
      : public boost::iterator_adaptor<
            segment_vector_iterator<T, BaseIter>, BaseIter
        >
    {
    private:
        typedef boost::iterator_adaptor<
                segment_vector_iterator<T, BaseIter>, BaseIter
            > base_type;

    public:
        segment_vector_iterator(BaseIter const& it, vector<T>* data = 0)
          : base_type(it), data_(data)
        {}

        vector<T>* get_data() { return data_; }
        vector<T> const* get_data() const { return data_; }

        bool is_at_end() const
        {
            return data_ == 0 ||
                this->base_type::base_reference() == data_->partitions_.end();
        }

    private:
        vector<T>* data_;
    };

    template <typename T, typename BaseIter>
    class const_segment_vector_iterator
      : public boost::iterator_adaptor<
            const_segment_vector_iterator<T, BaseIter>, BaseIter
        >
    {
    private:
        typedef boost::iterator_adaptor<
                const_segment_vector_iterator<T, BaseIter>, BaseIter
            > base_type;

    public:
        const_segment_vector_iterator(BaseIter const& it, vector<T> const* data = 0)
          : base_type(it), data_(data)
        {}

        vector<T> const* get_data() const { return data_; }

        bool is_at_end() const
        {
            return data_ == 0 ||
                this->base_type::base_reference() == data_->partitions_.end();
        }

    private:
        vector<T> const* data_;
    };

    ///////////////////////////////////////////////////////////////////////////
    namespace detail
    {
        template <typename BaseIterator>
        struct is_requested_locality
        {
            typedef typename std::iterator_traits<BaseIterator>::reference
                reference;

            is_requested_locality(boost::uint32_t locality_id =
                    naming::invalid_locality_id)
              : locality_id_(locality_id)
            {}

            bool operator()(reference val) const
            {
                return locality_id_ == naming::invalid_locality_id ||
                       locality_id_ == val.locality_id_;
            }

            boost::uint32_t locality_id_;
        };
    }

    /// This class implement the local segmented iterator for the hpx::vector
    template <typename T, typename BaseIter>
    class local_segment_vector_iterator
      : public boost::iterator_adaptor<
            local_segment_vector_iterator<T, BaseIter>, BaseIter,
            std::vector<T>, std::forward_iterator_tag
        >
    {
    private:
        typedef boost::iterator_adaptor<
                local_segment_vector_iterator<T, BaseIter>, BaseIter,
                std::vector<T>, std::forward_iterator_tag
            > base_type;
        typedef detail::is_requested_locality<BaseIter> predicate;

    public:
        local_segment_vector_iterator(BaseIter const& end)
          : base_type(end), predicate_(), end_(end)
        {}

        local_segment_vector_iterator(
                BaseIter const& it, BaseIter const& end,
                boost::uint32_t locality_id)
          : base_type(it), predicate_(locality_id), end_(end)
        {
            satisfy_predicate();
        }

        bool is_at_end() const
        {
            return !data_ || this->base() == end_;
        }

    private:
        friend class boost::iterator_core_access;

        typename base_type::reference dereference() const
        {
            HPX_ASSERT(!is_at_end());
            return data_->get_data();
        }

        void increment()
        {
            ++(this->base_reference());
            satisfy_predicate();
        }

        void satisfy_predicate()
        {
            while (this->base() != end_ && !predicate_(*this->base()))
                ++(this->base_reference());

            if (this->base() != end_)
                data_ = this->base()->local_data_;
            else
                data_.reset();
        }

    private:
        boost::shared_ptr<server::partition_vector<T> > data_;
        predicate predicate_;
        BaseIter end_;
    };

    ///////////////////////////////////////////////////////////////////////////
    /// This class implements the (global) iterator functionality for hpx::vector.
    template <typename T>
    class vector_iterator
      : public boost::iterator_facade<
            vector_iterator<T>, T, std::random_access_iterator_tag,
            detail::vector_value_proxy<T>
        >
    {
    private:
        typedef boost::iterator_facade<
                vector_iterator<T>, T, std::random_access_iterator_tag,
                detail::vector_value_proxy<T>
            > base_type;

    public:
        typedef std::size_t size_type;
        typedef typename vector<T>::segment_iterator segment_iterator;
        typedef typename vector<T>::local_segment_iterator local_segment_iterator;
        typedef typename vector<T>::local_iterator local_iterator;

        // constructors
        vector_iterator()
          : data_(0), global_index_(size_type(-1))
        {}

        vector_iterator(vector<T>* data, size_type global_index)
          : data_(data), global_index_(global_index)
        {}

        vector<T>* get_data() { return data_; }
        vector<T> const* get_data() const { return data_; }

        size_type get_global_index() const { return global_index_; }

    protected:
        friend class boost::iterator_core_access;

        bool equal(vector_iterator const& other) const
        {
            return data_ == other.data_ && global_index_ == other.global_index_;
        }

        typename base_type::reference dereference() const
        {
            HPX_ASSERT(data_);
            return detail::vector_value_proxy<T>(*data_, global_index_);
        }

        void increment()
        {
            HPX_ASSERT(data_);
            ++global_index_;
        }

        void decrement()
        {
            HPX_ASSERT(data_);
            --global_index_;
        }

        void advance(std::ptrdiff_t n)
        {
            HPX_ASSERT(data_);
            global_index_ += n;
        }

        std::ptrdiff_t distance_to(vector_iterator const& other) const
        {
            HPX_ASSERT(data_ && other.data_);
            HPX_ASSERT(data_ == other.data_);
            return other.global_index_ - global_index_;
        }

    protected:
        // refer to the vector
        vector<T>* data_;

        // global position in the referenced vector
        size_type global_index_;
    };

    ///////////////////////////////////////////////////////////////////////////
    template <typename T>
    class const_vector_iterator
      : public boost::iterator_facade<
            const_vector_iterator<T>, T const,
            std::random_access_iterator_tag, T const
        >
    {
    private:
        typedef boost::iterator_facade<
                const_vector_iterator<T>, T const,
                std::random_access_iterator_tag, T const
            > base_type;

    public:
        typedef std::size_t size_type;
        typedef typename vector<T>::const_segment_iterator segment_iterator;
        typedef typename vector<T>::const_local_segment_iterator local_segment_iterator;
        typedef typename vector<T>::const_local_iterator local_iterator;

        // constructors
        const_vector_iterator()
          : data_(0), global_index_(size_type(-1))
        {}

        const_vector_iterator(vector<T> const* data, size_type global_index)
          : data_(data), global_index_(global_index)
        {}

        vector<T> const* get_data() const { return data_; }
        size_type get_global_index() const { return global_index_; }

    protected:
        friend class boost::iterator_core_access;

        bool equal(const_vector_iterator const& other) const
        {
            return data_ == other.data_ && global_index_ == other.global_index_;
        }

        typename base_type::reference dereference() const
        {
            HPX_ASSERT(data_);
            return data_->get_value_sync(global_index_);
        }

        void increment()
        {
            HPX_ASSERT(data_);
            ++global_index_;
        }

        void decrement()
        {
            HPX_ASSERT(data_);
            --global_index_;
        }

        void advance(std::ptrdiff_t n)
        {
            HPX_ASSERT(data_);
            global_index_ += n;
        }

        std::ptrdiff_t distance_to(const_vector_iterator const& other) const
        {
            HPX_ASSERT(data_ && other.data_);
            HPX_ASSERT(data_ == other.data_);
            return other.global_index_ - global_index_;
        }

    protected:
        // refer to the vector
        vector<T> const* data_;

        // global position in the referenced vector
        size_type global_index_;
    };
}

///////////////////////////////////////////////////////////////////////////////
namespace hpx { namespace traits
{
    template <typename T>
    struct segmented_iterator_traits<vector_iterator<T> >
    {
        typedef std::true_type is_segmented_iterator;

        typedef vector_iterator<T> iterator;
        typedef typename iterator::segment_iterator segment_iterator;
        typedef typename iterator::local_segment_iterator local_segment_iterator;
        typedef typename iterator::local_iterator local_iterator;

        typedef typename local_iterator::local_raw_iterator local_raw_iterator;

        //  Conceptually this function is supposed to denote which segment
        //  the iterator is currently pointing to (i.e. just global iterator).
        static segment_iterator segment(iterator iter)
        {
            return iter.get_data()->get_segment_iterator(
                iter.get_global_index());
        }

        //  This function should specify which is the current segment and
        //  the exact position to which local iterator is pointing.
        static local_iterator local(iterator iter)
        {
            HPX_ASSERT(iter.get_data());    // avoid dereferencing end iterator
            return iter.get_data()->get_local_iterator(
                iter.get_global_index());
        }

        //  Build a full iterator from the segment and local iterators
        static iterator compose(segment_iterator seg_iter,
            local_iterator local_iter)
        {
            vector<T>* data = seg_iter.get_data();
            std::size_t index = local_iter.get_local_index();
            return iterator(data, data->get_global_index(seg_iter, index));
        }

        //  This function should specify the local iterator which is at the
        //  beginning of the partition.
        static local_iterator begin(segment_iterator seg_iter)
        {
            std::size_t offset = 0;
            if (seg_iter.is_at_end())
            {
                // return iterator to the end of last segment
                --seg_iter;
                offset = seg_iter.base()->size_;
            }

            return local_iterator(seg_iter.base()->partition_, offset,
                seg_iter.base()->local_data_);
        }

        //  This function should specify the local iterator which is at the
        //  end of the partition.
        static local_iterator end(segment_iterator seg_iter)
        {
            if (seg_iter.is_at_end())
                --seg_iter;     // return iterator to the end of last segment

            return local_iterator(seg_iter.base()->partition_,
                seg_iter.base()->size_, seg_iter.base()->local_data_);
        }

        //  This function should specify the local iterator which is at the
        //  beginning of the partition data.
        static local_raw_iterator begin(local_segment_iterator const& seg_iter)
        {
            return local_raw_iterator(seg_iter->begin(), seg_iter.base()->local_data_);
        }

        //  This function should specify the local iterator which is at the
        //  end of the partition data.
        static local_raw_iterator end(local_segment_iterator const& seg_iter)
        {
            return local_raw_iterator(seg_iter->end(), seg_iter.base()->local_data_);
        }

        // Extract the base id for the segment referenced by the given segment
        // iterator.
        static id_type get_id(segment_iterator const& iter)
        {
            return iter->get_id();
        }
    };

    template <typename T>
    struct segmented_iterator_traits<const_vector_iterator<T> >
    {
        typedef std::true_type is_segmented_iterator;

        typedef const_vector_iterator<T> iterator;
        typedef typename iterator::segment_iterator segment_iterator;
        typedef typename iterator::local_segment_iterator local_segment_iterator;
        typedef typename iterator::local_iterator local_iterator;

        typedef typename local_iterator::local_raw_iterator local_raw_iterator;

        //  Conceptually this function is supposed to denote which segment
        //  the iterator is currently pointing to (i.e. just global iterator).
        static segment_iterator segment(iterator iter)
        {
            return iter.get_data()->get_const_segment_iterator(
                iter.get_global_index());
        }

        //  This function should specify which is the current segment and
        //  the exact position to which local iterator is pointing.
        static local_iterator local(iterator const& iter)
        {
            HPX_ASSERT(iter.get_data());    // avoid dereferencing end iterator
            return iter.get_data()->get_const_local_iterator(
                iter.get_global_index());
        }

        //  Build a full iterator from the segment and local iterators
        static iterator compose(segment_iterator const& seg_iter,
            local_iterator const& local_iter)
        {
            vector<T> const* data = seg_iter.get_data();
            std::size_t index = local_iter.get_local_index();
            return iterator(data, data->get_global_index(seg_iter, index));
        }

        //  This function should specify the local iterator which is at the
        //  beginning of the partition.
        static local_iterator begin(segment_iterator seg_iter)
        {
            std::size_t offset = 0;
            if (seg_iter.is_at_end())
            {
                // return iterator to the end of last segment
                --seg_iter;
                offset = seg_iter.base()->size_;
            }

            return local_iterator(seg_iter.base()->partition_, offset,
                seg_iter.base()->local_data_);
        }

        //  This function should specify the local iterator which is at the
        //  end of the partition.
        static local_iterator end(segment_iterator seg_iter)
        {
            if (seg_iter.is_at_end())
                --seg_iter;     // return iterator to the end of last segment

            return local_iterator(seg_iter.base()->partition_,
                seg_iter.base()->size_, seg_iter.base()->local_data_);
        }

        //  This function should specify the local iterator which is at the
        //  beginning of the partition data.
        static local_raw_iterator begin(local_segment_iterator const& seg_iter)
        {
            return local_raw_iterator(seg_iter->cbegin(),
                seg_iter.base()->local_data_);
        }

        //  This function should specify the local iterator which is at the
        //  end of the partition data.
        static local_raw_iterator end(local_segment_iterator const& seg_iter)
        {
            return local_raw_iterator(seg_iter->cend(),
                seg_iter.base()->local_data_);
        }

        // Extract the base id for the segment referenced by the given segment
        // iterator.
        static id_type get_id(segment_iterator const& iter)
        {
            return iter->get_id();
        }
    };

    ///////////////////////////////////////////////////////////////////////////
    // Some 'remote' iterators need to be mapped before being applied to the
    // local algorithms.
    template <typename T>
    struct segmented_local_iterator_traits<local_vector_iterator<T> >
    {
        typedef std::true_type is_segmented_local_iterator;

        typedef vector_iterator<T> iterator;
        typedef local_vector_iterator<T> local_iterator;
        typedef typename local_iterator::local_raw_iterator local_raw_iterator;

        // Extract base iterator from local_iterator
        static local_raw_iterator local(local_iterator it)
        {
            return it.local();
        }

        // Construct remote local_iterator from local_raw_iterator
        static local_iterator remote(local_raw_iterator it)
        {
            return it.remote();
        }
    };

    template <typename T>
    struct segmented_local_iterator_traits<const_local_vector_iterator<T> >
    {
        typedef std::true_type is_segmented_local_iterator;

        typedef const_vector_iterator<T> iterator;
        typedef const_local_vector_iterator<T> local_iterator;
        typedef typename local_iterator::local_raw_iterator local_raw_iterator;

        // Extract base iterator from local_iterator
        static local_raw_iterator local(local_iterator it)
        {
            return it.local();
        }

        // Construct remote local_iterator from local_raw_iterator
        static local_iterator remote(local_raw_iterator it)
        {
            return it.remote();
        }
    };
}}

#endif //  SEGMENTED_ITERATOR_HPP

ext: hpp

max_stars_repo_name: akemp/hpx

max_issues_repo_name: akemp/hpx

max_forks_repo_name: akemp/hpx


text:
#ifndef PYTHONIC_INCLUDE_NUMPY_SIGN_HPP
#define PYTHONIC_INCLUDE_NUMPY_SIGN_HPP

#include "pythonic/include/utils/functor.hpp"
#include "pythonic/include/types/ndarray.hpp"
#include "pythonic/include/utils/numpy_traits.hpp"
#include <boost/simd/function/sign.hpp>

PYTHONIC_NS_BEGIN

namespace numpy
{
  namespace wrapper
  {
    template <class T>
    auto sign(T const &v) -> decltype(boost::simd::sign(v))
    {
      return boost::simd::sign(v);
    }

    template <class T>
    std::complex<T> sign(std::complex<T> const &v)
    {
      if (v.real())
        return {boost::simd::sign(v.real()), 0};
      else
        return {boost::simd::sign(v.imag()), 0};
    }
  }
#define NUMPY_NARY_FUNC_NAME sign
#define NUMPY_NARY_FUNC_SYM wrapper::sign
#include "pythonic/include/types/numpy_nary_expr.hpp"
}
PYTHONIC_NS_END

#endif

max_stars_repo_name: SylvainCorlay/pythran

max_issues_repo_name: SylvainCorlay/pythran

max_forks_repo_name: SylvainCorlay/pythran


text:
// (C) Copyright 2008 CodeRage, LLC (turkanis at coderage dot com)
// (C) Copyright 2003-2007 Jonathan Turkanis
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt.)

// See http://www.boost.org/libs/iostreams for documentation.

#ifndef BOOST_IOSTREAMS_OPERATIONS_FWD_HPP_INCLUDED
#define BOOST_IOSTREAMS_OPERATIONS_FWD_HPP_INCLUDED

#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif

#include <boost/mpl/not.hpp>
#include <boost/type_traits/is_base_and_derived.hpp>

namespace mars_boost {} namespace boost = mars_boost; namespace mars_boost { namespace iostreams {

template<typename T>
struct operations;

namespace detail {

struct custom_tag { };

template<typename T>
struct is_custom
    : mpl::not_<
          is_base_and_derived< custom_tag, operations<T> >
      >
    { };

} // End namespace detail.

template<typename T>
struct operations : detail::custom_tag { };

} } // End namespaces iostreams, boost.

#endif // #ifndef BOOST_IOSTREAMS_OPERATIONS_FWD_HPP_INCLUDED //--------------//

max_stars_repo_name: jonetomtom/mars

max_issues_repo_name: qingkouwei/mediaones

max_forks_repo_name: pengjinning/mars


text:
/*
 *            Copyright 2009-2017 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.
 *
 */

// Overload of uBLAS prod function with MKL/GSL implementations

#include <boost/filesystem.hpp>
#include <boost/format.hpp>
#include <votca/ctp/logger.h>
#include <votca/tools/constants.h>
#include <votca/tools/linalg.h>
#include <votca/xtp/gwbse.h>
#include <votca/xtp/numerical_integrations.h>
#include <votca/xtp/qmpackagefactory.h>

using boost::format;
using namespace boost::filesystem;

namespace votca {
namespace xtp {
namespace ub = boost::numeric::ublas;

// +++++++++++++++++++++++++++++ //
// GWBSE MEMBER FUNCTIONS        //
// +++++++++++++++++++++++++++++ //

void GWBSE::CleanUp() {}

void GWBSE::Initialize(Property *options) {

#if (GWBSE_DOUBLE)
  CTP_LOG(ctp::logDEBUG, *_pLog) << " Compiled with full double support"
                                 << flush;
#else
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << " Compiled with float/double mixture (standard)" << flush;
#endif

  string key = Identify();

  // getting level ranges
  _ranges = options->ifExistsReturnElseReturnDefault<string>(key + ".ranges",
                                                             "default");

  // now check validity, and get rpa, qp, and bse level ranges accordingly

  if (_ranges == "factor") {
    // get factors
    _rpamaxfactor = options->get(key + ".rpamax").as<double>();
    _qpminfactor = options->get(key + ".qpmin").as<double>();
    _qpmaxfactor = options->get(key + ".qpmax").as<double>();
    _bseminfactor = options->get(key + ".bsemin").as<double>();
    _bsemaxfactor = options->get(key + ".bsemax").as<double>();
  } else if (_ranges == "explicit") {
    // get explicit numbers
    _rpamax = options->get(key + ".rpamax").as<unsigned int>();
    _qpmin = options->get(key + ".qpmin").as<unsigned int>();
    _qpmax = options->get(key + ".qpmax").as<unsigned int>();
    _bse_vmin = options->get(key + ".bsemin").as<unsigned int>();
    _bse_cmax = options->get(key + ".bsemax").as<unsigned int>();
  } else if (_ranges == "" || _ranges == "default") {
    _ranges = "default";
  } else if (_ranges == "full") {
    _ranges = "full";
  } else {
    cerr << "\nSpecified range option " << _ranges << " invalid. ";
    throw std::runtime_error(
        "\nValid options are: default,factor,explicit,full");
  }

  _ignore_corelevels = options->ifExistsReturnElseReturnDefault<bool>(
      key + ".ignore_corelevels", false);

  _bse_nmax =
      options->ifExistsReturnElseReturnDefault<int>(key + ".exctotal", 25);
  _bse_nprint =
      options->ifExistsReturnElseReturnDefault<int>(key + ".print", 25);
  _fragA = options->ifExistsReturnElseReturnDefault<int>(key + ".fragment", -1);

  string BSEtype =
      options->ifExistsReturnElseReturnDefault<string>(key + ".BSEtype", "TDA");

  if (BSEtype == "full") {
    _do_full_BSE = true;
    CTP_LOG(ctp::logDEBUG, *_pLog) << " BSE type: full" << flush;
  } else {
    _do_full_BSE = false;
    CTP_LOG(ctp::logDEBUG, *_pLog) << " BSE type: TDA" << flush;
  }

  _openmp_threads =
      options->ifExistsReturnElseReturnDefault<int>(key + ".openmp", 0);

  if (options->exists(key + ".vxc")) {
    _doVxc =
        options->ifExistsReturnElseThrowRuntimeError<bool>(key + ".vxc.dovxc");
    if (_doVxc) {
      _functional = options->ifExistsReturnElseThrowRuntimeError<string>(
          key + ".vxc.functional");
      _grid = options->ifExistsReturnElseReturnDefault<string>(
          key + ".vxc.grid", "medium");
    }
  }

  _gwbasis_name =
      options->ifExistsReturnElseThrowRuntimeError<string>(key + ".gwbasis");
  _dftbasis_name =
      options->ifExistsReturnElseThrowRuntimeError<string>(key + ".dftbasis");

  _shift = options->ifExistsReturnElseThrowRuntimeError<double>(key + ".shift");
  _g_sc_limit = options->ifExistsReturnElseReturnDefault<double>(
      key + ".g_sc_limit",
      0.00001);  // convergence criteria for qp iteration [Hartree]]
  _g_sc_max_iterations = options->ifExistsReturnElseReturnDefault<int>(
      key + ".g_sc_max_iterations",
      40);  // convergence criteria for qp iteration [Hartree]]

  _gw_sc_max_iterations = options->ifExistsReturnElseReturnDefault<int>(
      key + ".gw_sc_max_iterations",
      20);  // convergence criteria for qp iteration [Hartree]]

  _gw_sc_limit = options->ifExistsReturnElseReturnDefault<double>(
      key + ".gw_sc_limit", 0.00001);  // convergence criteria for shift it
  _iterate_gw = false;
  string _shift_type =
      options->ifExistsReturnElseThrowRuntimeError<string>(key + ".shift_type");
  if (_shift_type != "fixed") {
    _iterate_gw = true;
  }
  CTP_LOG(ctp::logDEBUG, *_pLog) << " Shift: " << _shift_type << flush;
  CTP_LOG(ctp::logDEBUG, *_pLog) << " g_sc_limit [Hartree]: " << _g_sc_limit
                                 << flush;
  if (_iterate_gw) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " gw_sc_limit [Hartree]: " << _gw_sc_limit
                                   << flush;
    // if(_g_sc_max_iterations<20){
    //_g_sc_max_iterations=20;
    // CTP_LOG(ctp::logDEBUG, *_pLog) << " Setting G iterations to 20, as it
    // speeds up the GW iterations" << flush;
    // }
  }
  _min_print_weight = options->ifExistsReturnElseReturnDefault<double>(
      key + ".bse_print_weight",
      0.2);  // print exciton WF composition weight larger that thin minimum

  // setting some defaults
  _do_qp_diag = false;
  _do_bse_singlets = false;
  _do_bse_triplets = false;
  // possible tasks
  // diagQP, singlets, triplets, all, ibse
  string _tasks_string =
      options->ifExistsReturnElseThrowRuntimeError<string>(key + ".tasks");
  if (_tasks_string.find("all") != std::string::npos) {
    _do_qp_diag = true;
    _do_bse_singlets = true;
    _do_bse_triplets = true;
  }
  if (_tasks_string.find("qpdiag") != std::string::npos) _do_qp_diag = true;
  if (_tasks_string.find("singlets") != std::string::npos)
    _do_bse_singlets = true;
  if (_tasks_string.find("triplets") != std::string::npos)
    _do_bse_triplets = true;
  _store_eh_interaction = false;
  _do_bse_diag = true;
  // special construction for ibse mode
  if (_tasks_string.find("igwbse") != std::string::npos) {
    _do_qp_diag = false;   // no qp diagonalization
    _do_bse_diag = false;  // no diagonalization of BSE Hamiltonian
    _store_eh_interaction = true;
  }

  // possible storage
  // qpPert, qpdiag_energies, qp_diag_coefficients, bse_singlet_energies,
  // bse_triplet_energies, bse_singlet_coefficients, bse_triplet_coefficients
  _store_qp_pert = true;

  _store_qp_diag = false;
  _store_bse_triplets = false;
  _store_bse_singlets = false;
  string _store_string =
      options->ifExistsReturnElseThrowRuntimeError<string>(key + ".store");
  if ((_store_string.find("all") != std::string::npos) ||
      (_store_string.find("") != std::string::npos)) {
    // store according to tasks choice
    if (_do_qp_diag) _store_qp_diag = true;
    if (_do_bse_singlets && _do_bse_diag) _store_bse_singlets = true;
    if (_do_bse_triplets && _do_bse_diag) _store_bse_triplets = true;
  }
  if (_store_string.find("qpdiag") != std::string::npos) _store_qp_diag = true;
  if (_store_string.find("singlets") != std::string::npos)
    _store_bse_singlets = true;
  if (_store_string.find("triplets") != std::string::npos)
    _store_bse_triplets = true;
  if (_store_string.find("ehint") != std::string::npos)
    _store_eh_interaction = true;

  CTP_LOG(ctp::logDEBUG, *_pLog) << " Tasks: " << flush;
  if (_do_qp_diag) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " qpdiag " << flush;
  }
  if (_do_bse_singlets) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " singlets " << flush;
  }
  if (_do_bse_triplets) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " triplets " << flush;
  }
  CTP_LOG(ctp::logDEBUG, *_pLog) << " Store: " << flush;
  if (_store_qp_diag) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " qpdiag " << flush;
  }
  if (_store_bse_singlets) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " singlets " << flush;
  }
  if (_store_bse_triplets) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " triplets " << flush;
  }
  if (_store_eh_interaction) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << " ehint " << flush;
  }

  return;
}

void GWBSE::addoutput(Property *_summary) {

  const double hrt2ev = tools::conv::hrt2ev;
  Property *_gwbse_summary = &_summary->add("GWBSE", "");
  _gwbse_summary->setAttribute("units", "eV");
  _gwbse_summary->setAttribute("DeltaHLGap",
                               (format("%1$+1.6f ") % (_shift * hrt2ev)).str());

  _gwbse_summary->setAttribute(
      "DFTEnergy", (format("%1$+1.6f ") % _orbitals->getQMEnergy()).str());
  int printlimit = _bse_nprint;  // I use this to determine how much is printed,
                                 // I do not want another option to pipe through

  Property *_dft_summary = &_gwbse_summary->add("dft", "");
  _dft_summary->setAttribute("HOMO", _homo);
  _dft_summary->setAttribute("LUMO", _homo + 1);
  int begin = _homo - printlimit;
  int end = _homo + printlimit + 1;
  if (begin < 0) {
    begin = 0;
    end = 2 * _homo + 1;
  }
  for (int state = begin; state < end; state++) {

    Property *_level_summary = &_dft_summary->add("level", "");
    _level_summary->setAttribute("number", state);
    _level_summary->add("dft_energy",
                        (format("%1$+1.6f ") %
                         ((_orbitals->MOEnergies())(_qpmin + state) * hrt2ev))
                            .str());
    _level_summary->add(
        "gw_energy",
        (format("%1$+1.6f ") % (_qp_energies(_qpmin + state) * hrt2ev)).str());

    if (_do_qp_diag) {
      // cout << "_do_qp_diag" <<_do_qp_diag<<endl;
      _level_summary->add(
          "qp_energy",
          (format("%1$+1.6f ") % (_qp_diag_energies(_qpmin + state) * hrt2ev))
              .str());
    }
  }

  if (_do_bse_singlets) {
    Property *_singlet_summary = &_gwbse_summary->add("singlets", "");
    for (int state = 0; state < printlimit; ++state) {
      Property *_level_summary = &_singlet_summary->add("level", "");
      _level_summary->setAttribute("number", state + 1);
      _level_summary->add("omega", (format("%1$+1.6f ") %
                                    (_bse_singlet_energies(state) * hrt2ev))
                                       .str());
      if (_orbitals->hasTransitionDipoles()) {

        const tools::vec &dipoles = (_orbitals->TransitionDipoles())[state];
        double f = 2 * dipoles * dipoles * _bse_singlet_energies(state) / 3.0;

        _level_summary->add("f", (format("%1$+1.6f ") % f).str());
        Property *_dipol_summary = &_level_summary->add(
            "Trdipole", (format("%1$+1.4f %2$+1.4f %3$+1.4f") % dipoles.getX() %
                         dipoles.getY() % dipoles.getZ())
                            .str());
        _dipol_summary->setAttribute("unit", "e*bohr");
        _dipol_summary->setAttribute("gauge", "length");
      }
    }
  }
  if (_do_bse_triplets) {
    Property *_triplet_summary = &_gwbse_summary->add("triplets", "");
    for (int state = 0; state < printlimit; ++state) {

      Property *_level_summary = &_triplet_summary->add("level", "");
      _level_summary->setAttribute("number", state + 1);
      _level_summary->add("omega", (format("%1$+1.6f ") %
                                    (_bse_triplet_energies(state) * hrt2ev))
                                       .str());
    }
  }
  return;
}

/*
 *    Many-body Green's fuctions theory implementation
 *
 *  data required from orbitals file
 *  - atomic coordinates
 *  - DFT molecular orbitals (energies and coeffcients)
 *  - DFT exchange-correlation potential matrix in atomic orbitals
 *  - number of electrons, number of levels
 *

 */

bool GWBSE::Evaluate() {

// set the parallelization
#ifdef _OPENMP
  if (_openmp_threads > 0) {
    omp_set_num_threads(_openmp_threads);
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " Using "
                                   << omp_get_max_threads() << " threads"
                                   << flush;
  }
#endif
  /* check which QC program was used for the DFT run
   * -> implicit info about MO coefficient storage order
   */
  string _dft_package = _orbitals->getQMpackage();
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " DFT data was created by " << _dft_package
                                 << flush;

  std::vector<ctp::QMAtom *> _atoms;
  for (const auto &atom : _orbitals->QMAtoms()) {
    if (!atom->from_environment) {
      _atoms.push_back(atom);
    }
  }
  // load DFT basis set (element-wise information) from xml file
  BasisSet dftbs;

  if (_dftbasis_name != _orbitals->getDFTbasis()) {
    throw std::runtime_error(
        "Name of the Basisset from .orb file: " + _orbitals->getDFTbasis() +
        " and from GWBSE optionfile " + _dftbasis_name + " do not agree.");
  }

  dftbs.LoadBasisSet(_dftbasis_name);
  _orbitals->setDFTbasis(_dftbasis_name);
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " Loaded DFT Basis Set "
                                 << _dftbasis_name << flush;

  // fill DFT AO basis by going through all atoms

  _dftbasis.AOBasisFill(&dftbs, _atoms, _fragA);
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " Filled DFT Basis of size "
                                 << _dftbasis.AOBasisSize() << flush;
  if (_dftbasis._AOBasisFragB > 0) {
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " FragmentA size "
                                   << _dftbasis._AOBasisFragA << flush;
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " FragmentB size "
                                   << _dftbasis._AOBasisFragB << flush;
  }

  if (_do_full_BSE)
    _orbitals->setBSEtype("full");
  else {
    _orbitals->setBSEtype("TDA");
  }
  /* Preparation of calculation parameters:
   *  - number of electrons -> index of HOMO
   *  - number of levels
   *  - highest level considered in RPA
   *  - lowest and highest level considered in GWA
   *  - lowest and highest level considered in BSE
   *  - number of excitations calculates in BSE
   */

  // convert _rpamax if needed
  _homo = _orbitals->getNumberOfElectrons() - 1;  // indexed from 0

  unsigned int _ignored_corelevels = 0;
  if (_ignore_corelevels) {
    std::string _ecpsave = _orbitals->getECP();
    _orbitals->setECP("ecp");
    int _valence_levels =
        _orbitals->FragmentNuclearCharges(_atoms.size())(0) / 2;
    _orbitals->setECP(_ecpsave);
    _ignored_corelevels = _orbitals->getNumberOfElectrons() - _valence_levels;
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " Can ignore "
                                   << _ignored_corelevels << " core levels "
                                   << flush;
  }

  _rpamin = 0;  // lowest index occ min(gwa%mmin, screening%nsum_low) ! always 1
  if (_ranges == "default" || _ranges == "full") {
    _rpamax = _orbitals->getNumberOfLevels() - 1;  // total number of levels
  } else if (_ranges == "factor") {
    _rpamax = _rpamaxfactor * _orbitals->getNumberOfLevels() -
              1;  // total number of levels
  }
  if (_rpamax > _orbitals->getNumberOfLevels() - 1) {
    _rpamax = _orbitals->getNumberOfLevels() - 1;
  }
  // convert _qpmin and _qpmax if needed
  if (_ranges == "default") {
    _qpmin = 0;              // indexed from 0
    _qpmax = 2 * _homo + 1;  // indexed from 0
  } else if (_ranges == "factor") {
    if (_orbitals->getNumberOfElectrons() -
            int(_qpminfactor * _orbitals->getNumberOfElectrons()) - 1 <
        0) {
      _qpmin = 0;
    } else {
      _qpmin = _orbitals->getNumberOfElectrons() -
               int(_qpminfactor * _orbitals->getNumberOfElectrons()) - 1;
    }
    _qpmax = _orbitals->getNumberOfElectrons() +
             int(_qpmaxfactor * _orbitals->getNumberOfElectrons()) - 1;
  } else if (_ranges == "explicit") {
    _qpmin -= 1;
    _qpmax -= 1;
  } else if (_ranges == "full") {
    _qpmin = 0;
    _qpmax = _orbitals->getNumberOfLevels() - 1;
  }
  if (_qpmax > unsigned(_orbitals->getNumberOfLevels() - 1)) {
    _qpmax = _orbitals->getNumberOfLevels() - 1;
  }

  // autoignore core levels in QP
  if (_ignore_corelevels && (_qpmin < _ignored_corelevels)) {
    _qpmin = _ignored_corelevels;
  }

  // set BSE band range indices
  // anything else would be stupid!
  _bse_vmax = _homo;
  _bse_cmin = _homo + 1;

  if (_ranges == "default") {
    _bse_vmin = 0;              // indexed from 0
    _bse_cmax = 2 * _homo + 1;  // indexed from 0
  } else if (_ranges == "factor") {
    _bse_vmin = _orbitals->getNumberOfElectrons() -
                int(_bseminfactor * _orbitals->getNumberOfElectrons()) - 1;
    if (_orbitals->getNumberOfElectrons() -
            int(_bseminfactor * _orbitals->getNumberOfElectrons()) - 1 <
        0) {
      _bse_vmin = 0;
    }
    _bse_cmax = _orbitals->getNumberOfElectrons() +
                int(_bsemaxfactor * _orbitals->getNumberOfElectrons()) - 1;
  } else if (_ranges == "explicit") {
    _bse_vmin -= 1;
    _bse_cmax -= 1;
  } else if (_ranges == "full") {
    _bse_vmin = 0;
    _bse_cmax = _orbitals->getNumberOfLevels() - 1;
  }
  if (_bse_cmax > unsigned(_orbitals->getNumberOfLevels() - 1)) {
    _bse_cmax = _orbitals->getNumberOfLevels() - 1;
  }

  // autoignore core levels in BSE
  if (_ignore_corelevels && (_bse_vmin < _ignored_corelevels)) {
    _bse_vmin = _ignored_corelevels;
  }

  _bse_vtotal = _bse_vmax - _bse_vmin + 1;
  _bse_ctotal = _bse_cmax - _bse_cmin + 1;
  _bse_size = _bse_vtotal * _bse_ctotal;

  // indexing info BSE vector index to occupied/virtual orbital
  for (unsigned _v = 0; _v < _bse_vtotal; _v++) {
    for (unsigned _c = 0; _c < _bse_ctotal; _c++) {
      _index2v.push_back(_bse_vmin + _v);
      _index2c.push_back(_bse_cmin + _c);
    }
  }

  // some QP - BSE consistency checks are required
  if (_bse_vmin < _qpmin) _qpmin = _bse_vmin;
  if (_bse_cmax > _qpmax) _qpmax = _bse_cmax;

  _qptotal = _qpmax - _qpmin + 1;
  if (_bse_nmax > int(_bse_size) || _bse_nmax < 0) _bse_nmax = int(_bse_size);
  if (_bse_nprint > _bse_nmax) _bse_nprint = _bse_nmax;

  // store information in _orbitals for later use
  _orbitals->setRPAindices(_rpamin, _rpamax);
  _orbitals->setGWAindices(_qpmin, _qpmax);
  _orbitals->setBSEindices(_bse_vmin, _bse_vmax, _bse_cmin, _bse_cmax,
                           _bse_nmax);

  // information for hybrid DFT

  _ScaHFX = -1;

  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " Set RPA level range ["
                                 << _rpamin + 1 << ":" << _rpamax + 1 << "]"
                                 << flush;
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " Set QP  level range ["
                                 << _qpmin + 1 << ":" << _qpmax + 1 << "]"
                                 << flush;
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp() << " Set BSE level range occ[" << _bse_vmin + 1 << ":"
      << _bse_vmax + 1 << "]  virt[" << _bse_cmin + 1 << ":" << _bse_cmax + 1
      << "]" << flush;
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp() << " BSE Hamiltonian has size " << _bse_size << "x"
      << _bse_size << flush;

  // process the DFT data
  // a) form the expectation value of the XC functional in MOs

  _ScaHFX = _orbitals->getScaHFX();

  ub::matrix<double> _vxc_ao;
  if (_orbitals->hasAOVxc()) {
    if (_doVxc) {
      CTP_LOG(ctp::logDEBUG, *_pLog)
          << ctp::TimeStamp()
          << " There is already a Vxc matrix loaded from DFT, did you maybe "
             "run a DFT code with outputVxc?\n I will take the external "
             "implementation"
          << flush;
      _doVxc = false;
    }
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                   << " Loaded external Vxc matrix" << flush;
    _vxc_ao = _orbitals->AOVxc();
  } else if (_doVxc) {

    NumericalIntegration _numint;
    _numint.setXCfunctional(_functional);
    double ScaHFX_temp = _numint.getExactExchange(_functional);
    if (ScaHFX_temp != _ScaHFX) {
      throw std::runtime_error(
          (boost::format("GWBSE exact exchange a=%s differs from qmpackage "
                         "exact exchange a=%s, probably your functionals are "
                         "inconsistent") %
           ScaHFX_temp % _ScaHFX)
              .str());
    }
    _numint.GridSetup(_grid, &dftbs, _atoms, &_dftbasis);
    CTP_LOG(ctp::logDEBUG, *_pLog)
        << ctp::TimeStamp()
        << " Setup grid for integration with gridsize: " << _grid << " with "
        << _numint.getGridSize() << " points, divided into "
        << _numint.getBoxesSize() << " boxes" << flush;

    CTP_LOG(ctp::logDEBUG, *_pLog)
        << ctp::TimeStamp() << " Converted DFT orbital coefficient order from "
        << _dft_package << " to XTP" << flush;
    CTP_LOG(ctp::logDEBUG, *_pLog)
        << ctp::TimeStamp() << " Integrating Vxc in VOTCA with functional "
        << _functional << flush;
    ub::matrix<double> DMAT = _orbitals->DensityMatrixGroundState();

    _vxc_ao = _numint.IntegrateVXC(DMAT);
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                   << " Calculated Vxc in VOTCA" << flush;

  } else {
    throw std::runtime_error(
        "So your DFT data contains no Vxc, if you want to proceed use the "
        "dovxc option.");
  }

  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " Set hybrid exchange factor: " << _ScaHFX
                                 << flush;

  // now get expectation values but only for those in _qpmin:_qpmax range
  ub::matrix<double> _mos =
      ub::project(_dft_orbitals, ub::range(_qpmin, _qpmax + 1),
                  ub::range(0, _dftbasis.AOBasisSize()));

  ub::matrix<double> _temp = ub::prod(_vxc_ao, ub::trans(_mos));
  _vxc = ub::prod(_mos, _temp);
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp()
      << " Calculated exchange-correlation expectation values " << flush;
  _vxc_ao.resize(0, 0);

  /// ------- actual calculation begins here -------

  // load auxiliary GW basis set (element-wise information) from xml file
  BasisSet gwbs;
  gwbs.LoadBasisSet(_gwbasis_name);
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " Loaded GW Basis Set "
                                 << _gwbasis_name << flush;

  // fill auxiliary GW AO basis by going through all atoms
  AOBasis gwbasis;
  gwbasis.AOBasisFill(&gwbs, _atoms);
  _orbitals->setGWbasis(_gwbasis_name);
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " Filled GW Basis of size "
                                 << gwbasis.AOBasisSize() << flush;

  /*
   * for the representation of 2-point functions with the help of the
   * auxiliary GW basis, its AO overlap matrix is required.
   * cf. M. Rohlfing, PhD thesis, ch. 3
   */
  AOOverlap _gwoverlap;
  // Fill overlap
  _gwoverlap.Fill(gwbasis);

  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " Filled GW Overlap matrix of dimension: "
                                 << _gwoverlap.Matrix().size1() << flush;

  /*
   *  for the calculation of Coulomb and exchange term in the self
   *  energy and electron-hole interaction, the Coulomb interaction
   *  is represented using the auxiliary GW basis set.
   *  Here, we need to prepare the Coulomb matrix expressed in
   *  the AOs of the GW basis
   */

  // get Coulomb matrix as AOCoulomb
  AOCoulomb _gwcoulomb;

  // Fill Coulomb matrix
  _gwcoulomb.Fill(gwbasis);
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " Filled GW Coulomb matrix of dimension: "
                                 << _gwcoulomb.Matrix().size1() << flush;

  // PPM is symmetric, so we need to get the sqrt of the Coulomb matrix

  ub::matrix<double> _gwoverlap_cholesky = _gwoverlap.Matrix();
  linalg_cholesky_decompose(_gwoverlap_cholesky);

// remove L^T from Cholesky
#pragma omp parallel for
  for (unsigned i = 0; i < _gwoverlap_cholesky.size1(); i++) {
    for (unsigned j = i + 1; j < _gwoverlap_cholesky.size1(); j++) {
      _gwoverlap_cholesky(i, j) = 0.0;
    }
  }

  ub::matrix<double> _gwoverlap_cholesky_inverse;  // will also be needed in PPM
                                                   // itself
  int removed =
      linalg_invert_svd(_gwoverlap_cholesky, _gwoverlap_cholesky_inverse, 1e7);
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp() << " Removed " << removed
      << " functions from gwoverlap to avoid near linear dependencies" << flush;

  int removed_functions = _gwcoulomb.Symmetrize(_gwoverlap_cholesky);
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp() << " Prepared GW Coulomb matrix for symmetric PPM"
      << flush;
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp() << " Removed " << removed_functions
      << " functions from gwcoulomb to avoid near linear dependencies" << flush;
  /* calculate 3-center integrals,  convoluted with DFT eigenvectors
   *
   *  M_mn(beta) = \int{ \psi^DFT_m(r) \phi^GW_beta(r) \psi^DFT_n d3r  }
   *             = \sum_{alpha,gamma} { c_m,alpha c_n,gamma \int
   * {\phi^DFT_alpha(r) \phi^GW_beta(r) \phi^DFT_gamma(r) d3r}  }
   *
   *  cf. M. Rohlfing, PhD thesis, ch. 3.2
   *
   */

  // --- prepare a vector (gwdacay) of matrices (orbitals, orbitals) as
  // container => M_mn
  // prepare 3-center integral object

  TCMatrix _Mmn;
  _Mmn.Initialize(gwbasis.AOBasisSize(), _rpamin, _qpmax, _rpamin, _rpamax);
  _Mmn.Fill(gwbasis, _dftbasis, _dft_orbitals);
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp()
      << " Calculated Mmn_beta (3-center-repulsion x orbitals)  " << flush;

  // make _Mmn symmetric
  _Mmn.Symmetrize(_gwcoulomb.Matrix());
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp() << " Symmetrize Mmn_beta for self-energy  " << flush;

  // for use in RPA, make a copy of _Mmn with dimensions
  // (1:HOMO)(gwabasissize,LUMO:nmax)
  TCMatrix _Mmn_RPA;
  _Mmn_RPA.Initialize(gwbasis.AOBasisSize(), _rpamin, _homo, _homo + 1,
                      _rpamax);
  RPA_prepare_threecenters(_Mmn_RPA, _Mmn);
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " Prepared Mmn_beta for RPA  " << flush;

  // fix the frequencies for PPM
  _screening_freq = ub::zero_matrix<double>(2, 2);  // two frequencies
  // first one
  _screening_freq(0, 0) = 0.0;  // real part
  _screening_freq(0, 1) = 0.0;  // imaginary part
  // second one
  _screening_freq(1, 0) = 0.0;  // real part
  _screening_freq(1, 1) = 0.5;  // imaginary part  //hartree

  // one entry to epsilon for each frequency
  _epsilon.resize(_screening_freq.size1());

  /* for automatic iteration of _shift, we need to
   * - make a copy of _Mmn
   * - calculate eps
   * - construct ppm
   * - threecenters for sigma
   * - sigma_x
   * - sigma_c
   * - test for convergence
   *
   */

  // initialize _qp_energies;
  // shift unoccupied levels by the shift
  _qp_energies = ub::zero_vector<double>(_orbitals->getNumberOfLevels());
  for (size_t i = 0; i < _qp_energies.size(); ++i) {
    _qp_energies(i) = _orbitals->MOEnergies()(i);
    if (i > _homo) {
      _qp_energies(i) += _shift;
    }
  }

  _sigma_c.resize(_qptotal);
  _sigma_x.resize(_qptotal);

  TCMatrix _Mmn_backup;
  if (_iterate_gw) {

    // make copy of _Mmn, memory++

    _Mmn_backup.Initialize(gwbasis.AOBasisSize(), _rpamin, _qpmax, _rpamin,
                           _rpamax);
    int _mnsize = _Mmn_backup.get_mtot();
    for (int _i = 0; _i < _mnsize; _i++) {
      _Mmn_backup[_i] = _Mmn[_i];
    }
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                   << " Made backup of _Mmn  " << flush;
  } else {
    _gw_sc_max_iterations = 1;
  }

  const ub::vector<double> &_dft_energies = _orbitals->MOEnergies();
  for (unsigned gw_iteration = 0; gw_iteration < _gw_sc_max_iterations;
       ++gw_iteration) {

    ub::vector<double> _qp_old_rpa = _qp_energies;
    if (_iterate_gw) {
      CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp() << " GW Iteraton "
                                     << gw_iteration + 1 << " of "
                                     << _gw_sc_max_iterations << flush;
    }

    // for symmetric PPM, we can initialize _epsilon with the overlap matrix!
    for (unsigned _i_freq = 0; _i_freq < _screening_freq.size1(); _i_freq++) {
      _epsilon[_i_freq] = _gwoverlap.Matrix();
    }

    // determine epsilon from RPA
    RPA_calculate_epsilon(_Mmn_RPA);
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                   << " Calculated epsilon via RPA  " << flush;

    // construct PPM parameters
    PPM_construct_parameters(_gwoverlap_cholesky_inverse);
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                   << " Constructed PPM parameters  " << flush;

    // prepare threecenters for Sigma
    sigma_prepare_threecenters(_Mmn);
    CTP_LOG(ctp::logDEBUG, *_pLog)
        << ctp::TimeStamp() << " Prepared threecenters for sigma  " << flush;

    sigma_diag(_Mmn);
    CTP_LOG(ctp::logDEBUG, *_pLog)
        << ctp::TimeStamp() << " Calculated diagonal part of Sigma  " << flush;
    // iterative refinement of qp energies

    double _DFTgap = _dft_energies(_homo + 1) - _dft_energies(_homo);
    double _QPgap = _qp_energies(_homo + 1) - _qp_energies(_homo);
    _shift = _QPgap - _DFTgap;

    // qp energies outside the update range are simply shifted.
    for (unsigned i = _qpmax + 1; i < _dft_energies.size(); ++i) {
      _qp_energies(i) = _dft_energies(i) + _shift;
    }

    if (_iterate_gw) {
      bool _gw_converged = true;
      ub::vector<double> diff = _qp_old_rpa - _qp_energies;
      unsigned int _l_not_converged = 0;
      double E_max = 0;
      for (unsigned l = 0; l < diff.size(); l++) {
        if (std::abs(diff(l)) > std::abs(E_max)) {
          _l_not_converged = l;
          E_max = diff(l);
        }
        if (std::abs(diff(l)) > _gw_sc_limit) {
          _gw_converged = false;
        }
      }
      double alpha = 0.0;
      _qp_energies = alpha * _qp_old_rpa + (1 - alpha) * _qp_energies;
      if (tools::globals::verbose) {
        CTP_LOG(ctp::logDEBUG, *_pLog)
            << ctp::TimeStamp() << " GW_Iteration: " << gw_iteration + 1
            << " shift=" << _shift << " E_diff max=" << E_max
            << " StateNo:" << _l_not_converged << flush;
      }

      if (_gw_converged) {
        CTP_LOG(ctp::logDEBUG, *_pLog)
            << ctp::TimeStamp() << " Converged after " << gw_iteration + 1
            << " GW iterations" << flush;
        break;
      } else if (gw_iteration == _gw_sc_max_iterations - 1) {
        // continue regardless for now, but drop WARNING
        CTP_LOG(ctp::logDEBUG, *_pLog)
            << ctp::TimeStamp() << " WARNING! GWA spectrum not converged after "
            << _gw_sc_max_iterations << " iterations." << flush;
        CTP_LOG(ctp::logDEBUG, *_pLog)
            << ctp::TimeStamp() << "          GWA level " << _l_not_converged
            << " energy changed by " << diff(_l_not_converged) << flush;
        CTP_LOG(ctp::logDEBUG, *_pLog)
            << ctp::TimeStamp()
            << "          Run continues. Inspect results carefully!" << flush;
        break;
      }

      int _mnsize = _Mmn_backup.get_mtot();
      for (int _i = 0; _i < _mnsize; _i++) {
        _Mmn[_i] = _Mmn_backup[_i];
      }
    }
  }

  sigma_offdiag(_Mmn);
  CTP_LOG(ctp::logDEBUG, *_pLog)
      << ctp::TimeStamp() << " Calculated offdiagonal part of Sigma  " << flush;
  _gwoverlap.Matrix().resize(0, 0);
  _gwoverlap_cholesky_inverse.resize(0, 0);
  _Mmn_RPA.Cleanup();
  if (_iterate_gw) {
    _Mmn_backup.Cleanup();
    CTP_LOG(ctp::logDEBUG, *_pLog)
        << ctp::TimeStamp() << " Cleaned up Overlap, MmnRPA and Mmn_backup "
        << flush;
  } else {
    CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                   << " Cleaned up Overlap and MmnRPA" << flush;
  }
  // free no longer required three-center matrices in _Mmn
  // max required is _bse_cmax (could be smaller than _qpmax)
  _Mmn.Prune(gwbasis.AOBasisSize(), _bse_vmin, _bse_cmax);

  // Output of quasiparticle energies after all is done:

  CTP_LOG(ctp::logINFO, *_pLog)
      << (format(
              "  ====== Perturbative quasiparticle energies (Hartree) ====== "))
             .str()
      << flush;
  CTP_LOG(ctp::logINFO, *_pLog)
      << (format("   DeltaHLGap = %1$+1.6f Hartree") % _shift).str() << flush;

  for (unsigned _i = 0; _i < _qptotal; _i++) {
    if ((_i + _qpmin) == _homo) {
      CTP_LOG(ctp::logINFO, *_pLog)
          << (format("  HOMO  = %1$4d DFT = %2$+1.4f VXC = %3$+1.4f S-X = "
                     "%4$+1.4f S-C = %5$+1.4f GWA = %6$+1.4f") %
              (_i + _qpmin + 1) % _dft_energies(_i + _qpmin) % _vxc(_i, _i) %
              _sigma_x(_i, _i) % _sigma_c(_i, _i) % _qp_energies(_i + _qpmin))
                 .str()
          << flush;
    } else if ((_i + _qpmin) == _homo + 1) {
      CTP_LOG(ctp::logINFO, *_pLog)
          << (format("  LUMO  = %1$4d DFT = %2$+1.4f VXC = %3$+1.4f S-X = "
                     "%4$+1.4f S-C = %5$+1.4f GWA = %6$+1.4f") %
              (_i + _qpmin + 1) % _dft_energies(_i + _qpmin) % _vxc(_i, _i) %
              _sigma_x(_i, _i) % _sigma_c(_i, _i) % _qp_energies(_i + _qpmin))
                 .str()
          << flush;

    } else {
      CTP_LOG(ctp::logINFO, *_pLog)
          << (format("  Level = %1$4d DFT = %2$+1.4f VXC = %3$+1.4f S-X = "
                     "%4$+1.4f S-C = %5$+1.4f GWA = %6$+1.4f") %
              (_i + _qpmin + 1) % _dft_energies(_i + _qpmin) % _vxc(_i, _i) %
              _sigma_x(_i, _i) % _sigma_c(_i, _i) % _qp_energies(_i + _qpmin))
                 .str()
          << flush;
    }
  }

  // store perturbative QP energy data in orbitals object (DFT, S_x,S_c, V_xc,
  // E_qp)
  if (_store_qp_pert) {
    ub::matrix<double> &_qp_energies_store = _orbitals->QPpertEnergies();
    _qp_energies_store.resize(_qptotal, 5);
    for (unsigned _i = 0; _i < _qptotal; _i++) {
      _qp_energies_store(_i, 0) = _dft_energies(_i + _qpmin);
      _qp_energies_store(_i, 1) = _sigma_x(_i, _i);
      _qp_energies_store(_i, 2) = _sigma_c(_i, _i);
      _qp_energies_store(_i, 3) = _vxc(_i, _i);
      _qp_energies_store(_i, 4) = _qp_energies(_i + _qpmin);
    }
  }

  // constructing full quasiparticle Hamiltonian and diagonalize, if requested
  if (_do_qp_diag || _do_bse_singlets || _do_bse_triplets) {
    FullQPHamiltonian();
    if (_do_qp_diag) {
      CTP_LOG(ctp::logDEBUG, *_pLog)
          << ctp::TimeStamp() << " Full quasiparticle Hamiltonian  " << flush;
      CTP_LOG(ctp::logINFO, *_pLog)
          << (format("  ====== Diagonalized quasiparticle energies (Hartree) "
                     "====== "))
                 .str()
          << flush;
      for (unsigned _i = 0; _i < _qptotal; _i++) {
        if ((_i + _qpmin) == _homo) {
          CTP_LOG(ctp::logINFO, *_pLog)
              << (format("  HOMO  = %1$4d PQP = %2$+1.4f DQP = %3$+1.4f ") %
                  (_i + _qpmin + 1) % _qp_energies(_i + _qpmin) %
                  _qp_diag_energies(_i))
                     .str()
              << flush;
        } else if ((_i + _qpmin) == _homo + 1) {
          CTP_LOG(ctp::logINFO, *_pLog)
              << (format("  LUMO  = %1$4d PQP = %2$+1.4f DQP = %3$+1.4f ") %
                  (_i + _qpmin + 1) % _qp_energies(_i + _qpmin) %
                  _qp_diag_energies(_i))
                     .str()
              << flush;

        } else {
          CTP_LOG(ctp::logINFO, *_pLog)
              << (format("  Level = %1$4d PQP = %2$+1.4f DQP = %3$+1.4f ") %
                  (_i + _qpmin + 1) % _qp_energies(_i + _qpmin) %
                  _qp_diag_energies(_i))
                     .str()
              << flush;
        }
      }

      // free memory

      if (!_store_qp_diag) {
        _qp_diag_coefficients.resize(0, 0);
        _qp_diag_energies.resize(0);
      }
    }  // _do_qp_diag
  }    // constructing full quasiparticle Hamiltonian

  // proceed only if BSE requested
  if (_do_bse_singlets || _do_bse_triplets) {

    // calculate direct part of eh interaction, needed for singlets and triplets
    BSE_d_setup(_Mmn);
    // add qp part to Eh_d
    BSE_Add_qp2H(_eh_d);
    // Only if verbose do we use this really
    if (votca::tools::globals::verbose) {
      BSE_qp_setup();
    }

    CTP_LOG(ctp::logDEBUG, *_pLog)
        << ctp::TimeStamp() << " Direct part of e-h interaction " << flush;

    if (_do_full_BSE) {
      BSE_d2_setup(_Mmn);
      CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                     << " Direct part of e-h interaction RARC "
                                     << flush;
    }

    if (_do_bse_triplets && _do_bse_diag) {
      BSE_solve_triplets();
      CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                     << " Solved BSE for triplets " << flush;
      // analyze and report results
      BSE_analyze_triplets();

    }  // do_triplets

    // constructing electron-hole interaction for BSE
    if (_do_bse_singlets) {
      // calculate exchange part of eh interaction, only needed for singlets
      BSE_x_setup(_Mmn);
      CTP_LOG(ctp::logDEBUG, *_pLog)
          << ctp::TimeStamp() << " Exchange part of e-h interaction " << flush;
    }

    if (_do_bse_singlets && _do_bse_diag) {

      if (_do_full_BSE) {
        BSE_solve_singlets_BTDA();
        CTP_LOG(ctp::logDEBUG, *_pLog)
            << ctp::TimeStamp() << " Solved full BSE for singlets " << flush;
        BSE_analyze_singlets_BTDA();
      } else {
        BSE_solve_singlets();
        CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                       << " Solved BSE for singlets " << flush;
        BSE_analyze_singlets();
      }

      if (!_store_eh_interaction) {
        _eh_d.resize(0, 0);
        _eh_x.resize(0, 0);
      }
    }
  }
  CTP_LOG(ctp::logDEBUG, *_pLog) << ctp::TimeStamp()
                                 << " GWBSE calculation finished " << flush;
  return true;
}
}
};

ext: cc

max_stars_repo_name: choudarykvsp/xtp

max_issues_repo_name: choudarykvsp/xtp

max_forks_repo_name: choudarykvsp/xtp


text:
#pragma once

#include <cstdint>
#include <memory>
#include <string_view>
#include <functional>

#include <boost/signals2.hpp>

#define GLFW_INCLUDE_NONE
#include <GLFW/glfw3.h>

#include "platform.hxx"

namespace platform {
class window final {
public:

    window(std::string_view name, std::int32_t width, std::int32_t height);

    ~window();

    void update(std::function<void(window &)> &&callback);

    GLFWwindow *handle() const noexcept { return handle_; }

    void connect_handler(std::shared_ptr<platform::event_handler> handler);
    void connect_handler(std::shared_ptr<platform::input_handler> handler);

private:
    GLFWwindow *handle_;

    std::int32_t width_{0}, height_{0};

    std::string name_;

    boost::signals2::signal<void(std::int32_t, std::int32_t)> resize_callback_;
    boost::signals2::signal<void(platform::input::raw_data &)> input_update_callback_;

    void set_callbacks();
};
}

ext: hxx

max_stars_repo_name: Alabuta/RayTracingNextWeek

max_issues_repo_name: Alabuta/RayTracingNextWeek

max_forks_repo_name: Alabuta/RayTracingNextWeek


text:
#pragma once

#include "beast.hpp"
#include "shared_state.hpp"
#include <boost/beast/version.hpp>

// Handles an HTTP server connection
class http_session : public std::enable_shared_from_this<http_session>
{
    beast::tcp_stream stream_;
    beast::flat_buffer buffer_;
    std::shared_ptr<shared_state> state_;
    // The parser is stored in an optional container so we can
    // construct it from scratch it at the beginning of each new message.
    boost::optional<http::request_parser<http::string_body>> parser_;

public:
    http_session(
        tcp::socket&& socket,
        std::shared_ptr<shared_state> const& state);
    ~http_session() {}

    void run();

private:
    void do_read();
    void on_read(beast::error_code ec, std::size_t);
    void send_payload(std::string const& payload);
    template <class BodyType>
    void write(http::response<BodyType>&& res);
    void on_write(beast::error_code ec, std::size_t, bool close);
    void do_close();

    template<class Body, class Allocator>
    void handle_request(http::request<Body, http::basic_fields<Allocator>>&& req);
};

max_stars_repo_name: Gongmeda/local-webrtc-signaling

max_issues_repo_name: Gongmeda/local-webrtc-signaling

max_forks_repo_name: Gongmeda/local-webrtc-signaling


text:
/** Test cases for Segment class.
 *
 */

#include "tmfsc/device.h"

#ifndef LINK_STATIC
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE DeviceTest
#include <boost/test/unit_test.hpp>

#include <iostream>

using namespace quest::tmfsc;
using namespace std;

Device create_simple_device (const point& A, const point& B, const point& C,
const point& D) {
    Device dev;
    int ptA = dev.addPoint(A);
    int ptB = dev.addPoint(B);
    int ptC = dev.addPoint(C);
    int ptD = dev.addPoint(D);

    dev.addEdge(ptA, ptB);
    dev.addEdge(ptB, ptC);
    dev.addEdge(ptC, ptD);
    dev.addEdge(ptD, ptA);

    return dev;
}

BOOST_AUTO_TEST_CASE(basics)
{
    Device dev = create_simple_device ({0, 0}, {20, 0}, {20, 10}, {0, 10});
    dev.edgeType(1, Edge::EDGE_ABSORB);

    BOOST_CHECK(dev.numEdges() == 4);
    BOOST_CHECK(dev.numConts() == 1);
    BOOST_CHECK(dev.isAbsorbEdge(1));
    BOOST_CHECK(dev.isReflectEdge(2));
    BOOST_CHECK(dev.contToEdgeIndx(0) == 1);
}

BOOST_AUTO_TEST_CASE(vectors)
{
    Device dev = create_simple_device ({0, 0}, {20, 0}, {20, 10}, {0, 10});
    dev.edgeType(1, Edge::EDGE_ABSORB);

    svec u = dev.edgeUnitVect(0);
    BOOST_CHECK_CLOSE(u(0), 1, 1E-4);
    BOOST_CHECK_CLOSE(u(1), 0, 1E-4);

    u = dev.edgeUnitVect(1);
    BOOST_CHECK_CLOSE(u(0), 0, 1E-4);
    BOOST_CHECK_CLOSE(u(1), 1, 1E-4);

    u = dev.edgeUnitVect(2);
    BOOST_CHECK_CLOSE(u(0), -1, 1E-4);
    BOOST_CHECK_CLOSE(u(1), 0, 1E-4);

    u = dev.edgeUnitVect(3);
    BOOST_CHECK_CLOSE(u(0), 0, 1E-4);
    BOOST_CHECK_CLOSE(u(1), -1, 1E-4);

    u = dev.edgeVect(3);
    BOOST_CHECK_CLOSE(u(0), 0, 1E-4);
    BOOST_CHECK_CLOSE(u(1), -10, 1E-4);
}

BOOST_AUTO_TEST_CASE(intersection)
{
    Device dev = create_simple_device ({0, 0}, {20, 0}, {20, 10}, {0, 10});
    dev.edgeType(1, Edge::EDGE_ABSORB);

    point ctr = {10, 5};
    point q = {30, 5};

    BOOST_CHECK(dev.intersects(ctr, q) == 1);
    point ints = dev.intersection(1, ctr, q);
    BOOST_CHECK_CLOSE(ints(0), 20, 1E-4);
    BOOST_CHECK_CLOSE(ints(1), 5, 1E-4);

    q(0) = 20; q(1) = 5;
    BOOST_CHECK(dev.intersects(ctr, q) == 1);

    q(0) = 20-0.0001; q(1) = 5;
    BOOST_CHECK(dev.intersects(ctr, q) == -1);

    point p1 = {20-0.01, 10-0.01};
    point q1 = {20-0.01, 10+0.01};
    BOOST_CHECK(dev.intersects(p1, q1) == 2);
    ints = dev.intersection(2, p1, q1);
    BOOST_CHECK_CLOSE(ints(0), 20-0.01, 1E-4);
    BOOST_CHECK_CLOSE(ints(1), 10, 1E-4);
}

BOOST_AUTO_TEST_CASE(pointsOnEdge)
{
    Device dev = create_simple_device ({0, 0}, {20, 0}, {20, 10}, {0, 10});
    dev.edgeType(1, Edge::EDGE_ABSORB);

    Edge e = dev.edge(1);
    cout << "Edge#0 p=(" << e.p()[0] << "," << e.p()[1] << ")"
                 << " q=(" << e.q()[0] << "," << e.q()[1] << ")"
                 << endl << endl;

    vector<point> pts = dev.createPointsOnCont(0, 11);
    for (int i = 0; i < pts.size(); i += 1) {
        cout << "pt#" << i << pts[i];
    }
}

max_stars_repo_name: masumhabib/quest

max_issues_repo_name: masumhabib/quest

max_forks_repo_name: masumhabib/quest