text:
/*
 * This file is open source software, licensed to you under the terms
 * of the Apache License, Version 2.0 (the "License").  See the NOTICE file
 * distributed with this work for additional information regarding copyright
 * ownership.  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.
 */

/*
 * Copyright (C) 2015 Cloudius Systems, Ltd.
 */

#pragma once

#include <seastar/core/thread_impl.hh>
#include <seastar/core/future.hh>
#include <seastar/core/do_with.hh>
#include <seastar/core/future-util.hh>
#include <seastar/core/timer.hh>
#include <seastar/core/reactor.hh>
#include <seastar/core/scheduling.hh>
#include <memory>
#include <setjmp.h>
#include <type_traits>
#include <chrono>
#include <seastar/util/std-compat.hh>
#include <ucontext.h>
#include <boost/intrusive/list.hpp>

/// \defgroup thread-module Seastar threads
///
/// Seastar threads provide an execution environment where blocking
/// is tolerated; you can issue I/O, and wait for it in the same function,
/// rather then establishing a callback to be called with \ref future<>::then().
///
/// Seastar threads are not the same as operating system threads:
///   - seastar threads are cooperative; they are never preempted except
///     at blocking points (see below)
///   - seastar threads always run on the same core they were launched on
///
/// Like other seastar code, seastar threads may not issue blocking system calls.
///
/// A seastar thread blocking point is any function that returns a \ref future<>.
/// you block by calling \ref future<>::get(); this waits for the future to become
/// available, and in the meanwhile, other seastar threads and seastar non-threaded
/// code may execute.
///
/// Example:
/// \code
///    seastar::thread th([] {
///       sleep(5s).get();  // blocking point
///    });
/// \endcode
///
/// An easy way to launch a thread and carry out some computation, and return a
/// result from this execution is by using the \ref seastar::async() function.
/// The result is returned as a future, so that non-threaded code can wait for
/// the thread to terminate and yield a result.

/// Seastar API namespace
namespace seastar {

/// \addtogroup thread-module
/// @{

class thread;
class thread_attributes;
class thread_scheduling_group;

/// Class that holds attributes controling the behavior of a thread.
class thread_attributes {
public:
    thread_scheduling_group* scheduling_group = nullptr;  // FIXME: remove
    compat::optional<seastar::scheduling_group> sched_group;
};

/// \cond internal
extern thread_local jmp_buf_link g_unthreaded_context;

// Internal class holding thread state.  We can't hold this in
// \c thread itself because \c thread is movable, and we want pointers
// to this state to be captured.
class thread_context {
    struct stack_deleter {
        void operator()(char *ptr) const noexcept;
    };
    using stack_holder = std::unique_ptr<char[], stack_deleter>;
    static constexpr size_t base_stack_size = 128*1024;

    thread_attributes _attr;
#ifdef SEASTAR_THREAD_STACK_GUARDS
    const size_t _stack_size;
#else
    static constexpr size_t _stack_size = base_stack_size;
#endif
    stack_holder _stack{make_stack()};
    noncopyable_function<void ()> _func;
    jmp_buf_link _context;
    scheduling_group _scheduling_group;
    promise<> _done;
    bool _joined = false;
    timer<> _sched_timer{[this] { reschedule(); }};
    compat::optional<promise<>> _sched_promise;

    boost::intrusive::list_member_hook<> _preempted_link;
    using preempted_thread_list = boost::intrusive::list<thread_context,
        boost::intrusive::member_hook<thread_context, boost::intrusive::list_member_hook<>,
        &thread_context::_preempted_link>,
        boost::intrusive::constant_time_size<false>>;

    boost::intrusive::list_member_hook<> _all_link;
    using all_thread_list = boost::intrusive::list<thread_context,
        boost::intrusive::member_hook<thread_context, boost::intrusive::list_member_hook<>,
        &thread_context::_all_link>,
        boost::intrusive::constant_time_size<false>>;

    static thread_local preempted_thread_list _preempted_threads;
    static thread_local all_thread_list _all_threads;
private:
    static void s_main(int lo, int hi); // all parameters MUST be 'int' for makecontext
    void setup();
    void main();
    stack_holder make_stack();
public:
    thread_context(thread_attributes attr, noncopyable_function<void ()> func);
    ~thread_context();
    void switch_in();
    void switch_out();
    bool should_yield() const;
    void reschedule();
    void yield();
    friend class thread;
    friend void thread_impl::switch_in(thread_context*);
    friend void thread_impl::switch_out(thread_context*);
    friend scheduling_group thread_impl::sched_group(const thread_context*);
};

/// \endcond

/// \brief thread - stateful thread of execution
///
/// Threads allow using seastar APIs in a blocking manner,
/// by calling future::get() on a non-ready future.  When
/// this happens, the thread is put to sleep until the future
/// becomes ready.
class thread {
    std::unique_ptr<thread_context> _context;
    static thread_local thread* _current;
public:
    /// \brief Constructs a \c thread object that does not represent a thread
    /// of execution.
    thread() = default;
    /// \brief Constructs a \c thread object that represents a thread of execution
    ///
    /// \param func Callable object to execute in thread.  The callable is
    ///             called immediately.
    template <typename Func>
    thread(Func func);
    /// \brief Constructs a \c thread object that represents a thread of execution
    ///
    /// \param attr Attributes describing the new thread.
    /// \param func Callable object to execute in thread.  The callable is
    ///             called immediately.
    template <typename Func>
    thread(thread_attributes attr, Func func);
    /// \brief Moves a thread object.
    thread(thread&& x) noexcept = default;
    /// \brief Move-assigns a thread object.
    thread& operator=(thread&& x) noexcept = default;
    /// \brief Destroys a \c thread object.
    ///
    /// The thread must not represent a running thread of execution (see join()).
    ~thread() { assert(!_context || _context->_joined); }
    /// \brief Waits for thread execution to terminate.
    ///
    /// Waits for thread execution to terminate, and marks the thread object as not
    /// representing a running thread of execution.
    future<> join();
    /// \brief Voluntarily defer execution of current thread.
    ///
    /// Gives other threads/fibers a chance to run on current CPU.
    /// The current thread will resume execution promptly.
    static void yield();
    /// \brief Checks whether this thread ought to call yield() now.
    ///
    /// Useful where we cannot call yield() immediately because we
    /// Need to take some cleanup action first.
    static bool should_yield();

    /// \brief Yield if this thread ought to call yield() now.
    ///
    /// Useful where a code does long running computation and does
    /// not want to hog cpu for more then its share
    static void maybe_yield();

    static bool running_in_thread() {
        return thread_impl::get() != nullptr;
    }
private:
    friend class reactor;
    // To be used by seastar reactor only.
    static bool try_run_one_yielded_thread();
};

/// An instance of this class can be used to assign a thread to a particular scheduling group.
/// Threads can share the same scheduling group if they hold a pointer to the same instance
/// of this class.
///
/// All threads that belongs to a scheduling group will have a time granularity defined by \c period,
/// and can specify a fraction \c usage of that period that indicates the maximum amount of time they
/// expect to run. \c usage, is expected to be a number between 0 and 1 for this to have any effect.
/// Numbers greater than 1 are allowed for simplicity, but they just have the same meaning of 1, alas,
/// "the whole period".
///
/// Note that this is not a preemptive runtime, and a thread will not exit the CPU unless it is scheduled out.
/// In that case, \c usage will not be enforced and the thread will simply run until it loses the CPU.
/// This can happen when a thread waits on a future that is not ready, or when it voluntarily call yield.
///
/// Unlike what happens for a thread that is not part of a scheduling group - which puts itself at the back
/// of the runqueue everytime it yields, a thread that is part of a scheduling group will only yield if
/// it has exhausted its \c usage at the call to yield. Therefore, threads in a schedule group can and
/// should yield often.
///
/// After those events, if the thread has already run for more than its fraction, it will be scheduled to
/// run again only after \c period completes, unless there are no other tasks to run (the system is
/// idle)
///
/// \deprecated This class has been replaced by \ref scheduling_group, which can be used to schedule
/// both threads and continuations.
class [[deprecated("Use seastar::scheduling_group instead")]] thread_scheduling_group {
    std::chrono::nanoseconds _period;
    std::chrono::nanoseconds _quota;
    std::chrono::time_point<thread_clock> _this_period_ends = {};
    std::chrono::time_point<thread_clock> _this_run_start = {};
    std::chrono::nanoseconds _this_period_remain = {};
public:
    /// \brief Constructs a \c thread_scheduling_group object
    ///
    /// \param period a duration representing the period
    /// \param usage which fraction of the \c period to assign for the scheduling group. Expected between 0 and 1.
    thread_scheduling_group(std::chrono::nanoseconds period, float usage);
    /// \brief changes the current maximum usage per period
    ///
    /// \param new_usage The new fraction of the \c period (Expected between 0 and 1) during which to run
    void update_usage(float new_usage) {
        _quota = std::chrono::duration_cast<std::chrono::nanoseconds>(new_usage * _period);
    }
private:
    void account_start();
    void account_stop();
    compat::optional<thread_clock::time_point> next_scheduling_point() const;
    friend class thread_context;
};

template <typename Func>
inline
thread::thread(thread_attributes attr, Func func)
        : _context(std::make_unique<thread_context>(std::move(attr), std::move(func))) {
}

template <typename Func>
inline
thread::thread(Func func)
        : thread(thread_attributes(), std::move(func)) {
}

inline
future<>
thread::join() {
    _context->_joined = true;
    return _context->_done.get_future();
}

/// Executes a callable in a seastar thread.
///
/// Runs a block of code in a threaded context,
/// which allows it to block (using \ref future::get()).  The
/// result of the callable is returned as a future.
///
/// \param attr a \ref thread_attributes instance
/// \param func a callable to be executed in a thread
/// \param args a parameter pack to be forwarded to \c func.
/// \return whatever \c func returns, as a future.
///
/// Example:
/// \code
///    future<int> compute_sum(int a, int b) {
///        thread_attributes attr = {};
///        attr.scheduling_group = some_scheduling_group_ptr;
///        return seastar::async(attr, [a, b] {
///            // some blocking code:
///            sleep(1s).get();
///            return a + b;
///        });
///    }
/// \endcode
template <typename Func, typename... Args>
inline
futurize_t<std::result_of_t<std::decay_t<Func>(std::decay_t<Args>...)>>
async(thread_attributes attr, Func&& func, Args&&... args) {
    using return_type = std::result_of_t<std::decay_t<Func>(std::decay_t<Args>...)>;
    struct work {
        thread_attributes attr;
        Func func;
        std::tuple<Args...> args;
        promise<return_type> pr;
        thread th;
    };
    return do_with(work{std::move(attr), std::forward<Func>(func), std::forward_as_tuple(std::forward<Args>(args)...)}, [] (work& w) mutable {
        auto ret = w.pr.get_future();
        w.th = thread(std::move(w.attr), [&w] {
            futurize<return_type>::apply(std::move(w.func), std::move(w.args)).forward_to(std::move(w.pr));
        });
        return w.th.join().then([ret = std::move(ret)] () mutable {
            return std::move(ret);
        });
    });
}

/// Executes a callable in a seastar thread.
///
/// Runs a block of code in a threaded context,
/// which allows it to block (using \ref future::get()).  The
/// result of the callable is returned as a future.
///
/// \param func a callable to be executed in a thread
/// \param args a parameter pack to be forwarded to \c func.
/// \return whatever \c func returns, as a future.
template <typename Func, typename... Args>
inline
futurize_t<std::result_of_t<std::decay_t<Func>(std::decay_t<Args>...)>>
async(Func&& func, Args&&... args) {
    return async(thread_attributes{}, std::forward<Func>(func), std::forward<Args>(args)...);
}
/// @}

}

ext: hh

max_stars_repo_name: elazarl/seastar

max_issues_repo_name: elazarl/seastar

max_forks_repo_name: elazarl/seastar


text:
/*
Author: Gudakov Ramil Sergeevich a.k.a. Gauss
Гудаков Рамиль Сергеевич
Contacts: [ramil2085@mail.ru, ramil2085@gmail.com]
See for more information LICENSE.md.
*/

#include <boost/asio/ip/impl/address_v4.ipp>
#include <boost/lexical_cast.hpp>

#include "HandlerMMO_Master.h"
#include "MMOEngine/include/BaseEvent.h"
#include "MMOEngine/include/Events.h"
#include "MMOEngine/include/Master.h"
#include "CommonParam.h"
#include "ResolverSelf_IP_v4.h"
#include "Logger.h"

using namespace std;

THandlerMMO_Master::THandlerMMO_Master( nsMMOEngine::TBase* pBase ) : THandlerMMO( pBase, eMaster )
{
  std::string slaveLogin = SLAVE_LOGIN;
  mSHA256.FastCalc( (void*)slaveLogin.data(), slaveLogin.length(), mSlaveHashLogin );
}
//-----------------------------------------------------------------------------------
void THandlerMMO_Master::HandleFromMMOEngine( nsEvent::TEvent* pEvent )
{
  nsMMOEngine::TBaseEvent* pBE = ( nsMMOEngine::TBaseEvent* )pEvent->pContainer->GetPtr();
  nsMMOEngine::TMaster* pMaster = ( nsMMOEngine::TMaster* )pEvent->pSrc;

  string sEvent;
  switch ( pBE->mType )
  {
    case nsMMOEngine::eTryConnectDown:
    {
      auto pTryConnectDown = ( nsMMOEngine::TTryConnectDownEvent* ) pBE;
      sEvent = "TryConnectDown";
      auto sessionID = pTryConnectDown->sessionID;
      std::string password = CLIENT_PASSWORD;

      if ( memcmp( mSlaveHashLogin.GetPtr(), pTryConnectDown->hashLogin, pTryConnectDown->hashLoginSize ) == 0 )
        password = SLAVE_PASSWORD;
      else
        AddTryConnectClientToMaster( sessionID );

      pMaster->Accept( sessionID, password );
    }
    break;
    case nsMMOEngine::eConnectDown:
      sEvent = "ConnectDown";
      AddConnection( ( ( nsMMOEngine::TConnectDownEvent* ) pBE )->sessionID );
      break;
    case nsMMOEngine::eDisconnectDown:
      sEvent = "DisconnectDown";
      RemoveConnection( ( ( nsMMOEngine::TDisconnectDownEvent* ) pBE )->sessionID );
      break;
    case nsMMOEngine::eConnectUp:
      sEvent = "ConnectUp";
      AddConnection( ( ( nsMMOEngine::TConnectUpEvent* ) pBE )->sessionID );
      break;
    case nsMMOEngine::eDisconnectUp:
      sEvent = "DisconnectUp";
      RemoveConnection( ( ( nsMMOEngine::TDisconnectUpEvent* ) pBE )->sessionID );
      break;
    case nsMMOEngine::eError:
    {
      nsMMOEngine::TErrorEvent* pEr = ( nsMMOEngine::TErrorEvent* )pBE;
      sEvent = nsMMOEngine::GetStrError( pEr->code );
    }
    break;
    case nsMMOEngine::eRecvFromDown:
    {
      sEvent = "RecvFromDown";
      //nsMMOEngine::TRecvFromDownEvent* pR = (nsMMOEngine::TRecvFromDownEvent*)pBE;
      //char s[200];
      //memcpy( s, pR->data, pR->dataSize );
      //s[pR->dataSize] = '\0';
      //sEvent += " msg: ";
      //sEvent += s;
      //// получили пакет от Slave с clientKey
      //std::list<unsigned int> l_id;
      //l_id.push_front( boost::lexical_cast<int>(s) );
      //char* sMsgFromMaster = "Master say hello!";
      //int sizeMsg = strlen( sMsgFromMaster );
      //pMaster->SendByClientKey( l_id, sMsgFromMaster, sizeMsg );
    }
    break;
    case nsMMOEngine::eRecvFromUp:
    {
      sEvent = "RecvFromUp";
      nsMMOEngine::TRecvFromUpEvent* pR = ( nsMMOEngine::TRecvFromUpEvent* )pBE;
      char s[200];
      memcpy( s, pR->GetData(), pR->GetSize() );
      s[pR->GetSize()] = '\0';
      sEvent += " msg: ";
      sEvent += s;
    }
    break;
    case nsMMOEngine::eTryLogin:
    {
      sEvent = "TryLogin";
      nsMMOEngine::TTryLoginEvent* pETL = ( nsMMOEngine::TTryLoginEvent* )pBE;
      char lenLogin = *( pETL->c.GetPtr() );
      char s[100];
      memcpy( s, pETL->c.GetPtr() + 1, lenLogin );
      s[lenLogin] = '\0';
      unsigned int clientKey = boost::lexical_cast<int>( s );
      pMaster->SetResultLogin( true, pETL->sessionID, clientKey, ( void* )"Welcome", strlen( "Welcome" ) );
    }
    break;
    case nsMMOEngine::eLogin:
    {
      sEvent = "Login";
      auto pLoginEvent = ( nsMMOEngine::TLoginEvent* )pBE;
      unsigned int clientKey = pLoginEvent->clientKey;
      AddClient( clientKey );
    }
    break;
    case nsMMOEngine::eLogoff:
    {
      nsMMOEngine::TLogoffEvent* pLogoff = ( nsMMOEngine::TLogoffEvent* )pBE;
      char s[100];
      sprintf( s, "%u", pLogoff->clientKey );
      sEvent = "Logoff ";
      sEvent += s;
      RemoveClient( pLogoff->clientKey );
    }
    break;
    case nsMMOEngine::eDestroyGroup:
      sEvent = "DestroyGroup";
      break;
    default:BL_FIX_BUG();
  }

  if ( pBE->mType == nsMMOEngine::eError )
  {
    GetLogger( ServerLog )->WriteF( "MMOEngine M (0x%p): %s.\t", pMaster, sEvent.data() );
    PrintCC( ServerLog );
  }
}
//---------------------------------------------------------------------------------------------

ext: cpp

max_stars_repo_name: RamilGauss/MMO-Framework

max_issues_repo_name: RamilGauss/MMO-Framework

max_forks_repo_name: RamilGauss/MMO-Framework


text:
#include "DirectedStringGraph.h"
#include <fstream>
#include <ostream>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include "BidirectedStringGraph.h"
#include <math.h>
#include <lemon/network_simplex.h>
#include <lemon/smart_graph.h>

const char DirectedStringGraph::magic[] =
	{'D', 'i', 'g', 'r', 'a', 'p', 'h', '\0', '\0', '\0'};

//
// Perform transitive edge reduction on a directed string graph:
// Remove all edges v -> x where there exist edges v -> w -> x.
//
// See:
//    "The fragment string assembly graph", Eugene W. Myers, 2005
//
// There are, however, some changes that have to be made to make the algorithm
// work in practice.  In particular, the edge v -> x should only be removed if
// it is actually labeled by the same sequence as v -> w -> x.
//
// It's important to keep in mind that there can only at most one edge v -> w
// for any v and w, since this represents one kind of overlap.  But, it's
// possible for there to be an edge w -> v as well.
//
void DirectedStringGraph::transitive_reduction()
{
	info("Performing transitive reduction on directed string graph with "
	     "%zu vertices and %zu edges", this->num_vertices(), this->num_edges());

	std::vector<DirectedStringGraphVertex> & vertices = this->vertices();
	std::vector<DirectedStringGraphEdge> & edges = this->edges();

	info("Sorting adjacency lists of vertices by edge length");
	this->sort_adjlists_by_edge_len();

	static const unsigned char VACANT = 0;
	static const unsigned char INPLAY = 1;
	static const unsigned char ELIMINATED = 2;

	// Per-vertex mark.  Initially set to VACANT.
	std::vector<unsigned char> vertex_marks(this->num_vertices(), VACANT);

	// Map from vertex index to index of edges back to the vertex currently
	// under consideration
	std::vector<edge_idx_t> v_idx_to_back_edge_idx(this->num_vertices(),
						       std::numeric_limits<edge_idx_t>::max());

	// Per-edge boolean indicating whether each edge is to be deleted as
	// part of the transitive reduction or not.
	std::vector<bool> reduce_edge(this->num_edges(), false);

	info("Looking for transitive edges based at each of %zu vertices",
	     vertices.size());

	// Iterate through every vertex @v in the graph that has outgoing edges.
	for (size_t v_idx = 0; v_idx < vertices.size(); v_idx++) {
		const DirectedStringGraphVertex & v = vertices[v_idx];

		if (v.out_degree() == 0)
			continue;

		// Mark each vertex adjacent to @v as INPLAY, and initialize the
		// map from the adjacent vertices' indices to the back edges'
		// indices.
		foreach(const edge_idx_t edge_idx, v.edge_indices()) {
			const DirectedStringGraphEdge & e = edges[edge_idx];
			const v_idx_t w_idx = e.get_v2_idx();
			assert(edge_idx != std::numeric_limits<edge_idx_t>::max());
			assert(v_idx_to_back_edge_idx[w_idx] ==
			       std::numeric_limits<edge_idx_t>::max());
			v_idx_to_back_edge_idx[w_idx] = edge_idx;

			vertex_marks[w_idx] = INPLAY;
		}

		// Length of the longest sequence label on the edges leaving
		// vertex @v.
		const BaseVec::size_type longest = edges[v.edge_indices().back()].length();

		// For each outgoing edge from v -> w in order of labeled
		// sequence length, consider each vertex w that is still marked
		// INPLAY.
		foreach(const edge_idx_t edge_idx, v.edge_indices()) {
			const DirectedStringGraphEdge & e = edges[edge_idx];
			const v_idx_t w_idx = e.get_v2_idx();

			if (vertex_marks[w_idx] != INPLAY)
				continue;

			// The edge v -> w must be an irreducible edge if w is
			// still marked INPLAY at this point, since all shorter
			// edges were already considered.
			//
			// Now, consider the edges leaving vertex w.  Each such
			// edge that goes to a vertex marked INPLAY must be
			// directly reachable from v, and therefore the edge
			// must be removed, unless its sequence does not
			// actually match the sequence from the edges
			// v -> w -> x.
			const DirectedStringGraphVertex & w = vertices[w_idx];
			foreach(const edge_idx_t w_edge_idx, w.edge_indices()) {
				const DirectedStringGraphEdge & e2 = edges[w_edge_idx];
				if (e.length() + e2.length() > longest)
					break;

				const v_idx_t x_idx = e2.get_v2_idx();
				if (vertex_marks[x_idx] != INPLAY)
					continue;

				const edge_idx_t back_edge_idx =
						v_idx_to_back_edge_idx[x_idx];

				assert(back_edge_idx !=
				       std::numeric_limits<edge_idx_t>::max());

				const DirectedStringGraphEdge & back_edge =
					edges[back_edge_idx];

				assert(back_edge.get_v1_idx() == v_idx);
				assert(back_edge.get_v2_idx() == x_idx);

				if (e.length() + e2.length() != back_edge.length())
					continue;

				for (BaseVec::size_type i = 0; i < e.length(); i++)
					if (e.get_seq()[i] != back_edge.get_seq()[i])
						goto next_edge;

				for (BaseVec::size_type i = 0; i < e2.length(); i++)
					if (e2.get_seq()[i] !=
					    back_edge.get_seq()[i + e.length()])
						goto next_edge;

				vertex_marks[x_idx] = ELIMINATED;
				next_edge:
				;
			}
		}

		// Once again, go through the outgoing edges from v.  For each
		// neighboring vertex marked ELIMINATED, mark the corresponding
		// edge(s) for reduction.  Return both INPLAY and ELIMINATED
		// vertices to VACANT status.
		foreach(const edge_idx_t edge_idx, v.edge_indices()) {
			const DirectedStringGraphEdge & e = edges[edge_idx];
			const v_idx_t w_idx = e.get_v2_idx();
			if (vertex_marks[w_idx] == ELIMINATED)
				reduce_edge[edge_idx] = true;
			v_idx_to_back_edge_idx[w_idx] = std::numeric_limits<edge_idx_t>::max();
			vertex_marks[w_idx] = VACANT;
		}
	}

	info("Transitive reduction algorithm complete.  Now updating the string graph");

	for (size_t i = 0; i < edges.size(); i++) {
		if (reduce_edge[i]) {
			const DirectedStringGraphEdge & e = _edges[i];
			const edge_idx_t j = locate_edge(e.get_v2_idx() ^ 1,
							 e.get_v1_idx() ^ 1);
			if (!reduce_edge[j]) {
				std::cout << "The following 2 edges are opposites "
					<< "but were not both reduced:" << std::endl;
				e.print(std::cout, 0, true);
				std::cout << std::endl;
				_edges[j].print(std::cout, 0, true);
				std::cout << std::endl;
			}
		}
	}

	// Update the directed string graph to remove the edges marked %true in
	// the @reduce_edge array.

	// Map from the old edge indices to the new edge indices.
	std::vector<edge_idx_t> new_edge_indices(edges.size());
	size_t i, j;
	for (i = 0, j = 0; i < edges.size(); i++) {
		if (reduce_edge[i]) {
			new_edge_indices[i] = std::numeric_limits<edge_idx_t>::max();
		} else {
			new_edge_indices[i] = j;
			edges[j++] = edges[i];
		}
	}
	const size_t num_original_edges = edges.size();
	const size_t num_removed_edges = num_original_edges - j;
	const size_t num_remaining_edges = j;
	edges.resize(num_remaining_edges);

	info("Removing %zu of %zu edges (%.2f%%)",
	     num_removed_edges, num_original_edges,
	     TO_PERCENT(num_removed_edges, num_original_edges));

	// Re-number the edge indices list of each vertex, and remove any
	// indices that correspond to edges that were removed.
	foreach(DirectedStringGraphVertex & v, vertices) {
		std::vector<edge_idx_t> & edge_indices = v.edge_indices();
		size_t j = 0;
		foreach(const edge_idx_t edge_idx, edge_indices) {
			const edge_idx_t new_edge_idx = new_edge_indices[edge_idx];
			if (new_edge_idx != std::numeric_limits<edge_idx_t>::max()) {
				edge_indices[j++] = new_edge_idx;
			}
		}
		edge_indices.resize(j);
	}

	info("Done removing transitive edges");
}

void DirectedStringGraph::follow_unbranched_path(DirectedStringGraphEdge & e,
						 std::vector<bool> & remove_edge,
						 std::vector<bool> & remove_vertex,
						 const std::vector<bool> & v_inner)
{
	BaseVecVec::size_type new_seq_len = e.length();
	v_idx_t vi_idx = e.get_v2_idx();
	assert(v_inner[vi_idx]);

	v_idx_t num_inner_vertices = 0;
	// Found beginning of unbranched path.  Walk along it until
	// the end to get the total sequence length.
	do {
		const DirectedStringGraphVertex &vi = _vertices[vi_idx];
		assert(vi.out_degree() == 1);
		const DirectedStringGraphEdge &ei_i1 = _edges[vi.first_edge_idx()];
		if (new_seq_len + ei_i1.length() < new_seq_len)
			fatal_error("Edge too long");
		new_seq_len += ei_i1.length();
		num_inner_vertices++;
		vi_idx = ei_i1.get_v2_idx();
	} while (v_inner[vi_idx]);

	// Now build the new sequence, update the first edge in the path to skip
	// to the last vertex in the path, and mark the other edges for removal
	// in the @remove_edge array.
	BaseVec & new_seq = e.get_seq();
	BaseVec::size_type seq_idx = new_seq.length();
	new_seq.resize(new_seq_len);
	vi_idx = e.get_v2_idx();

	e.set_num_inner_vertices(num_inner_vertices);
	while (num_inner_vertices--) {
		const DirectedStringGraphVertex &vi = _vertices[vi_idx];
		const edge_idx_t ei_i1_idx = vi.first_edge_idx();
		const DirectedStringGraphEdge &ei_i1 = _edges[ei_i1_idx];
		const BaseVec & ei_i1_seq = ei_i1.get_seq();
		for (BaseVec::size_type i = 0; i < ei_i1_seq.length(); i++) {
			assert2(seq_idx < new_seq_len);
			new_seq.set(seq_idx++, ei_i1_seq[i]);
		}
		e.increment_mapped_read_count(ei_i1.get_mapped_read_count());
		remove_edge[ei_i1_idx] = true;
		remove_vertex[vi_idx] = true;
		vi_idx = ei_i1.get_v2_idx();
	}
	assert(!v_inner[vi_idx]);
	assert(new_seq_len == seq_idx);
	e.set_v2_idx(vi_idx);
}

void DirectedStringGraph::collapse_unbranched_paths()
{
	const v_idx_t n_verts = num_vertices();
	const edge_idx_t n_edges = num_edges();

	info("Collapsing unbranched paths in directed string graph");
	info("Original graph has %zu vertices and %zu edges", n_verts, n_edges);

	size_t num_inner_vertices = 0;
	// Find whether each vertex is inner or not.  A vertex is inner iff it
	// has indegree 1 and outdegree 1.
	std::vector<bool> v_inner(n_verts, false);
	{
		std::vector<unsigned char> v_in_degrees(n_verts, 0);
		std::vector<unsigned char> v_out_degrees(n_verts, 0);

		foreach (const DirectedStringGraphEdge & e, _edges) {
			v_idx_t v1_idx, v2_idx;
			e.get_v_indices(v1_idx, v2_idx);
			if (v_out_degrees[v1_idx] < 2)
				v_out_degrees[v1_idx]++;
			if (v_in_degrees[v2_idx] < 2)
				v_in_degrees[v2_idx]++;
		}
		for (v_idx_t v_idx = 0; v_idx < n_verts; v_idx++) {
			if (v_in_degrees[v_idx] == 1 && v_out_degrees[v_idx] == 1) {
				v_inner[v_idx] = true;
				num_inner_vertices++;
			}
		}
	}

	info("Found %zu inner vertices (%.2f%% of all vertices)",
	     num_inner_vertices, TO_PERCENT(num_inner_vertices, n_verts));

	// Go through each non-inner vertex and look for any neighboring inner
	// vertices.  These are the starts of unbranched paths that will be
	// collapsed.  For each such path, follow it until its end and determine
	// the total length of the edges.  Then, set the first edge's label to
	// the concatenation of the labels of all the edges in the unbranched
	// path, then make the first edge point directly to the last vertex in
	// the path and mark the other edges for deletion.
	//
	// Note: smooth rings are not collapsed yet.
	//
	size_t num_unbranched_paths = 0;
	std::vector<bool> remove_edge(n_edges, false);
	std::vector<bool> remove_vertex(n_verts, false);
	for (v_idx_t v_idx = 0; v_idx < n_verts; v_idx++) {
		if (!v_inner[v_idx]) {
			const DirectedStringGraphVertex & v = _vertices[v_idx];
			foreach(edge_idx_t edge_idx, v.edge_indices()) {
				DirectedStringGraphEdge & e = _edges[edge_idx];
				v_idx_t v2_idx = e.get_v2_idx();
				if (v_inner[v2_idx] && !remove_vertex[v2_idx]) {
					num_unbranched_paths++;
					follow_unbranched_path(e, remove_edge,
							       remove_vertex, v_inner);
				}
			}
		}
	}

	for (v_idx_t v_idx = 0; v_idx < n_verts; v_idx++) {
		if (v_inner[v_idx] && !remove_vertex[v_idx]) {
			std::cerr << "Graph contains a smooth ring!" << std::endl;
			unimplemented();
		}
	}

	info("Found %zu unbranched paths", num_unbranched_paths);

	// Compute the new vertex indices.
	std::vector<v_idx_t> old_to_new_v_indices(n_verts,
						  std::numeric_limits<v_idx_t>::max());
	v_idx_t new_v_idx = 0;
	for (v_idx_t old_v_idx = 0; old_v_idx < n_verts; old_v_idx++)
		if (!remove_vertex[old_v_idx])
			old_to_new_v_indices[old_v_idx] = new_v_idx++;

	info("Updated vertices are indexed [0, %lu)", new_v_idx);

	info("Updating edges");
	// Compute the new edge indices, set the new vertex indices in each edge,
	// and move the edges
	std::vector<edge_idx_t> old_to_new_edge_indices(n_edges,
							std::numeric_limits<edge_idx_t>::max());
	edge_idx_t new_edge_idx = 0;
	for (edge_idx_t old_edge_idx = 0; old_edge_idx < n_edges; old_edge_idx++) {
		if (!remove_edge[old_edge_idx]) {
			DirectedStringGraphEdge & e = _edges[old_edge_idx];
			old_to_new_edge_indices[old_edge_idx] = new_edge_idx;
			v_idx_t v1_idx, v2_idx;
			e.get_v_indices(v1_idx, v2_idx);
			e.set_v_indices(old_to_new_v_indices[v1_idx],
					old_to_new_v_indices[v2_idx]);
			assert2(e.get_v1_idx() != std::numeric_limits<v_idx_t>::max() &&
				e.get_v2_idx() != std::numeric_limits<v_idx_t>::max());
			_edges[new_edge_idx++] = _edges[old_edge_idx];
		}
	}
	info("Updated edges are indexed [0, %lu)", new_edge_idx);
	info("%zu edges were removed (%f%% of total)",
	     _edges.size() - new_edge_idx,
	     TO_PERCENT(_edges.size() - new_edge_idx, _edges.size()));

	_edges.resize(new_edge_idx);

	info("Updating vertices");
	// Set new edge indices in each vertex and move the vertices
	new_v_idx = 0;
	for (v_idx_t old_v_idx = 0; old_v_idx < n_verts; old_v_idx++) {
		if (!v_inner[old_v_idx]) {
			DirectedStringGraphVertex & v = _vertices[old_v_idx];
			std::vector<edge_idx_t> & edge_indices = v.edge_indices();
			for (size_t i = 0; i < v.out_degree(); i++) {
				edge_indices[i] = old_to_new_edge_indices[edge_indices[i]];
				assert2(edge_indices[i] != std::numeric_limits<edge_idx_t>::max());
			}
			_vertices[new_v_idx++] = _vertices[old_v_idx];
		}
	}
	assert(new_v_idx == n_verts - num_inner_vertices);
	_vertices.resize(new_v_idx);
	info("Done collapsing unbranched paths in directed string graph");
}

void DirectedStringGraph::mark_component(const v_idx_t v_idx,
					 std::vector<bool> & visited,
					 v_idx_t & component_size) const
{
	assert(!visited[v_idx]);
	visited[v_idx] = true;
	component_size++;
	foreach (edge_idx_t edge_idx, _vertices[v_idx].edge_indices()) {
		const v_idx_t w_idx = _edges[edge_idx].get_v2_idx();
		if (!visited[w_idx]) {
			mark_component(w_idx, visited, component_size);
		}
	}
}

void DirectedStringGraph::print_stats(std::ostream & os) const
{
	os << "DirectedStringGraph {" << std::endl;
	os << "    Number of vertices: " << num_vertices() << std::endl;
	os << "    Number of edges: " << num_edges() << std::endl;

	std::vector<unsigned char> out_degrees(num_vertices(), 0);
	std::vector<unsigned char> in_degrees(num_vertices(), 0);
	double total_mapped_count = 0.0;
	size_t more_than_one_mapped_count = 0;
	foreach(const DirectedStringGraphEdge & e, _edges) {
		if (out_degrees[e.get_v1_idx()] + 1 != 0) {
			out_degrees[e.get_v1_idx()]++;
		}
		if (in_degrees[e.get_v2_idx()] + 1 != 0) {
			in_degrees[e.get_v2_idx()]++;
		}
		assert(e.get_mapped_read_count() >= 1.0);
		total_mapped_count += e.get_mapped_read_count();
		if (e.get_mapped_read_count() > 1.0)
			more_than_one_mapped_count++;
	}
	std::vector<v_idx_t> out_degree_hist(0xff, 0);
	std::vector<v_idx_t> in_degree_hist(0xff, 0);
	std::vector<v_idx_t> in_out_degree_hist(0xffff, 0);
	v_idx_t v_in_neq_out = 0;
	for (v_idx_t v_idx = 0; v_idx < num_vertices(); v_idx++) {
		out_degree_hist[out_degrees[v_idx]]++;
		in_degree_hist[in_degrees[v_idx]]++;
		if (out_degrees[v_idx] != in_degrees[v_idx])
			v_in_neq_out++;
		in_out_degree_hist[(v_idx_t(in_degrees[v_idx]) << 8) + out_degrees[v_idx]]++;
	}
	v_idx_t max_out_degree = 0, max_in_degree = 0;
	for (size_t i = 0; i < out_degree_hist.size(); i++)
		if (out_degree_hist[i] != 0)
			max_out_degree = i;
	for (size_t i = 0; i < in_degree_hist.size(); i++)
		if (in_degree_hist[i] != 0)
			max_in_degree = i;
	os << "    Number of isolated vertices: "
	   << in_out_degree_hist[0x000] << std::endl;
	os << "    Number of inner vertices: "
	   << in_out_degree_hist[0x101] << std::endl;
	os << "    Number of branch beginning vertices: "
	   << in_out_degree_hist[0x001] << std::endl;
	os << "    Number of branch ending vertices: "
	   << in_out_degree_hist[0x100] << std::endl;
	os << "    Number of vertices with unequal in degree and out degree: "
	   << v_in_neq_out << std::endl;
	os << "    Max in degree: " << max_in_degree
	   << (max_in_degree == 0xff ? '+' : ' ') << std::endl;
	os << "    Max out degree: " << max_out_degree
	   << (max_in_degree == 0xff ? '+' : ' ') << std::endl;

	os << "    Number of edges that one or more contained reads map onto: "
	   << more_than_one_mapped_count << std::endl;
	os << "    Average number of contained reads that map onto each edge: "
	   << (num_edges() ? total_mapped_count / num_edges() : 0) << std::endl;

	std::vector<bool> visited(num_vertices(), false);
	std::vector<v_idx_t> component_sizes;

	for (v_idx_t v_idx = 0; v_idx < num_vertices(); v_idx++) {
		if (!visited[v_idx] && in_degrees[v_idx] == 0) {
			v_idx_t component_size = 0;
			mark_component(v_idx, visited, component_size);
			component_sizes.push_back(component_size);
		}
	}
	for (v_idx_t v_idx = 0; v_idx < num_vertices(); v_idx++) {
		if (!visited[v_idx]) {
			v_idx_t component_size = 0;
			mark_component(v_idx, visited, component_size);
			component_sizes.push_back(component_size);
		}
	}
	std::sort(component_sizes.begin(), component_sizes.end());
	os << "    Number of components: " << component_sizes.size() << std::endl;
	os << "    Component sizes:" << std::endl;
	foreach (v_idx_t component_size, component_sizes) {
		os << "        " << component_size << ' '
		   << (component_size == 1 ? "vertex" : "vertices") << std::endl;
	}
	os << "}" << std::endl;
}

//
// Consider a read f contained in another read g:
//
//
//                  f.B ---------> f.E
//          g.B ------------------------> g.E
//
// where the relevant portion of the actual graph might look like:
//
// 	v1.? -> v2.? -> v3.? -> g.E
//
// Or:
//
//                  f.B ---------> f.E
//          g.E <------------------------ g.B
//
//  where the relevant portion of the actual graph might look like:
//
//            v1.? -> v2.? -> v3.? -> g.B
//
// In reality, there are no vertices f.B and f.E because they were removed from
// the graph.  The goal here is to start at vertex g.E (or g.B) and walk
// backwards (e.g. towards the v?.? shown above) to find all possible edges into
// which the f.E vertex would be located, so that the number of reads that map
// onto those edges can be incremented.  Ideally, there would be only 1 such
// edge--- but due to branching, there may be many.  At most 100 mapped edges
// are considered before throwing away this contained read.
//
// @v_idx:
// 	Next branchpoint vertex.  (initially: g.E for forward case; g.B for
// 	reverse-complement case).
//
// @overhang_len:
// 	Number of bp of uncontained read remaining after threading the contained
// 	read all the way to the vertex of index @v_idx.
//
// @mapped_edges:
// 	Array of indices of edges onto which the contained read would map.
//
// @num_mapped_edges:
// 	Number of filled members of @mapped_edges.
//
// @max_mapped_edges:
// 	Size of max_mapped_edges.  The search is terminated when
// 	@num_mapped_edges reaches @max_mapped_edges.
//
void DirectedStringGraph::walk_back_edges(const v_idx_t v_idx,
					  size_t overhang_len,
					  edge_idx_t mapped_edges[],
					  size_t & num_mapped_edges,
					  size_t max_mapped_edges)
{
	assert(_back_edges.size() == _vertices.size());

	// Map the read onto any edges that go into the current vertex with
	// length less than the remaining read length.
	foreach (edge_idx_t edge_idx, _back_edges[v_idx]) {
		const DirectedStringGraphEdge & e = _edges[edge_idx];
		assert(e.get_v2_idx() == v_idx);
		if (overhang_len < e.length()) {
			mapped_edges[num_mapped_edges++] = edge_idx;
			if (num_mapped_edges == max_mapped_edges)
				return;
		}
	}

	// For each edge going into the current vertex with length greater than
	// or equal to the remaining read length, recursively call this
	// procedure on them with the remaining read length decremented by the
	// legnth of the edge.
	//
	// If upon return of the recursive call, the number of mapped edges has
	// reached the maximum, return early.
	foreach (edge_idx_t edge_idx, _back_edges[v_idx]) {
		const DirectedStringGraphEdge & e = _edges[edge_idx];
		if (overhang_len >= e.length()) {
			walk_back_edges(e.get_v1_idx(),
					overhang_len - e.length(),
					mapped_edges,
					num_mapped_edges,
					max_mapped_edges);
			if (num_mapped_edges == max_mapped_edges)
				return;
		}
	}
}

//
//
void DirectedStringGraph::map_contained_read(const v_idx_t downstream_read_idx,
					     const v_idx_t downstream_read_dir,
					     const BaseVec::size_type overhang_len)
{
	assert(downstream_read_idx < num_vertices() / 2);
	assert(downstream_read_dir < 2);

	const v_idx_t downstream_v_idx = downstream_read_idx * 2 + downstream_read_dir;

	// We need to be able to walk the graph in the opposite direction that
	// the edges are going--- so index the back edges, so that it's possible
	// to get a list of edges entering a vertex, given the vertex index.
	if (_back_edges.size() == 0) {
		info("Indexing back edges (num_vertices = %zu, num_edges = %zu)",
				num_vertices(), num_edges());
		_back_edges.resize(num_vertices());
		foreach (const DirectedStringGraphVertex & v, _vertices)
			foreach (const edge_idx_t edge_idx, v.edge_indices())
				_back_edges[_edges[edge_idx].get_v2_idx()].push_back(edge_idx);
	}

	// Walk the tree backwards starting at vertex @downstream_v_idx to get a
	// list of edge indices into which the end of the contained read maps.
	edge_idx_t mapped_edges[100];
	size_t num_mapped_edges = 0;
	walk_back_edges(downstream_v_idx, overhang_len,
			mapped_edges, num_mapped_edges, 100);

	// If there was at least one mapped edge but the maximum was not
	// reached, increment the mapped read count on the corresponding edges,
	// but divide by the number of mapped locations to weight the mapping
	// appropriately.
	if (num_mapped_edges < 100 && num_mapped_edges > 0) {
		double div = 1.0f / double(num_mapped_edges);
		for (size_t i = 0; i < num_mapped_edges; i++) {
			assert(mapped_edges[i] < num_edges());
			_edges[mapped_edges[i]].increment_mapped_read_count(div);
		}
	}
}

//
// Given a bidirected string graph, build the corresponding directed string
// graph.
//
void DirectedStringGraph::build_from_bidigraph(const BidirectedStringGraph & bidigraph)
{
	assert(num_vertices() == bidigraph.num_vertices() * 2);

	// Go through each bidirected edge in the bidirected string graph and
	// add the corresponding 2 directed edges to the directed string graph.
	foreach (const BidirectedStringGraphEdge & e, bidigraph.edges()) {
		add_edge_pair(e.get_v1_idx(),
			      e.get_v2_idx(),
			      e.get_dirs(),
			      e.get_seq_1_to_2(), 0, e.length() - 1, false,
			      e.get_seq_2_to_1(), 0, e.length() - 1, false);
	}
	_orig_num_reads = bidigraph._orig_num_reads;
}

void DirectedStringGraph::calculate_A_statistics()
{
	size_t num_reads = _orig_num_reads;
	size_t genome_len;

	size_t bootstrap_genome_len = 0;
	size_t bootstrap_num_reads = num_reads;
	size_t n_edges = num_edges();
	for (size_t i = 0; i < n_edges; i++)
		if (_edges[i].get_v1_idx() & 1)
			bootstrap_genome_len += _edges[i].length();

	float global_arrival_rate = FLOAT_DIV_NONZERO(bootstrap_num_reads,
						      bootstrap_genome_len);

	info("Bootstrapping with bootstream_genome_len = %zu, "
	     "bootstrap_num_reads = %zu, global_arrival_rate = %f",
	     bootstrap_genome_len, bootstrap_num_reads, global_arrival_rate);

	const size_t NUM_BOOTSTRAP_ITERATIONS = 3;
	const float SINGLE_COPY_THRESHOLD = 17.0;

	for (size_t i = 0; i < NUM_BOOTSTRAP_ITERATIONS; i++) {
		size_t num_unique_edges = 0;
		size_t num_optional_edges = 0;
		size_t num_required_edges = 0;

		bootstrap_genome_len = 0;
		bootstrap_num_reads = 0;
		for (size_t j = 0; j < n_edges; j++) {
			DirectedStringGraphEdge & e = _edges[j];
			size_t edge_len = _edges[j].length();
			unsigned edge_reads = _edges[j].get_mapped_read_count();
			float A_statistic = (global_arrival_rate * edge_len) -
					     (edge_reads * M_LN2);
			e.set_A_statistic(A_statistic);
			//info("edge_len = %zu, edge_reads = %u, A_statistic = %f",
			     //edge_len, edge_reads, A_statistic);
			if (A_statistic >= SINGLE_COPY_THRESHOLD) {
				bootstrap_genome_len += edge_len;
				bootstrap_num_reads += edge_reads;
				num_unique_edges++;
			} else if (edge_reads == 0) {
				num_optional_edges++;
			} else {
				num_required_edges++;
			}
		}
		global_arrival_rate = FLOAT_DIV_NONZERO(bootstrap_num_reads,
							bootstrap_genome_len);
		genome_len = DIV_NONZERO(num_reads, global_arrival_rate);
		info("Iteration %zu of %zu:  Estimated genome length "
		     "%zu", i, NUM_BOOTSTRAP_ITERATIONS, genome_len);
		info("num_unique_edges = %zu", num_unique_edges);
		info("num_optional_edges = %zu", num_optional_edges);
		info("num_required_edges = %zu", num_required_edges);
	}
}

void DirectedStringGraph::min_cost_circulation()
{
	static const int INFINITE_FLOW = 1000000;
	static const int HUGE_COST = 1000000;
	static const float SINGLE_COPY_THRESHOLD = 17.0;

	info("Adding special vertex and edges");
	v_idx_t n_verts = num_vertices();
	_vertices.resize(n_verts + 2);
	_vertices[n_verts].set_special();
	_vertices[n_verts + 1].set_special();

	for (v_idx_t v_idx = 0; v_idx < n_verts; v_idx++) {
		DirectedStringGraphEdge * e;

		//v_idx_t special_v_idx = n_verts | (v_idx & 1);
		for (v_idx_t special_v_idx = n_verts;
		     special_v_idx < n_verts + 2;
		     special_v_idx++)
		{
			e = &add_unlabeled_edge(v_idx, special_v_idx);
			e->set_special();

			e = &add_unlabeled_edge(special_v_idx, v_idx);
			e->set_special();
		}
	}

	n_verts += 2;
	edge_idx_t n_edges = num_edges();

	info("Creating lemon::SmartDigraph with %zu nodes and %zu arcs "
	     "and initializing network flow parameters",
	     n_verts, n_edges);
	lemon::SmartDigraph G;
	G.reserveNode(n_verts);
	G.reserveArc(n_edges);
	lemon::SmartDigraph::ArcMap<int> lower_map(G);
	lemon::SmartDigraph::ArcMap<int> upper_map(G);
	lemon::SmartDigraph::ArcMap<int> cost_map(G);
	lemon::SmartDigraph::NodeMap<int> supply_map(G);
	for (v_idx_t i = 0; i < n_verts; i++) {
		supply_map[G.addNode()] = 0;
	}
	foreach (DirectedStringGraphEdge & e, _edges) {
		lemon::SmartDigraph::Node node1 = G.nodeFromId(e.get_v1_idx());
		lemon::SmartDigraph::Node node2 = G.nodeFromId(e.get_v2_idx());
		lemon::SmartDigraph::Arc arc = G.addArc(node1, node2);
		int flow_lower_bound;
		int flow_upper_bound;
		int cost_per_unit_flow;
		if (e.is_special()) {
			flow_lower_bound = 0;
			flow_upper_bound = INFINITE_FLOW;
			cost_per_unit_flow = HUGE_COST;
		} else {
			if (e.get_A_statistic() >= SINGLE_COPY_THRESHOLD) {
				flow_lower_bound = 1;
				flow_upper_bound = 1;
			} else {
				if (e.get_num_inner_vertices() > 0)
					flow_lower_bound = 1;
				else
					flow_lower_bound = 0;
				flow_upper_bound = INFINITE_FLOW;
			}
			cost_per_unit_flow = 1;
		}
		lower_map[arc] = flow_lower_bound;
		upper_map[arc] = flow_upper_bound;
		cost_map[arc] = cost_per_unit_flow;
	}
	info("Initializing lemon::NetworkSimplex<lemon::SmartDigraph>");

	typedef lemon::NetworkSimplex<lemon::SmartDigraph> simplex_t;
	simplex_t simplex(G);
	simplex.lowerMap(lower_map);
	simplex.upperMap(upper_map);
	simplex.supplyMap(supply_map);
	simplex.costMap(cost_map);
	info("Running network simplex algorithm");
	simplex_t::ProblemType res = simplex.run();
	if (res == simplex_t::INFEASIBLE) {
		fatal_error("No feasible solution to min-cost circulation");
	} else if (res == simplex_t::UNBOUNDED) {
		fatal_error("Objective function unbounded");
	}
	info("Extracting network flow solution (total cost: %zu)",
	     simplex.totalCost());
	for (edge_idx_t edge_idx = 0; edge_idx < n_edges; edge_idx++) {
		lemon::SmartDigraph::Arc arc = G.arcFromId(edge_idx);
		_edges[edge_idx].set_traversal_count(simplex.flow(arc));
	}
	info("Done");
}

ext: cc


text:
// g++ cdf.cpp  -lboost_unit_test_framework

#include <cmath>
#include <boost/random.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/accumulators/numeric/functional/vector.hpp>
#include <boost/accumulators/numeric/functional/complex.hpp>
#include <boost/accumulators/numeric/functional/valarray.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/p_square_cumul_dist.hpp>

#include <vector>
#include <utility>

using namespace boost;
using namespace boost::accumulators;
//using namespace unit_test;

typedef iterator_range<std::vector<std::pair<double, double> >::iterator > histogram_type;

std::vector<std::pair<double, int> mapdata(histogram_type& h)
{
  const double mapToBlackPct = 0.2;
  const int otherSteps = 15;
  const double stepSizePct = (1-mapToBlackPct) / otherSteps;
  std::cout << "stepSizePct = " << stepSizePct << std::endl;

  std::vector<std::pair<double, int> mapping;

  for (std::size_t i = 0; i < h.size(); ++i)
    {
      int step = -1;
      if(h[i].second < mapToBlackPct)
	step = 0;
      else
	step = std::min(otherSteps, ceil((h[i].second - mapToBlackPct) / stepSizePct));
      std::cout << "[" << h[i].first << "," << h[i].second << "]"  << " ==> " << step << std::endl;

      // NOTE: Using -1 here so that the searches are easier later
      mapping.push_back(std::make_pair(h[i].first, step-1));
    }

  return mapping;
}

int main(void)
{

  // tolerance in %
  double epsilon = 3;

  typedef accumulator_set<double, stats<tag::p_square_cumulative_distribution> > accumulator_t;

  accumulator_t acc(tag::p_square_cumulative_distribution::num_cells = 100);

  // two random number generators
  boost::lagged_fibonacci607 rng;
  boost::normal_distribution<> mean_sigma(0,1);
  boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal(rng, mean_sigma);

  for (std::size_t i=0; i<100000; ++i)
    {
      acc(normal());
    }

  histogram_type histogram = p_square_cumulative_distribution(acc);
  //histogram_type histogram;

  for (std::size_t i = 0; i < histogram.size(); ++i)
    {
      std::cout << "[" << histogram[i].first << "," << histogram[i].second << "]" << std::endl;

      /*
      // problem with small results: epsilon is relative (in percent), not absolute!
      if ( histogram[i].second > 0.001 )
        BOOST_CHECK_CLOSE( 0.5 * (1.0 + erf( histogram[i].first / sqrt(2.0) )), histogram[i].second, epsilon );
      */
    }

  std::cout << "*************************" << std::endl;
  std::vector<std::pair<double, int> mapping (mapdata(histogram));

  std::vector<double> vals;
  vals.push_back(-0.75);
  vals.push_back(-0.25);
  vals.push_back(0);
  vals.push_back(0.5);
  vals.push_back(2);
  //vals.push_back(4.2);

  for(std::vector<double>::const_iterator iter; iter!=vals.end(); ++iter)
    {

    }

  return 0;
}

max_stars_repo_name: klaricmn/snippets

max_issues_repo_name: klaricmn/snippets

max_forks_repo_name: klaricmn/snippets


text:
#include "stddef.h"
#include <cstddef>
#include "utils.h"
#include "ps/ps.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <iterator>
#include <vector>
#include <unordered_set>
#include <thread>
#include <numeric>
#include <boost/program_options.hpp>
#include <limits>
#include <sstream>
#include <string>
#include <iostream>
#include <unistd.h>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <ctime>
#include <random>
#include <dirent.h>

#define MAX_STRING 100
#define EXP_TABLE_SIZE 1000
#define MAX_EXP 6
#define MAX_SENTENCE_LENGTH 1000
#define MAX_CODE_LENGTH 40

const int vocab_hash_size = 30000000, table_size = 1e8;  // Maximum 30 * 0.7 = 21M words in the vocabulary

// Precision switchable here
typedef float ValT;

struct vocab_word {
  long long cn;
  int *point;
  char *word, *code, codelen;
};

string in_file, out_file, vocab_save, vocab_retrieve;

struct vocab_word *vocab;
int binary = 0, debug_mode = 2, window = 5, min_reduce = 1, min_count,hs = 0, negative = 25;
int *vocab_hash,*table;
long long vocab_max_size = 1000, vocab_size = 0;
long long train_words = 0, word_count_actual = 0, file_size = 0;
ValT alpha, starting_alpha, sample;
long long int embed_dim,data_words;

bool shuffle_b;
bool write_results;
bool sync_push;
bool clustered_input;
bool localize_positives;
int prep_context_ahead;
bool downsample_replicated_negative_samples;
ValT *expTable;
double neg_power;
string init_model;
unsigned long long model_seed;
long max_runtime;

//shuffled vector acts as a scramble mapping of words - > keys
vector<ps::Key> forwards; // word -> key
vector<long long> backwards; // key -> word

// number of parameters to replicate
size_t replicate_n;

size_t peek_ahead_n;

using namespace ps;
using namespace std;

// get the key for the syn0 parameter of a word
inline Key syn0KeyResolver(long long word) {
  if (shuffle_b) {
    return (forwards[word] * 2);
  } else {
    return (word * 2);
  }
}
// get the key for the syn1 parameter of a word
inline Key syn1KeyResolver(long long word) {
  if (shuffle_b) {
    return (forwards[word] * 2 + 1);
  } else {
    return (word * 2 + 1);
  }
}
// get the word from a syn1 key
inline long long syn1Reverse(Key key) {
  if (shuffle_b) {
    return backwards[(key - 1) / 2];
  } else {
    return (key - 1)/2;
  }
}

typedef DefaultColoServerHandle <ValT> HandleT;
typedef ColoKVServer <ValT, HandleT> ServerT;
typedef ColoKVWorker <ValT, HandleT> WorkerT;

// Config
uint num_workers = -1;
size_t num_iterations = 0;
size_t num_threads = 0;
Key num_keys = 0;

// calculates weight for negative sampling distribution
inline double unigram_pow(int a, unordered_set<Key>& replicated, int num_servers) {
  double pw = pow(vocab[a].cn, neg_power);

  // adjust weight of the negative samples that are replicated
  // (because such negative samples will locally available 100% of the time)
  if (downsample_replicated_negative_samples) {
    if (replicated.find(syn1KeyResolver(a)) != replicated.end()) {
      pw = pw / num_servers;
    }
  }
  return pw;
}

// construct a table of all negative samples for efficient sampling
void InitUnigramTable(const vector<Key>& replicated_parameters) {
  // replicated parameters as set
  unordered_set<Key> replicated (replicated_parameters.begin(), replicated_parameters.end());
  int num_servers = atoi(Environment::Get()->find("DMLC_NUM_SERVER"));

  int a, i;
  double train_words_pow = 0;
  double d1;
  table = (int *) malloc(table_size * sizeof(int));
  for (a = 0; a < vocab_size; a++) train_words_pow += unigram_pow(a, replicated, num_servers);
  i = 0;
  d1 = unigram_pow(i, replicated, num_servers) / train_words_pow;
  for (a = 0; a < table_size; a++) {
    table[a] = i;
    if (a / (double) table_size > d1) {
      if (i == 1 || i == 10 || i == 100 || i == 1000 || i == 10000) {
        ADLOG(setw(5) << setfill(' ') << i << " most frequent words have " << d1 << " sampling probability");
      }
      i++;
      d1 += unigram_pow(i, replicated, num_servers) / train_words_pow;
    }
    if (i >= vocab_size) i = vocab_size - 1;
  }
}

// Reads a single word from a file, assuming space + tab + EOL to be word boundaries
void ReadWord(char *word, FILE *fin, char *eof) {
  int a = 0, ch;
  while (1) {
    ch = fgetc_unlocked(fin);
    if (ch == EOF) {
      *eof = 1;
      break;
    }
    if (ch == 13) continue; // carriage return
    if ((ch == ' ') || (ch == '\t') || (ch == '\n')) {
      if (a > 0) {
        if (ch == '\n') ungetc(ch, fin);
        break;
      }
      if (ch == '\n') {
        std::strcpy(word, (char *) "</s>");
        return;
      } else continue;
    }
    word[a] = ch;
    a++;
    if (a >= MAX_STRING - 1) a--;   // Truncate too long words
  }
  word[a] = 0;
}

// Returns hash value of a word
int GetWordHash(char *word) {
  unsigned long long a, hash = 0;
  for (a = 0; a < strlen(word); a++) hash = hash * 257 + word[a];
  hash = hash % vocab_hash_size;
  return hash;
}

// Returns position of a word in the vocabulary; if the word is not found, returns -1
int SearchVocab(char *word) {
  unsigned int hash = GetWordHash(word);
  while (1) {
    if (vocab_hash[hash] == -1) return -1;
    if (!strcmp(word, vocab[vocab_hash[hash]].word)) return vocab_hash[hash];
    hash = (hash + 1) % vocab_hash_size;
  }
  return -1;
}

// Reads a word and returns its index in the vocabulary
int ReadWordIndex(FILE *fin, char *eof) {
  char word[MAX_STRING], eof_l = 0;
  ReadWord(word, fin, &eof_l);

  if (eof_l) {
    *eof = 1;
    return -1;
  }

  return SearchVocab(word);
}

// Adds a word to the vocabulary
int AddWordToVocab(char *word) {
  unsigned int hash, length = strlen(word) + 1;
  if (length > MAX_STRING) length = MAX_STRING;
  vocab[vocab_size].word = (char *) calloc(length, sizeof(char));
  std::strcpy(vocab[vocab_size].word, word);
  vocab[vocab_size].cn = 0;
  vocab_size++;
  // Reallocate memory if needed
  if (vocab_size + 2 >= vocab_max_size) {
    vocab_max_size += 1000;
    vocab = (struct vocab_word *) realloc(vocab, vocab_max_size * sizeof(struct vocab_word));
  }
  hash = GetWordHash(word);
  while (vocab_hash[hash] != -1) hash = (hash + 1) % vocab_hash_size;
  vocab_hash[hash] = vocab_size - 1;
  return vocab_size - 1;
}

// Used later for sorting by word counts
int VocabCompare(const void *a, const void *b) {
  long long l = ((struct vocab_word *) b)->cn - ((struct vocab_word *) a)->cn;
  if (l > 0) return 1;
  if (l < 0) return -1;
  return 0;
}

// Sorts the vocabulary by frequency using word counts
void SortVocab() {
  int a, size;
  unsigned int hash;
  // Sort the vocabulary and keep </s> at the first position
  qsort(&vocab[1], vocab_size - 1, sizeof(struct vocab_word), VocabCompare);
  for (a = 0; a < vocab_hash_size; a++) vocab_hash[a] = -1;
  size = vocab_size;
  train_words = 0;
  for (a = 0; a < size; a++) {
    // Words occuring less than min_count times will be discarded from the vocab
    if ((vocab[a].cn < min_count) && (a != 0)) {
      vocab_size--;
      free(vocab[a].word);
    } else {
      // Hash will be re-computed, as after the sorting it is not actual
      hash = GetWordHash(vocab[a].word);
      while (vocab_hash[hash] != -1) hash = (hash + 1) % vocab_hash_size;
      vocab_hash[hash] = a;
      train_words += vocab[a].cn;
    }
  }
  vocab = (struct vocab_word *) realloc(vocab, (vocab_size + 1) * sizeof(struct vocab_word));

}

// Reduces the vocabulary by removing infrequent tokens
void ReduceVocab() {
  int a, b = 0;
  unsigned int hash;
  for (a = 0; a < vocab_size; a++)
    if (vocab[a].cn > min_reduce) {
      vocab[b].cn = vocab[a].cn;
      vocab[b].word = vocab[a].word;
      b++;
    } else free(vocab[a].word);
  vocab_size = b;
  for (a = 0; a < vocab_hash_size; a++) vocab_hash[a] = -1;
  for (a = 0; a < vocab_size; a++) {
    // Hash will be re-computed, as it is not actual
    hash = GetWordHash(vocab[a].word);
    while (vocab_hash[hash] != -1) hash = (hash + 1) % vocab_hash_size;
    vocab_hash[hash] = a;
  }
  fflush(stdout);
  min_reduce++;
}

void LearnVocabFromTrainFile() {

  char word[MAX_STRING], eof = 0;
  FILE *fin;
  long long a, i, wc = 0;
  for (a = 0; a < vocab_hash_size; a++) vocab_hash[a] = -1;

  fin = fopen(in_file.c_str(), "rb");

  if (fin == NULL) {
    printf("ERROR: training data file not found!\n");
    exit(1);
  }
  vocab_size = 0;
  AddWordToVocab((char *) "</s>");
  while (1) {
    ReadWord(word, fin, &eof);
    if (eof) break;
    train_words++;
    wc++;
    if ((debug_mode > 1) && (wc >= 1000000)) {
      printf("%lldM%c", train_words / 1000000, 13);
      fflush(stdout);
      wc = 0;
    }
    i = SearchVocab(word);
    if (i == -1) {
      a = AddWordToVocab(word);
      vocab[a].cn = 1;
    } else vocab[i].cn++;
    if (vocab_size > vocab_hash_size * 0.7) ReduceVocab();
  }
  SortVocab();
  if (debug_mode > 0) {
    printf("Vocab size: %lld\n", vocab_size);
    printf("Words in train file: %lld\n", train_words);
  }
  file_size = ftell(fin);
  fclose(fin);
}

void SaveVocab() {
  long long i;
  FILE *fo = fopen(vocab_save.c_str(), "wb");
  if (fo == NULL) {
    perror("FAILED to open file to save vocab: ");
}
  for (i = 0; i < vocab_size; i++) fprintf(fo, "%s %lld\n", vocab[i].word, vocab[i].cn);
  fclose(fo);
}

void ReadVocab() {
  long long a, i = 0;
  char c, eof = 0;
  int irrelevant_returnval = 0; // only here so the compiler does not warn anymore, as fscanf has a unused returnvalue
  irrelevant_returnval++;
  irrelevant_returnval--;
  char word[MAX_STRING];
  FILE *fin = fopen(vocab_retrieve.c_str(), "rb");
  if (fin == NULL) {
    printf("Vocabulary file not found\n");
    exit(1);
  }
  for (a = 0; a < vocab_hash_size; a++) vocab_hash[a] = -1;
  vocab_size = 0;
  while (1) {
    ReadWord(word, fin, &eof);
    if (eof) break;
    a = AddWordToVocab(word);
    irrelevant_returnval = fscanf(fin, "%lld%c", &vocab[a].cn, &c);
    i++;
  }
  SortVocab();
  if (debug_mode > 0) {
    printf("Vocab size: %lld\n", vocab_size);
    printf("Words in train file: %lld\n", train_words);
  }
  fin = fopen(in_file.c_str(), "rb");
  if (fin == NULL) {
    printf("ERROR: training data file not found!\n");
    exit(1);
  }
  fseek(fin, 0, SEEK_END);
  file_size = ftell(fin);
  fclose(fin);
}

// writes current embeddings either in binary or human readable format
void write_checkpoint(string output, WorkerT &kv, const bool write_syn1=false) {
  util::Stopwatch sw; sw.start();

  // pull entire model
  vector<Key> syn0_keys (vocab_size);
  vector<Key> syn1_keys (vocab_size);
  for (long a = 0; a < vocab_size; a++) {
    syn0_keys[a] = syn0KeyResolver(a);
    syn1_keys[a] = syn1KeyResolver(a);
  }

  vector<ValT> syn0 (embed_dim * vocab_size);
  vector<ValT> syn1 (embed_dim * vocab_size);

  kv.WaitReplicaSync();
  kv.Wait(kv.Pull(syn0_keys, &syn0));
  if (write_syn1) kv.Wait(kv.Pull(syn1_keys, &syn1));

  if (binary == 1) { // write binary files
    FILE *fo_syn0;
    FILE *fo_syn1;

    fo_syn0 = fopen(output.c_str(), "wb");
    if (write_syn1) fo_syn1 = fopen((output+".syn1").c_str(), "wb");
    fprintf(fo_syn0, "%lld %lld\n", vocab_size, embed_dim);
    for (long a = 0; a < vocab_size; a++) {
      fprintf(fo_syn0, "%s ", vocab[a].word);
      fwrite(&syn0[a*embed_dim], sizeof(ValT), embed_dim, fo_syn0);
      if (write_syn1) fwrite(&syn1[a*embed_dim], sizeof(ValT), embed_dim, fo_syn1);
      fprintf(fo_syn0, "\n");
    }
    fclose(fo_syn0);
    if (write_syn1) fclose(fo_syn1);

  } else { // write a text file
    ofstream file;
    file.open(output.c_str());

    file << vocab_size << " " << embed_dim << endl;
    for (long a = 0; a < vocab_size; a++) {
      file << vocab[a].word;
      for (long b = 0; b < embed_dim; b++) {
        file << syn0[a*embed_dim+b] << endl;
      }
    }
    file.close();
  }
  sw.stop();
  ADLOG("Current embeddings have been written to '" << output << "' (" << sw << ")");
}

unsigned long long negs_next_random;
// draw one negative sample (i.e., returns one syn1 key, sampled according to a specific distribution)
inline Key DrawNegSample() {
  negs_next_random = negs_next_random * (unsigned long long) 25214903917 + 11;
  unsigned long long target = table[(negs_next_random >> 16) % table_size];
  if (target == 0) {
    target = negs_next_random % (vocab_size - 1) + 1;
  }
  return syn1KeyResolver(target);
}

// A "peekable" random number generator. I.e., this generator allows for
// finding out random numbers that will be generated in the future (by "peeking")
template<typename T>
struct PeekableRandom {
  PeekableRandom(T seed): next_random{seed} {}

  // get the next random number
  T next() {
    // sample new items if necessary
    if (q.size() < 1) {
      sample(1);
    }
    // pop first element
    auto e = q[0];
    q.pop_front();
    return e;
  }

  // peek `num_steps_ahead` into the future
  T peek(size_t num_steps_ahead) {
    if (q.size() <= num_steps_ahead) {
      sample(num_steps_ahead);
    }
    return q[num_steps_ahead];
  }

private:
  // generate a batch of random values
  void sample(size_t min) {
    size_t to_sample = max(min, sample_batch);
    for(size_t i=0; i!=to_sample; ++i) {
      next_random = next_random * (T) 25214903917 + 11;
      q.push_back(next_random);
    }
  }

  // debug output
  std::string print() {
    std::stringstream s;
    for(size_t i=0; i!=q.size(); ++i) {
      s << q[i] << " ";
    }
    return s.str();
  }

  std::deque<T> q; // future random numbers
  T next_random; // random state
  const size_t sample_batch = 100; // minimum size of random numbers to generate
};

//training/computation happens here, equivalent to the original training-thread.
void training_thread(WorkerT &kv, int customer_id, int worker_id) {
  long long a, b, d, word, last_word, sentence_length = 0, sentence_position = 0;
  long long word_count = 0, last_word_count = 0, sen[MAX_SENTENCE_LENGTH + 1];
  long long c, target, label, local_iter = num_iterations;
  PeekableRandom<unsigned long long> context_size_rand((long long) 7+worker_id*13); // rng for determining context size
  unsigned long long next_random_read = (long long) 11+worker_id*17; // rng for frequent word skips
  char eof = 0;

  ValT f, g;
  ADLOG("[w" << worker_id << "] begins work" << endl);

  // allocate memory for local parameters
  vector <ValT> syn0_vec     (embed_dim);
  vector <ValT> syn1neg_vec  (embed_dim);
  vector <ValT> neu1_vec     (embed_dim);
  vector <ValT> neu1e_vec    (embed_dim);
  vector <ValT> syn1neg_push (embed_dim);
  vector <Key>  syn0_key     (1);
  vector <Key>  syn1neg_key  (1);

  // groups of negative samples (one per context)
  std::unordered_map<long long, SampleID> sample_ids {};
  long long num_context = 0;
  long long future_context = 0;

  FILE *fi = fopen(in_file.c_str(), "rb");

  //partitions file for threads
  if (clustered_input) {
    fseek(fi, (file_size / (long long) num_threads) * (long long) (customer_id - 1), SEEK_SET);

  } else {
    auto pos = (file_size / (long long) num_workers) * (long long) worker_id;
    ADLOG("Worker " << worker_id << ": Start at position " << pos << " of " << file_size);
    fseek(fi, pos, SEEK_SET);
  }

  // wait until all workers are ready to start
  kv.Barrier();

  util::Stopwatch sw_epoch;
  util::Stopwatch sw_epoch_all;
  util::Stopwatch sw_train;
  util::Stopwatch sw_worker;
  sw_train.start();
  sw_worker.start();
  sw_epoch.start();
  sw_epoch_all.start();
  util:: Stopwatch sw_wait_syn0;
  util:: Stopwatch sw_wait_syn1_pos;
  util:: Stopwatch sw_find_neg;
  //train loop
  while (1) { //loop ends when thread has reached the end of its partition during its last iteration.

    //adjusts learning rate (alpha)
    if (word_count - last_word_count > 10000) {
      word_count_actual += word_count -last_word_count;
      last_word_count = word_count;

      alpha = starting_alpha *
        (1 - (word_count_actual * ps::NumServers()) / (ValT) (num_iterations * train_words + 1));

      if (alpha < starting_alpha * 0.0001) alpha = starting_alpha * 0.0001;
    }

    // read next sentence
    if (sentence_length == 0) {
      std::unordered_set<Key> to_localize {};

      while (1) {

        word = ReadWordIndex(fi, &eof);  // gets the position in the vocab; common words come first, uncommon words last

        if (eof) break;
        if (word == -1) continue; // word not in vocab
        word_count++;
        if (word == 0) break;
        // The subsampling randomly discards frequent words while keeping the ranking same
        if (sample > 0) {
          ValT ran = (sqrt(vocab[word].cn / (sample * train_words)) + 1) * (sample * train_words) / vocab[word].cn;
          next_random_read= next_random_read * (unsigned long long) 25214903917 + 11;
          if (ran < (next_random_read & 0xFFFF) / (ValT) 65536) continue;
        }

        // localize parameters for the first couple of words of the sentence
        // (later words will be localized by the context prep)
        if (sentence_length < prep_context_ahead+window) {
          to_localize.insert(syn0KeyResolver(word));
          if (sentence_length < prep_context_ahead) {
            to_localize.insert(syn1KeyResolver(word));
          }
        }

        sen[sentence_length] = word;
        sentence_length++;
        if (sentence_length >= MAX_SENTENCE_LENGTH) break;
      }
      sentence_position = 0;

      // "peek ahead": localize a few words from the beginning of the next sentence
      if (peek_ahead_n != 0 && localize_positives) {
        auto read_pos = ftell(fi); // store current read position (so we can jump back below)
        auto next_random_read_copy = next_random_read; // fork the generator state so we can look into the future
        for (size_t i=0; i!=peek_ahead_n; ++i) {
          long long next_word = ReadWordIndex(fi, &eof);
          // skip conditions (copied from above to achieve identical behavior)
          if (eof) break;
          if (next_word == -1) { --i; continue; }; // word not in vocab
          if (next_word == 0)  { break; }; // sentence end
          if (sample > 0) { // discard frequent words from time to time
            ValT ran = (sqrt(vocab[next_word].cn / (sample * train_words)) + 1) * (sample * train_words) / vocab[next_word].cn;
            next_random_read_copy = next_random_read_copy * (unsigned long long) 25214903917 + 11;
            if (ran < (next_random_read_copy & 0xFFFF) / (ValT) 65536) {
              --i;
              continue;
            }
          }
          to_localize.insert(syn0KeyResolver(next_word));
          to_localize.insert(syn1KeyResolver(next_word));
        }
        fseek(fi, read_pos, SEEK_SET); // reset read position
      }

      // localize the words that we just read
      if (localize_positives) {
        std::vector<Key> to_localize_vec(to_localize.begin(), to_localize.end());
        kv.Localize(to_localize_vec);
      }
    }

    //finish iteration/epoch when either eof is reached or the # word_count
    if (eof || (word_count > min(train_words,data_words) / num_workers)) {
      word_count_actual += word_count - last_word_count;
      local_iter--;
      auto fepoch = num_iterations - local_iter; //finished epochs

      sw_epoch.stop();
      ADLOG("Worker " << worker_id << " finished epoch " << fepoch << " (" << sw_epoch << ")." );
      kv.Barrier();
      sw_epoch_all.stop();
      sw_train.stop();
      if (worker_id == 0) {
        ADLOG("All workers finished epoch " << fepoch << " (epoch: " << sw_epoch_all << ", total: " << sw_train << ")");
      }

      if (write_results) {
        kv.WaitAll();
        kv.Barrier();
        if (customer_id == 1 && ps::MyRank() == ps::NumServers()-1) {// last rank saves (usually this is the first node)
          util::Stopwatch sw_write; sw_write.start();
          ADLOG("Write epoch " << fepoch << " embeddings");
          write_checkpoint(out_file + ".epoch." + to_string(fepoch), kv);
        }
      }
      kv.Barrier();

      if (local_iter == 0) {
        ADLOG("[w" << worker_id << "] Wait syn0: " << sw_wait_syn0 << "  Wait syn1 positive: " << sw_wait_syn1_pos << "  Find negative: " << sw_find_neg);
        break;
      }

      // maximum time
      if (sw_train.elapsed_s() > max_runtime ||
          sw_train.elapsed_s() + sw_epoch_all.elapsed_s() > max_runtime * 1.05) {
        ADLOG("Worker " << worker_id << " stops after epoch " << fepoch << " because max. time is reached: " << sw_train.elapsed_s() << "s (+1 epoch) > " << max_runtime << "s (epoch: " << sw_epoch_all.elapsed_s() << "s)");
        break;
      }

      sw_epoch.start();
      sw_epoch_all.start();
      sw_train.resume();
      //variable reset for next epoch, worker starts at its determined starting-point again.
      word_count = 0;
      last_word_count = 0;
      sentence_length = 0;
      eof = 0;

      if (clustered_input) {
        fseek(fi, (file_size / (long long) num_threads) * (long long) (customer_id - 1), SEEK_SET);
        continue;
      } else {
        fseek(fi, (file_size / (long long) num_workers) * (long long) worker_id, SEEK_SET); // SEEK_SET begins at file start; computes offset via middle value
        continue;
      }
    }

    word = sen[sentence_position]; //word is extracted here again
    if (word == -1) continue;

    std::fill(neu1_vec.begin(), neu1_vec.end(), 0);
    std::fill(neu1e_vec.begin(), neu1e_vec.end(), 0);

    // prepare a future context
    while (future_context <= num_context + prep_context_ahead) {
      // find out the window size for the future context
      int b_future = context_size_rand.peek(future_context-num_context) % window;

      // prep a negative sample of adequate size for this future context
      sample_ids[future_context] = kv.PrepareSample((window-b_future) * 2 * negative);

      // localize words for this future context
      auto future_position = sentence_position + (future_context-num_context);
      if (localize_positives && future_position < sentence_length) {
        std::vector<Key> to_localize {};
        to_localize.reserve((window-b_future) * 2 + 1);
        to_localize.push_back(syn1KeyResolver(sen[future_position]));
        for (a = b_future; a < window * 2 + 1 - b_future; a++) {
          c = future_position - window + a;
          if (a != window && c >= 0 && c < sentence_length) {
            to_localize.push_back(syn0KeyResolver(sen[c]));
          }
        }
        kv.Localize(to_localize);
      }
      ++future_context;
    }

    // determine window size for the current context
    b = context_size_rand.next() % window;

    // Iterate through one context (number of word pairs = (window-b)*2)
    for (a = b; a < window * 2 + 1 - b; a++) {
      if (a != window) {
        c = sentence_position - window + a;

        if (c < 0) continue;
        if (c >= sentence_length) continue;
        last_word = sen[c];
        if (last_word == -1) continue;

        syn0_key[0] = syn0KeyResolver(last_word);
        sw_wait_syn0.resume();
        kv.Wait(kv.Pull(syn0_key, &syn0_vec));
        sw_wait_syn0.stop();

        std::fill(neu1e_vec.begin(), neu1e_vec.end(), 0);

        //negative sampling
        for (d = 0; d < negative + 1; d++) {

          if (d == 0) { // positive (data point)
            label = 1;
            target = word;
            syn1neg_key[0] = syn1KeyResolver(target); // precomputed in preload neg_samples

            sw_wait_syn1_pos.resume();
            kv.Localize(syn1neg_key);
            kv.Wait(kv.Pull(syn1neg_key, &syn1neg_vec));
            sw_wait_syn1_pos.stop();
          } else { // negative sample
            label = 0;

            // Retrieve a negative sample
            sw_find_neg.resume();
            kv.Wait(kv.PullSample(sample_ids[num_context], syn1neg_key, syn1neg_vec));
            target = syn1Reverse(syn1neg_key[0]);
            sw_find_neg.stop();

            if (target == word)
              continue;
          }

          // retrieve output layer of negative-sampled word
          f = 0;
          for (c = 0; c < embed_dim; c++)f += syn0_vec[c] * syn1neg_vec[c];

          if (f > MAX_EXP) g = (label - 1) * alpha;
          else if (f < -MAX_EXP) g = (label - 0) * alpha;
          else g = (label - expTable[(int) ((f + MAX_EXP) * (EXP_TABLE_SIZE / MAX_EXP / 2))]) * alpha;

          for (c = 0; c < embed_dim; c++) neu1e_vec[c] += g * syn1neg_vec[c];

          for (c = 0; c < embed_dim; c++) syn1neg_push[c] = g * syn0_vec[c];

          // publish/push updates;  hidden->output
          if (sync_push) {
            kv.Wait(kv.Push(syn1neg_key, syn1neg_push));
          } else {
            kv.Push(syn1neg_key, syn1neg_push);
          }

        }
        // Learn weights input -> hidden
        if (sync_push) {
          kv.Wait(kv.Push(syn0_key, neu1e_vec));
        } else {
          kv.Push(syn0_key, neu1e_vec);
        }
      }
    }
    kv.FinishSample(sample_ids[num_context]);
    sample_ids.erase(num_context);
    sentence_position++;
    num_context++;

    if (sentence_position >= sentence_length) {
      sentence_length = 0;
      continue;
    }

  }

  sw_worker.stop();
  ADLOG("[w" << worker_id << "] finished training (" << sw_worker << "). Processed " << word_count << " words.");

  kv.Barrier();
  sw_train.stop();
  if (worker_id == 0) {
    ADLOG("All workers finished training (" << sw_train << ")");
  }

  kv.WaitAll(); // wait until all requests of this worker are answered
  kv.Barrier(); // make sure all requests in the system have been answered

  fclose(fi);
  ADLOG("[w" << worker_id << "] has passed the training barrier ");
}

// handles loading vocab into memory, found in the original TrainModel()
void load_vocab(int worker_id) {
  if (vocab_retrieve.size() != 0) {
    ReadVocab();

  } else {
    LearnVocabFromTrainFile();

    ADLOG("[w" << worker_id << "] ______ done loading vocab from'" << in_file << "' into memory" << endl);

  }

  if (vocab_save.size() != 0) {
    SaveVocab();
    ADLOG("______Vocab was written to::'" << vocab_save << "'");
  }
}

int ArgPos(char *str, int argc, char **argv) {
  int a;
  for (a = 1; a < argc; a++)
    if (!strcmp(str, argv[a])) {
      if (a == argc - 1) {
        printf("Argument missing for %s\n", str);
        exit(1);
      }
      return a;
    }
  return -1;
}

// Initialize first layer (syn0) of network
void initialize_model(WorkerT &kv) {
  if (init_model.compare("none") == 0) { // no model init
    ALOG("No parameter initialization");
  } else if (init_model.compare("random") == 0) { // init syn0 randomly, don't init syn1
    ALOG("Initialize model randomly (seed " << model_seed << ")");
    vector <ValT> syn_vec(embed_dim * vocab_size);
    vector <Key> syn_key(vocab_size);
    unsigned long long next_random = model_seed;
    for (int e = 0; e < vocab_size; ++e) {
      syn_key[e] = syn0KeyResolver(e);
      for (int j = 0; j < embed_dim; ++j) {
        next_random = next_random * (unsigned long long) 25214903917 + 11; // taken from original w2v
        syn_vec[e * embed_dim + j] =
          (((next_random & 0xFFFF) / (ValT) 65536) - 0.5) / embed_dim;
      }
    }
    kv.Wait(kv.Push(syn_key, syn_vec));
    ALOG("Model init finished: " << syn_vec[0] << " " << syn_vec[1] << " " << syn_vec[2] << " .. ");
  } else { // load syn0 and syn1 from a checkpoint
    ALOG("Initialize model from checkpoint:" << init_model);
    // create keys
    vector<Key> syn0_keys (vocab_size);
    vector<Key> syn1_keys (vocab_size);
    for (long a = 0; a < vocab_size; a++) {
      syn0_keys[a] = syn0KeyResolver(a);
      syn1_keys[a] = syn1KeyResolver(a);
    }

    vector <ValT> syn0 (embed_dim * vocab_size);
    vector<ValT> syn1 (embed_dim * vocab_size);

    // load files
    FILE *fo_syn0;
    FILE *fo_syn1;
    long long read_vocab_size;
    long long read_embed_dim;
    // long read_word;
    char read_word[MAX_STRING];
    fo_syn0 = fopen(init_model.c_str(), "rb");
    fo_syn1 = fopen((init_model+".syn1").c_str(), "rb");
    if (fo_syn0 == NULL || fo_syn1 == NULL) {
      ALOG("Error opening checkpoint files (" << init_model << ")");
      abort();
    }
    fscanf(fo_syn0, "%lld %lld\n", &read_vocab_size, &read_embed_dim);
    assert(read_vocab_size == vocab_size);
    assert(read_embed_dim == embed_dim);

    // read word by word
    for (long a = 0; a < vocab_size; a++) {
      fscanf(fo_syn0, "%s", read_word);
      fgetc(fo_syn0); // read space
      assert(strcmp(read_word, vocab[a].word) == 0);

      fread(&syn0[a*embed_dim], sizeof(ValT), embed_dim, fo_syn0);
      fread(&syn1[a*embed_dim], sizeof(ValT), embed_dim, fo_syn1);
      fgetc(fo_syn0); // read line break
    }
    fclose(fo_syn0);
    fclose(fo_syn1);

    // push parameters into the server
    kv.Wait(kv.Push(syn0_keys, syn0), kv.Push(syn1_keys, syn1));
  }
  kv.WaitReplicaSync();
}

// returns a vector of files inside given directory
void read_directory(const std::string &name, vector <string> &v) {
  DIR *dirp = opendir(name.c_str());
  struct dirent *dp;
  while ((dp = readdir(dirp)) != NULL) {
    v.push_back(dp->d_name);
  }
  closedir(dirp);
}

// returns "file.0x.of.08.txt" -type files(given "file"-pattern) for clustered inputs.
string clustered_ingest(string target_file) {
  vector <string> dirlist;
  string dirname = "";
  string result = "";

  // strips filename to get path to directory
  const size_t last_slash_idx = target_file.find_last_of("\\/");
  if (std::string::npos != last_slash_idx) {
    dirname = target_file.substr(0, last_slash_idx);
    target_file.erase(0, last_slash_idx + 1);
  }
  read_directory(dirname, dirlist);

  string sub_fid;
  string target = "";

  //selects filename which substring matches the given filename pattern

  for (auto file_in_dir: dirlist) {
    sub_fid = file_in_dir.substr(0, target_file.length());
    if (target_file.compare(sub_fid) == 0) {

      target = file_in_dir;
      break;
    }

  }
  dirname = dirname + "/";
  result = dirname + target;
  CHECK(result.length() > target_file.length()) << " no suitable Candidate found, check in " << dirname
                                                << "-directory" << endl;
  return result;
}

// initializes data structures that are shared among worker threads
void init_shared_datastructures () {

  // adapt input file name for reading pre-clustered data (if desired)
  if (clustered_input) {
    string slash = "/";
    std::size_t found = in_file.find(slash);
    CHECK(found != std::string::npos) << " full path needed for clustered input" << endl;
    in_file = clustered_ingest(in_file);
    ADLOG("Clustered Input:: Inputfile will be loaded as:: '" << in_file << "' " << endl);
  }

  // load vocabulary and initialize exponent table
  vocab = (struct vocab_word *) calloc(vocab_max_size, sizeof(struct vocab_word));
  vocab_hash = (int *) calloc(vocab_hash_size, sizeof(int));
  expTable = (ValT *) malloc((EXP_TABLE_SIZE + 1) * sizeof(ValT));
  for (int i = 0; i < EXP_TABLE_SIZE; i++) {
    expTable[i] = exp((i / (ValT) EXP_TABLE_SIZE * 2 - 1) *
                      MAX_EXP); // Precompute the exp() table
    expTable[i] = expTable[i] / (expTable[i] + 1); // Precompute f(x) = x / (x + 1)
  }
  // load vocabulary
  load_vocab(0);

  // initialize "shuffling tables"
  // we use this to assign random keys to parameters. this results in more even
  // load on the different servers
  forwards.resize(num_keys / 2 - 2);
  iota(forwards.begin(), forwards.end(), 0);
  unsigned seed = 2;
  shuffle(forwards.begin(), forwards.end(), default_random_engine(seed));
  backwards.resize(forwards.size());
  for(size_t i=0; i!=forwards.size(); ++i) {
    backwards[forwards[i]] = i;
  }

  // sanity checks on maximum key (given as program option)
  unsigned int min_key = vocab_size * 2 + 5, max_key = vocab_size * 2 + 10;
  //  +10 offset kept for possible keys being used as debug flags
  CHECK(num_keys >= min_key && num_keys <= max_key) << "number of keys should be set between " << min_key << " and "
                                                    << max_key << " (vocab size " << vocab_size << ")";

  ADLOG("______Global datastructures have been set up.______" << endl);
}

// sort a collection of pairs by the first element, in descending order
bool sortdesc(const pair<double, Key>& a, const pair<double, Key>& b) {
  return (a.first > b.first);
}

// determine which parameters to replicate.
// we calculate the expected access frequency for each parameter
// and then pick the N with the highest expectations for replication
std::vector<Key> determine_hotspot_keys(size_t N_hotspots) {
  std::vector<Key> to_replicate {};
  if (N_hotspots == 0) return to_replicate;

  if (N_hotspots > vocab_size*2) N_hotspots = vocab_size*2; // cannot replicate more than all keys

  to_replicate.reserve(N_hotspots);

  // calculate the expected access frequency for each parameter
  std::vector<std::pair<double,Key>> expected_frequency (num_keys);
  // calculate total count and pow
  double total_pow = 0;
  long total_cn = 0;
  for (long long a = 0; a != vocab_size; ++a) {
    total_cn += vocab[a].cn;
    total_pow += pow(vocab[a].cn, neg_power);
  }
  for (long long a=0; a != vocab_size; ++a) {
    // for syn0, access frequency depends on data frequency of the corresponding word
    expected_frequency[syn0KeyResolver(a)] = std::make_pair(1.0*vocab[a].cn/total_cn, syn0KeyResolver(a));
    // for syn1, access frequency depends on (1) sampling distribution (majority of accesses) and
    // (2) data frequency
    expected_frequency[syn1KeyResolver(a)] =
      std::make_pair(pow(vocab[a].cn, neg_power)/total_pow*negative+
                     1.0*vocab[a].cn/total_cn, syn1KeyResolver(a));
    // note: over all, parameters of layer syn1 are accessed much more frequently than parameters in syn0.
    // this is reflected here
  }

  // replicate the `N_hotspot` parameters that we expect to be accessed most frequently
  std::sort(expected_frequency.begin(), expected_frequency.end(), sortdesc);
  for (size_t i=0; i != N_hotspots; ++i) {
    // ADLOG("Replicate " << expected_frequency[i].second << ", expected accesses: " << expected_frequency[i].first);
    to_replicate.push_back(expected_frequency[i].second);
  }
  assert(to_replicate.size() == N_hotspots);

  return to_replicate;
}

void RunWorker(int customer_id, ServerT *server = nullptr) {

  Start(customer_id);
  std::unordered_map <std::string, util::Stopwatch> sw{};
  WorkerT kv(0, customer_id, *server);

  int worker_id = ps::MyRank() * num_threads + customer_id - 1; // a unique id for this worker thread

  // initialize the model with random parameter values
  kv.BeginSetup();
  if (customer_id == 1 && ps::MyRank() == ps::NumServers()-1) { // last rank initializes
    initialize_model(kv);
  }
  kv.EndSetup();

  training_thread(kv, customer_id, worker_id);

  kv.Finalize();
  Finalize(customer_id, false);
}

//boost program options
int process_program_options(const int argc, const char *const argv[]) {

  namespace po = boost::program_options;
  po::options_description desc("Allowed options");

  desc.add_options()
    ("help,h", "produce help message")
    ("num_threads,t", po::value<size_t>(&num_threads)->default_value(2),
     "number of worker threads to run (per process)")
    ("num_iterations,i", po::value<size_t>(&num_iterations)->default_value(15), "number of iterations to run")
    ("input_file,f", po::value<string>(&in_file), "name of the training file")
    ("output_file,o", po::value<string>(&out_file)->default_value("vectors.bin"), "output file (to store word vectors)")
    ("vocab_save", po::value<string>(&vocab_save), "name of the resulting vocab-file")
    ("vocab_retrieve", po::value<string>(&vocab_retrieve), "name of the source vocab-file")
    ("shuffle", po::bool_switch(&shuffle_b)->default_value(true), "boolean to scramble words on keys randomly")
    ("debug_mode,d", po::value<int>(&debug_mode)->default_value(2), "disables debug mode")
    ("window,w", po::value<int>(&window)->default_value(5), "adjusts sizing of word-window, default is 5")
    ("embed_dim,v", po::value<long long int>(&embed_dim)->default_value(200),
     "number of values per key; so embed_dim")
    ("num_keys,k", po::value<Key>(&num_keys)->default_value(10), "number of parameters")
    ("negative", po::value<int>(&negative)->default_value(25),
     "number of negative samples per context word pair")
    ("clustered_input", po::bool_switch(&clustered_input)->default_value(false),
     "flag to utilize separate files for each server in a distributed setting")
    ("write_results", po::value<bool>(&write_results)->default_value(false), "write out results")
    ("localize_pos", po::value<bool>(&localize_positives)->default_value(true), "localize contextual data beforehand (default: yes)")
    ("sync_push", po::value<bool>(&sync_push)->default_value(true), "use synchronous pushes? (default: yes)")
    ("data_words", po::value<long long int>(&data_words)->default_value(numeric_limits<long long int>::max()), "use synchronous pushes? (default: yes)")
    ("min_count", po::value<int>(&min_count)->default_value(5), "learn embeddings for all words with count larger than min_cout")
    ("neg_power", po::value<double>(&neg_power)->default_value(0.75), "power for negative sampling")
    ("starting_alpha", po::value<ValT>(&alpha)->default_value(0.025), "starting alpha (learning rate)")
    ("subsample", po::value<ValT>(&sample)->default_value(1e-4), "subsample frequent words")
    ("replicate", po::value<size_t>(&replicate_n)->default_value(0), "number of parameters to replicate")
    ("peek_ahead", po::value<size_t>(&peek_ahead_n)->default_value(0), "peek ahead into the next sentence and localize N words. 0 to disable")
    ("prep_context_ahead", po::value<int>(&prep_context_ahead)->default_value(10), "How many contexts ahead of time to localize positives and prepare negatives")
    ("downsample_replicated_negatives", po::value<bool>(&downsample_replicated_negative_samples)->default_value(false), "downsample replicated negative samples according to the number of servers")
    ("binary", po::value<int>(&binary)->default_value(1), "output in binary human-readable format(default)")
    ("init_model", po::value<string>(&init_model)->default_value("random"), "how to init the parameters. options: 'none', 'random', '[PATH TO CHECKPOINT]'")
    ("model_seed", po::value<unsigned long long>(&model_seed)->default_value(134827), "seed for model generation")
    ("max_runtime", po::value<long>(&max_runtime)->default_value(std::numeric_limits<long>::max()), "set a maximum run tim, after which the job will be terminated (in seconds)")
    ;

  // add system options
  ServerT::AddSystemOptions(desc);

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

  // Warnings about nonsensical settings
  if (replicate_n != 0 && SamplingSupport<ValT, WorkerT>::sampling_strategy == SamplingSupportType::OnlyLocal &&
      !downsample_replicated_negative_samples) {
    ALOG("[WARNING] You are combining the 'only local' sampling strategy with replication, but downsampling of replicated parameters is turned off. This might deteriorate the sampling distribution even further.");
  }

  if (vm.count("help")) {
    cout << desc << "\n";
    return 1;
  }
  return 0;
}

int main(int argc, char *argv[]) {
  // Read cmd arguments
  int po_error = process_program_options(argc, argv);
  if (po_error) return 1;

  Postoffice::Get()->enable_dynamic_allocation(num_keys, num_threads);
  std::string role = std::string(getenv("DMLC_ROLE"));

  std::cout << "Word2vec: Starting " << role << ": running " << num_iterations << " iterations on " << num_keys
            << " keys in " << num_threads << " threads\n"
            << "embed_dim " << embed_dim << ", sync_push: " << sync_push << ", min_count: " << min_count << ", neg_power " << neg_power << "\n"
            << "file " << in_file << ", vocab: " << vocab_retrieve << ", data_words " << data_words << "\n"
            << "window " << window << ", negative " << negative << "\n"
            << "localize_pos " << localize_positives << ", " << "replicate_n " << replicate_n << "\n"
            << ReplicaManager<ValT,HandleT>::PrintOptions() << "\n";

  if (role.compare("scheduler") == 0) {
    Start(0);
    Finalize(0, true);
  } else if (role.compare("server") == 0) {

    // initialize data structures that are shared among threads
    init_shared_datastructures();
    starting_alpha = alpha;

    // replication
    vector<Key> hotspot_keys = determine_hotspot_keys(replicate_n);
    // std::sort(hotspot_keys.begin(), hotspot_keys.end());

    // construct unigram table (for efficient negative sampling)
    InitUnigramTable(hotspot_keys);

    // Start the server system
    int server_customer_id = 0; // server gets customer_id=0, workers 1..n
    Start(server_customer_id);
    HandleT handle(num_keys, embed_dim); // the handle specifies how the server handles incoming Push() and Pull() calls
    auto server = new ServerT(server_customer_id, handle, &hotspot_keys);
    RegisterExitCallback([server]() { delete server; });
    num_workers = ps::NumServers() * num_threads;

    // generate vocab out of the clustered files first, else nobody is on the same page
    CHECK(!clustered_input || (clustered_input && vocab_retrieve.size() > 0))
      << "ERROR_________clustered computing with several files requires a pregenerated vocab created out of those files,"
      << endl << " so create one first and relaunch accordingly" << endl;

    // make sure all servers are set up
    server->Barrier();

    // sampling support
    negs_next_random = (long long) 17+ps::MyRank()*123;
    server->enable_sampling_support(&DrawNegSample);

    // run worker(s)
    std::vector <std::thread> workers{};
    for (size_t i = 0; i != num_threads; ++i) {
      workers.push_back(std::thread(RunWorker, i + 1, server));
      std::string name = std::to_string(ps::MyRank())+"-worker-"+std::to_string(ps::MyRank()*num_threads + i);
      SET_THREAD_NAME((&workers[workers.size()-1]), name.c_str());
    }

    // wait for the workers to finish
    for (auto &w : workers)
      w.join();

    // stop the server
    server->shutdown();
    Finalize(server_customer_id, true);
  }
}