Upload everything

blossom.cpp

+#include "blossom.h"
+
+pair < double, vector < int > > solve_cpp ( Graph < double > g, vector < int > in_set ) {
+	int n = g.get_n (), m = g.get_m ();
+	vector<int> bts ( n + 2, 0 ), stb ( n + 2 );
+	int N = 0;
+	for ( auto x : in_set ) {
+		++N;
+		bts[x] = N;
+		stb[N] = x;
+	}
+	int M = 0;
+	vector<MWM::InputEdge> edges;
+	vector<int> ou ( N + 2 ), ov ( N + 2 );
+	for ( int i = 1; i <= N; ++i ) {
+		for ( auto v : g.getedges ( stb[i] ) ) {
+			if ( bts[v.v] && i < bts[v.v] ) {
+				edges.push_back ( { i, bts[v.v], 500000000 - static_cast < int > ( round ( v.w * 100000000 ) ) } );
+				ou[i + 1] += 1; ov[bts[v.v] + 1] += 1;
+				++M;
+			}
+		}
+	}
+	edges.resize ( M * 2 );
+	for ( int i = 1; i <= N + 1; ++i ) ov[i] += ov[i - 1];
+	for ( int i = 0; i < M; ++i ) edges[M + ( ov[edges[i].to]++ )] = edges[i];
+	for ( int i = 1; i <= N + 1; ++i ) ou[i] += ou[i - 1];
+	for ( int i = 0; i < M; ++i ) edges[ou[edges[i + M].from]++] = edges[i + M];
+	edges.resize ( M );
+	auto ans = MWM ( N, edges ).maximum_weighted_matching ();
+	vector < int > method ( n + 2 );
+	for ( int i = 1; i <= N; ++i ) method[stb[i]] = stb[ans.second[i]];
+	return make_pair ( ( 500000000LL * ( N / 2 ) - ans.first ) / 100000000.0, method );
+}
\ No newline at end of file

blossom.h

+#pragma once
+
+#include <cstdio>
+#include <cassert>
+#include <cmath>
+
+#include <algorithm>
+#include <vector>
+#include <set>
+#include <queue>
+#include <functional>
+
+#include "graph.h"
+
+using namespace std;
+
+template <typename CostType, typename TotalCostType = int64_t>
+class MaximumWeightedMatching
+{
+	/*
+	Maximum Weighted Matching in General Graphs.
+	- O(nm log(n)) time
+	- O(n + m) space
+
+	Note: each vertex is 1-indexed.
+	*/
+public:
+	using cost_t = CostType;
+	using tcost_t = TotalCostType;
+
+private:
+	enum Label { kSeparated = -2, kInner = -1, kFree = 0, kOuter = 1 };
+	static constexpr cost_t Inf = cost_t ( 1 ) << ( sizeof ( cost_t ) * 8 - 2 );
+
+private:
+	template <typename T>
+	class BinaryHeap
+	{
+	public:
+		struct Node
+		{
+			bool operator < ( const Node& rhs ) const { return value < rhs.value; }
+			T value; int id;
+		};
+		BinaryHeap () {}
+		BinaryHeap ( int N ) : size_ ( 0 ), node ( N + 1 ), index ( N, 0 ) {}
+		int size () const { return size_; }
+		bool empty () const { return size_ == 0; }
+		void clear () { while ( size_ > 0 ) index[node[size_--].id] = 0; }
+		T min () const { return node[1].value; }
+		int argmin () const { return node[1].id; } // argmin ?
+		T get_val ( int id ) const { return node[index[id]].value; }
+		void pop () { if ( size_ > 0 ) pop ( 1 ); }
+		void erase ( int id ) { if ( index[id] ) pop ( index[id] ); }
+		bool has ( int id ) const { return index[id] != 0; }
+		void update ( int id, T v ) {
+			if ( !has ( id ) ) return push ( id, v );
+			bool up = ( v < node[index[id]].value );
+			node[index[id]].value = v;
+			if ( up ) up_heap ( index[id] );
+			else down_heap ( index[id] );
+		}
+		void decrease_key ( int id, T v ) {
+			if ( !has ( id ) ) return push ( id, v );
+			if ( v < node[index[id]].value ) node[index[id]].value = v, up_heap ( index[id] );
+		}
+		void push ( int id, T v ) {
+			// assert(!has(id));
+			index[id] = ++size_; node[size_] = { v, id };
+			up_heap ( size_ );
+		}
+	private:
+		void pop ( int pos ) {
+			index[node[pos].id] = 0;
+			if ( pos == size_ ) { --size_; return; }
+			bool up = ( node[size_].value < node[pos].value );
+			node[pos] = node[size_--]; index[node[pos].id] = pos;
+			if ( up ) up_heap ( pos );
+			else down_heap ( pos );
+		}
+		void swap_node ( int a, int b ) {
+			swap ( node[a], node[b] ); index[node[a].id] = a; index[node[b].id] = b;
+		}
+		void down_heap ( int pos ) {
+			for ( int k = pos, nk = k; 2 * k <= size_; k = nk ) {
+				if ( node[2 * k] < node[nk] ) nk = 2 * k;
+				if ( 2 * k + 1 <= size_ && node[2 * k + 1] < node[nk] ) nk = 2 * k + 1;
+				if ( nk == k ) break;
+				swap_node ( k, nk );
+			}
+		}
+		void up_heap ( int pos ) {
+			for ( int k = pos; k > 1 && node[k] < node[k >> 1]; k >>= 1 ) swap_node ( k, k >> 1 );
+		}
+		int size_;
+		vector<Node> node;
+		vector<int> index;
+	};
+
+	template <typename Key>
+	class PairingHeaps
+	{
+	private:
+		struct Node
+		{
+			Node () : prev ( -1 ) {} // "prev < 0" means the node is unused.
+			Node ( Key v ) : key ( v ), child ( 0 ), next ( 0 ), prev ( 0 ) {}
+			Key key; int child, next, prev;
+		};
+	public:
+		PairingHeaps ( int H, int N ) : heap ( H ), node ( N ) {
+			// It consists of `H` Pairing heaps.
+			// Each heap-node ID can appear at most 1 time(s) among heaps
+			// and should be in [1, N).
+		}
+
+		void clear ( int h ) { if ( heap[h] ) clear_rec ( heap[h] ), heap[h] = 0; }
+		void clear_all () {
+			for ( size_t i = 0; i < heap.size (); ++i ) heap[i] = 0;
+			for ( size_t i = 0; i < node.size (); ++i ) node[i] = Node ();
+		}
+		bool empty ( int h ) const { return !heap[h]; }
+		bool used ( int v ) const { return node[v].prev >= 0; }
+		Key min ( int h ) const { return node[heap[h]].key; }
+		int argmin ( int h ) const { return heap[h]; }
+
+		void pop ( int h ) {
+			// assert(!empty(h));
+			erase ( h, heap[h] );
+		}
+		void push ( int h, int v, Key key ) {
+			// assert(!used(v));
+			node[v] = Node ( key );
+			heap[h] = merge ( heap[h], v );
+		}
+		void erase ( int h, int v ) {
+			if ( !used ( v ) ) return;
+			int w = two_pass_pairing ( node[v].child );
+			if ( !node[v].prev ) heap[h] = w;
+			else {
+				cut ( v );
+				heap[h] = merge ( heap[h], w );
+			}
+			node[v].prev = -1;
+		}
+		void decrease_key ( int h, int v, Key key ) {
+			if ( !used ( v ) ) return push ( h, v, key );
+			if ( !node[v].prev ) node[v].key = key;
+			else {
+				cut ( v ); node[v].key = key;
+				heap[h] = merge ( heap[h], v );
+			}
+		}
+
+	private:
+		void clear_rec ( int v ) {
+			for ( ; v; v = node[v].next ) {
+				if ( node[v].child ) clear_rec ( node[v].child );
+				node[v].prev = -1;
+			}
+		}
+
+		inline void cut ( int v ) {
+			auto& n = node[v]; int pv = n.prev, nv = n.next;
+			auto& pn = node[pv];
+			if ( pn.child == v ) pn.child = nv;
+			else pn.next = nv;
+			node[nv].prev = pv;
+			n.next = n.prev = 0;
+		}
+
+		int merge ( int l, int r ) {
+			if ( !l ) return r;
+			if ( !r ) return l;
+			if ( node[l].key > node[r].key ) swap ( l, r );
+			int lc = node[r].next = node[l].child;
+			node[l].child = node[lc].prev = r;
+			return node[r].prev = l;
+		}
+
+		int two_pass_pairing ( int root ) {
+			if ( !root ) return 0;
+			int a = root; root = 0;
+			while ( a ) {
+				int b = node[a].next, na = 0;
+				node[a].prev = node[a].next = 0;
+				if ( b ) na = node[b].next, node[b].prev = node[b].next = 0;
+				a = merge ( a, b );
+				node[a].next = root; root = a; a = na;
+			}
+			int s = node[root].next; node[root].next = 0;
+			while ( s ) {
+				int t = node[s].next; node[s].next = 0;
+				root = merge ( root, s );
+				s = t;
+			}
+			return root;
+		}
+
+	private:
+		vector<int> heap;
+		vector<Node> node;
+	};
+
+	template <typename T>
+	struct PriorityQueue : public priority_queue< T, vector<T>, greater<T> >
+	{
+		PriorityQueue () {}
+		PriorityQueue ( int N ) { this->c.reserve ( N ); }
+		T min () { return this->top (); }
+		void clear () { this->c.clear (); }
+	};
+
+	template <typename T>
+	struct Queue
+	{
+		Queue () {}
+		Queue ( int N ) : qh ( 0 ), qt ( 0 ), data ( N ) {}
+		T operator [] ( int i ) const { return data[i]; }
+		void enqueue ( int u ) { data[qt++] = u; }
+		int dequeue () { return data[qh++]; }
+		bool empty () const { return qh == qt; }
+		void clear () { qh = qt = 0; }
+		int size () const { return qt; }
+		int qh, qt;
+		vector<T> data;
+	};
+
+public:
+	struct InputEdge { int from, to; cost_t cost; };
+
+private:
+	template <typename T> using ModifiableHeap = BinaryHeap<T>;
+	template <typename T> using ModifiableHeaps = PairingHeaps<T>;
+	template <typename T> using FastHeap = PriorityQueue<T>;
+
+	struct Edge { int to; cost_t cost; };
+	struct Link { int from, to; };
+	struct Node
+	{
+		struct NodeLink { int b, v; };
+		Node () {}
+		Node ( int u ) : parent ( 0 ), size ( 1 ) { link[0] = link[1] = { u, u }; }
+		int next_v () const { return link[0].v; }
+		int next_b () const { return link[0].b; }
+		int prev_v () const { return link[1].v; }
+		int prev_b () const { return link[1].b; }
+		int parent, size;
+		NodeLink link[2];
+	};
+	struct Event
+	{
+		Event () {}
+		Event ( cost_t time, int id ) : time ( time ), id ( id ) {}
+		bool operator < ( const Event& rhs ) const { return time < rhs.time; }
+		bool operator > ( const Event& rhs ) const { return time > rhs.time; }
+		cost_t time; int id;
+	};
+	struct EdgeEvent
+	{
+		EdgeEvent () {}
+		EdgeEvent ( cost_t time, int from, int to ) : time ( time ), from ( from ), to ( to ) {}
+		bool operator > ( const EdgeEvent& rhs ) const { return time > rhs.time; }
+		bool operator < ( const EdgeEvent& rhs ) const { return time < rhs.time; }
+		cost_t time; int from, to;
+	};
+
+public:
+	MaximumWeightedMatching ( int N, const vector<InputEdge>& in )
+		: N ( N ), B ( ( N - 1 ) / 2 ), S ( N + B + 1 ), ofs ( N + 2 ), edges ( in.size () * 2 ),
+		heap2 ( S ), heap2s ( S, S ), heap3 ( edges.size () ), heap4 ( S ) {
+
+		for ( auto& e : in ) ofs[e.from + 1]++, ofs[e.to + 1]++;
+		for ( int i = 1; i <= N + 1; ++i ) ofs[i] += ofs[i - 1];
+		for ( auto& e : in ) {
+			edges[ofs[e.from]++] = { e.to, e.cost * 2 };
+			edges[ofs[e.to]++] = { e.from, e.cost * 2 };
+		}
+		for ( int i = N + 1; i > 0; --i ) ofs[i] = ofs[i - 1];
+		ofs[0] = 0;
+	}
+
+	pair< tcost_t, vector<int> > maximum_weighted_matching ( bool init_matching = false ) {
+		initialize ();
+		set_potential ();
+		if ( init_matching ) find_maximal_matching ();
+		for ( int u = 1; u <= N; ++u ) if ( !mate[u] ) do_edmonds_search ( u );
+		tcost_t ret = compute_optimal_value ();
+		return make_pair ( ret, mate );
+	}
+
+private:
+	tcost_t compute_optimal_value () const {
+		tcost_t ret = 0;
+		for ( int u = 1; u <= N; ++u ) if ( mate[u] > u ) {
+			cost_t max_c = 0;
+			for ( int eid = ofs[u]; eid < ofs[u + 1]; ++eid ) {
+				if ( edges[eid].to == mate[u] ) max_c = max ( max_c, edges[eid].cost );
+			}
+			ret += max_c;
+		}
+		return ret >> 1;
+	}
+
+	inline tcost_t reduced_cost ( int u, int v, const Edge& e ) const {
+		return tcost_t ( potential[u] ) + potential[v] - e.cost;
+	}
+
+	void rematch ( int v, int w ) {
+		int t = mate[v]; mate[v] = w;
+		if ( mate[t] != v ) return;
+		if ( link[v].to == surface[link[v].to] ) {
+			mate[t] = link[v].from;
+			rematch ( mate[t], t );
+		} else {
+			int x = link[v].from, y = link[v].to;
+			rematch ( x, y ); rematch ( y, x );
+		}
+	}
+
+	void fix_mate_and_base ( int b ) {
+		if ( b <= N ) return;
+		int bv = base[b], mv = node[bv].link[0].v, bmv = node[bv].link[0].b;
+		int d = ( node[bmv].link[1].v == mate[mv] ) ? 0 : 1;
+		while ( 1 ) {
+			int mv = node[bv].link[d].v, bmv = node[bv].link[d].b;
+			if ( node[bmv].link[1 ^ d].v != mate[mv] ) break;
+			fix_mate_and_base ( bv ); fix_mate_and_base ( bmv );
+			bv = node[bmv].link[d].b;
+		}
+		fix_mate_and_base ( base[b] = bv );
+		mate[b] = mate[bv];
+	}
+
+	void reset_time () {
+		time_current_ = 0; event1 = { Inf, 0 };
+	}
+
+	void reset_blossom ( int b ) {
+		label[b] = kFree; link[b].from = 0; slack[b] = Inf; lazy[b] = 0;
+	}
+
+	void reset_all () {
+		label[0] = kFree; link[0].from = 0;
+		for ( int v = 1; v <= N; ++v ) { // should be optimized for sparse graphs.
+			if ( label[v] == kOuter ) potential[v] -= time_current_;
+			else {
+				int bv = surface[v];
+				potential[v] += lazy[bv];
+				if ( label[bv] == kInner ) potential[v] += time_current_ - time_created[bv];
+			}
+			reset_blossom ( v );
+		}
+		for ( int b = N + 1, r = B - unused_bid_idx_; r > 0 && b < S; ++b ) if ( base[b] != b ) {
+			if ( surface[b] == b ) {
+				fix_mate_and_base ( b );
+				if ( label[b] == kOuter ) potential[b] += ( time_current_ - time_created[b] ) << 1;
+				else if ( label[b] == kInner ) fix_blossom_potential<kInner> ( b );
+				else fix_blossom_potential<kFree> ( b );
+			}
+			heap2s.clear ( b );
+			reset_blossom ( b ); --r;
+		}
+
+		que.clear ();
+		reset_time (); heap2.clear ();
+		heap3.clear (); heap4.clear ();
+	}
+
+	void do_edmonds_search ( int root ) {
+		if ( potential[root] == 0 ) return;
+		link_blossom ( surface[root], { 0, 0 } );
+		push_outer_and_fix_potentials ( surface[root], 0 );
+		for ( bool augmented = false; !augmented; ) {
+			augmented = augment ( root );
+			if ( augmented ) break;
+			augmented = adjust_dual_variables ( root );
+		}
+		reset_all ();
+	}
+
+	template <Label Lab>
+	inline cost_t fix_blossom_potential ( int b ) {
+		// Return the amount.
+		// (If v is an atom, the potential[v] will not be changed.)
+		cost_t d = lazy[b]; lazy[b] = 0;
+		if ( Lab == kInner ) {
+			cost_t dt = time_current_ - time_created[b];
+			if ( b > N ) potential[b] -= dt << 1;
+			d += dt;
+		}
+		return d;
+	}
+
+	template <Label Lab>
+	inline void update_heap2 ( int x, int y, int by, cost_t t ) {
+		if ( t >= slack[y] ) return;
+		slack[y] = t; best_from[y] = x;
+		if ( y == by ) {
+			if ( Lab != kInner ) heap2.decrease_key ( y, EdgeEvent ( t + lazy[y], x, y ) );
+		} else {
+			int gy = group[y];
+			if ( gy != y ) {
+				if ( t >= slack[gy] ) return;
+				slack[gy] = t;
+			}
+			heap2s.decrease_key ( by, gy, EdgeEvent ( t, x, y ) );
+			if ( Lab == kInner ) return;
+			EdgeEvent m = heap2s.min ( by );
+			heap2.decrease_key ( by, EdgeEvent ( m.time + lazy[by], m.from, m.to ) );
+		}
+	}
+
+	void activate_heap2_node ( int b ) {
+		if ( b <= N ) {
+			if ( slack[b] < Inf ) heap2.push ( b, EdgeEvent ( slack[b] + lazy[b], best_from[b], b ) );
+		} else {
+			if ( heap2s.empty ( b ) ) return;
+			EdgeEvent m = heap2s.min ( b );
+			heap2.push ( b, EdgeEvent ( m.time + lazy[b], m.from, m.to ) );
+		}
+	}
+
+	void swap_blossom ( int a, int b ) {
+		// Assume that `b` is a maximal blossom.
+		swap ( base[a], base[b] ); if ( base[a] == a ) base[a] = b;
+		swap ( heavy[a], heavy[b] ); if ( heavy[a] == a ) heavy[a] = b;
+		swap ( link[a], link[b] );
+		swap ( mate[a], mate[b] );
+		swap ( potential[a], potential[b] ); swap ( lazy[a], lazy[b] );
+		swap ( time_created[a], time_created[b] );
+		for ( int d = 0; d < 2; ++d ) node[node[a].link[d].b].link[1 ^ d].b = b;
+		swap ( node[a], node[b] );
+	}
+
+	void set_surface_and_group ( int b, int sf, int g ) {
+		surface[b] = sf, group[b] = g;
+		if ( b <= N ) return;
+		for ( int bb = base[b]; surface[bb] != sf; bb = node[bb].next_b () ) {
+			set_surface_and_group ( bb, sf, g );
+		}
+	}
+
+	void merge_smaller_blossoms ( int bid ) {
+		int lb = bid, largest_size = 1;
+		for ( int beta = base[bid], b = beta; ;) {
+			if ( node[b].size > largest_size ) largest_size = node[b].size, lb = b;
+			if ( ( b = node[b].next_b () ) == beta ) break;
+		}
+		for ( int beta = base[bid], b = beta; ;) {
+			if ( b != lb ) set_surface_and_group ( b, lb, b );
+			if ( ( b = node[b].next_b () ) == beta ) break;
+		}
+		group[lb] = lb;
+		if ( largest_size > 1 ) {
+			surface[bid] = heavy[bid] = lb;
+			swap_blossom ( lb, bid );
+		} else heavy[bid] = 0;
+	}
+
+	void contract ( int x, int y, int eid ) {
+		int bx = surface[x], by = surface[y]; assert ( bx != by );
+		const int h = -( eid + 1 );
+		link[surface[mate[bx]]].from = link[surface[mate[by]]].from = h;
+
+		int lca = -1;
+		while ( 1 ) {
+			if ( mate[by] != 0 ) swap ( bx, by );
+			bx = lca = surface[link[bx].from];
+			if ( link[surface[mate[bx]]].from == h ) break;
+			link[surface[mate[bx]]].from = h;
+		}
+
+		const int bid = unused_bid[--unused_bid_idx_]; assert ( unused_bid_idx_ >= 0 );
+		int tree_size = 0;
+		for ( int d = 0; d < 2; ++d ) {
+			for ( int bv = surface[x]; bv != lca; ) {
+				int mv = mate[bv], bmv = surface[mv], v = mate[mv];
+				int f = link[v].from, t = link[v].to;
+				tree_size += node[bv].size + node[bmv].size;
+				link[mv] = { x, y };
+
+				if ( bv > N ) potential[bv] += ( time_current_ - time_created[bv] ) << 1;
+				if ( bmv > N ) heap4.erase ( bmv );
+				push_outer_and_fix_potentials ( bmv, fix_blossom_potential<kInner> ( bmv ) );
+
+				node[bv].link[d] = { bmv, mv };
+				node[bmv].link[1 ^ d] = { bv, v }; node[bmv].link[d] = { bv = surface[f], f };
+				node[bv].link[1 ^ d] = { bmv, t };
+			}
+			node[surface[x]].link[1 ^ d] = { surface[y], y };
+			swap ( x, y );
+		}
+		if ( lca > N ) potential[lca] += ( time_current_ - time_created[lca] ) << 1;
+		node[bid].size = tree_size + node[lca].size;
+		base[bid] = lca; link[bid] = link[lca]; mate[bid] = mate[lca];
+		label[bid] = kOuter;
+		surface[bid] = bid; time_created[bid] = time_current_;
+		potential[bid] = 0; lazy[bid] = 0;
+
+		merge_smaller_blossoms ( bid ); // O(n log n) time / Edmonds search
+	}
+
+	void link_blossom ( int v, Link l ) {
+		link[v] = { l.from, l.to };
+		if ( v <= N ) return;
+		int b = base[v]; link_blossom ( b, l );
+		int pb = node[b].prev_b ();
+		l = { node[pb].next_v (), node[b].prev_v () };
+		for ( int bv = b; ; ) {
+			int bw = node[bv].next_b ();
+			if ( bw == b ) break;
+			link_blossom ( bw, l );
+			Link nl = { node[bw].prev_v (), node[bv].next_v () };
+			bv = node[bw].next_b ();
+			link_blossom ( bv, nl );
+		}
+	}
+
+	void push_outer_and_fix_potentials ( int v, cost_t d ) {
+		label[v] = kOuter;
+		if ( v > N ) {
+			for ( int b = base[v]; label[b] != kOuter; b = node[b].next_b () ) {
+				push_outer_and_fix_potentials ( b, d );
+			}
+		} else {
+			potential[v] += time_current_ + d;
+			if ( potential[v] < event1.time ) event1 = { potential[v], v };
+			que.enqueue ( v );
+		}
+	}
+
+	bool grow ( int root, int x, int y ) {
+		int by = surface[y];
+		bool visited = ( label[by] != kFree );
+		if ( !visited ) link_blossom ( by, { 0, 0 } );
+		label[by] = kInner; time_created[by] = time_current_; heap2.erase ( by );
+		if ( y != by ) heap4.update ( by, time_current_ + ( potential[by] >> 1 ) );
+		int z = mate[by];
+		if ( z == 0 && by != surface[root] ) {
+			rematch ( x, y ); rematch ( y, x );
+			return true;
+		}
+		int bz = surface[z];
+		if ( !visited ) link_blossom ( bz, { x, y } );
+		else link[bz] = link[z] = { x, y };
+		push_outer_and_fix_potentials ( bz, fix_blossom_potential<kFree> ( bz ) );
+		time_created[bz] = time_current_; heap2.erase ( bz );
+		return false;
+	}
+
+	void free_blossom ( int bid ) {
+		unused_bid[unused_bid_idx_++] = bid;
+		base[bid] = bid;
+	}
+
+	int recalculate_minimum_slack ( int b, int g ) {
+		// Return the destination of the best edge of blossom `g`.
+		if ( b <= N ) {
+			if ( slack[b] >= slack[g] ) return 0;
+			slack[g] = slack[b]; best_from[g] = best_from[b];
+			return b;
+		}
+		int v = 0;
+		for ( int beta = base[b], bb = beta; ; ) {
+			int w = recalculate_minimum_slack ( bb, g );
+			if ( w != 0 ) v = w;
+			if ( ( bb = node[bb].next_b () ) == beta ) break;
+		}
+		return v;
+	}
+
+	void construct_smaller_components ( int b, int sf, int g ) {
+		surface[b] = sf, group[b] = g; // `group[b] = g` is unneeded.
+		if ( b <= N ) return;
+		for ( int bb = base[b]; surface[bb] != sf; bb = node[bb].next_b () ) {
+			if ( bb == heavy[b] ) {
+				construct_smaller_components ( bb, sf, g );
+			} else {
+				set_surface_and_group ( bb, sf, bb );
+				int to = 0;
+				if ( bb > N ) slack[bb] = Inf, to = recalculate_minimum_slack ( bb, bb );
+				else if ( slack[bb] < Inf ) to = bb;
+				if ( to > 0 ) heap2s.push ( sf, bb, EdgeEvent ( slack[bb], best_from[bb], to ) );
+			}
+		}
+	}
+
+	void move_to_largest_blossom ( int bid ) {
+		const int h = heavy[bid];
+		cost_t d = ( time_current_ - time_created[bid] ) + lazy[bid]; lazy[bid] = 0;
+		for ( int beta = base[bid], b = beta; ;) {
+			time_created[b] = time_current_;
+			lazy[b] = d;
+			if ( b != h ) construct_smaller_components ( b, b, b ), heap2s.erase ( bid, b );
+			if ( ( b = node[b].next_b () ) == beta ) break;
+		}
+		if ( h > 0 ) swap_blossom ( h, bid ), bid = h;
+		free_blossom ( bid );
+	}
+
+	void expand ( int bid ) {
+		int mv = mate[base[bid]];
+		move_to_largest_blossom ( bid ); // O(n log n) time / Edmonds search
+		Link old_link = link[mv];
+		int old_base = surface[mate[mv]], root = surface[old_link.to];
+		int d = ( mate[root] == node[root].link[0].v ) ? 1 : 0;
+		for ( int b = node[old_base].link[d ^ 1].b; b != root; ) {
+			label[b] = kSeparated; activate_heap2_node ( b ); b = node[b].link[d ^ 1].b;
+			label[b] = kSeparated; activate_heap2_node ( b ); b = node[b].link[d ^ 1].b;
+		}
+		for ( int b = old_base; ; b = node[b].link[d].b ) {
+			label[b] = kInner;
+			int nb = node[b].link[d].b;
+			if ( b == root ) link[mate[b]] = old_link;
+			else link[mate[b]] = { node[b].link[d].v, node[nb].link[d ^ 1].v };
+			link[surface[mate[b]]] = link[mate[b]]; // fix tree links
+			if ( b > N ) {
+				if ( potential[b] == 0 ) expand ( b );
+				else heap4.push ( b, time_current_ + ( potential[b] >> 1 ) );
+			}
+			if ( b == root ) break;
+			push_outer_and_fix_potentials ( nb, fix_blossom_potential<kInner> ( b = nb ) );
+		}
+	}
+
+	bool augment ( int root ) {
+		// Return true if an augmenting path is found.
+		while ( !que.empty () ) {
+			int x = que.dequeue (), bx = surface[x];
+			if ( potential[x] == time_current_ ) {
+				if ( x != root ) rematch ( x, 0 );
+				return true;
+			}
+			for ( int eid = ofs[x]; eid < ofs[x + 1]; ++eid ) {
+				auto& e = edges[eid]; int y = e.to, by = surface[y];
+				if ( bx == by ) continue;
+				Label l = label[by];
+				if ( l == kOuter ) {
+					cost_t t = reduced_cost ( x, y, e ) >> 1; // < 2 * Inf
+					if ( t == time_current_ ) {
+						contract ( x, y, eid ); bx = surface[x];
+					} else if ( t < event1.time ) {
+						heap3.emplace ( t, x, eid );
+					}
+				} else {
+					tcost_t t = reduced_cost ( x, y, e ); // < 3 * Inf
+					if ( t >= Inf ) continue;
+					if ( l != kInner ) {
+						if ( cost_t ( t ) + lazy[by] == time_current_ ) {
+							if ( grow ( root, x, y ) ) return true;
+						} else update_heap2<kFree> ( x, y, by, t );
+					} else {
+						if ( mate[x] != y ) update_heap2<kInner> ( x, y, by, t );
+					}
+				}
+			}
+		}
+		return false;
+	}
+
+	bool adjust_dual_variables ( int root ) {
+		// delta1 : rematch
+		cost_t time1 = event1.time;
+
+		// delta2 : grow
+		cost_t time2 = Inf;
+		if ( !heap2.empty () ) time2 = heap2.min ().time;
+
+		// delta3 : contract : O(m log n) time / Edmonds search [ bottleneck (?) ]
+		cost_t time3 = Inf;
+		while ( !heap3.empty () ) {
+			EdgeEvent e = heap3.min ();
+			int x = e.from, y = edges[e.to].to; // e.to is some edge id.
+			if ( surface[x] != surface[y] ) {
+				time3 = e.time;
+				break;
+			} else heap3.pop ();
+		}
+
+		// delta4 : expand
+		cost_t time4 = Inf;
+		if ( !heap4.empty () ) time4 = heap4.min ();
+
+		// -- events --
+		cost_t time_next = min ( min ( time1, time2 ), min ( time3, time4 ) );
+		assert ( time_current_ <= time_next && time_next < Inf );
+		time_current_ = time_next;
+
+		if ( time_current_ == event1.time ) {
+			int x = event1.id;
+			if ( x != root ) rematch ( x, 0 );
+			return true;
+		}
+		while ( !heap2.empty () && heap2.min ().time == time_current_ ) {
+			int x = heap2.min ().from, y = heap2.min ().to;
+			if ( grow ( root, x, y ) ) return true; // `grow` function will call `heap2.erase(by)`.
+		}
+		while ( !heap3.empty () && heap3.min ().time == time_current_ ) {
+			int x = heap3.min ().from, eid = heap3.min ().to;
+			int y = edges[eid].to; heap3.pop ();
+			if ( surface[x] == surface[y] ) continue;
+			contract ( x, y, eid );
+		}
+		while ( !heap4.empty () && heap4.min () == time_current_ ) {
+			int b = heap4.argmin (); heap4.pop ();
+			expand ( b );
+		}
+		return false;
+	}
+
+private:
+	void initialize () {
+		que = Queue<int> ( N );
+		mate.assign ( S, 0 );
+		link.assign ( S, { 0, 0 } );
+		label.assign ( S, kFree );
+		base.resize ( S ); for ( int u = 1; u < S; ++u ) base[u] = u;
+		surface.resize ( S ); for ( int u = 1; u < S; ++u ) surface[u] = u;
+
+		potential.resize ( S );
+		node.resize ( S ); for ( int b = 1; b < S; ++b ) node[b] = Node ( b );
+
+		unused_bid.resize ( B ); for ( int i = 0; i < B; ++i ) unused_bid[i] = N + B - i;
+		unused_bid_idx_ = B;
+
+		// for O(nm log n) implementation
+		reset_time ();
+		time_created.resize ( S );
+		slack.resize ( S ); for ( int i = 0; i < S; ++i ) slack[i] = Inf;
+		best_from.assign ( S, 0 );
+		heavy.assign ( S, 0 );
+		lazy.assign ( S, 0 );
+		group.resize ( S ); for ( int i = 0; i < S; ++i ) group[i] = i;
+	}
+
+	void set_potential () {
+		for ( int u = 1; u <= N; ++u ) {
+			cost_t max_c = 0;
+			for ( int eid = ofs[u]; eid < ofs[u + 1]; ++eid ) {
+				max_c = max ( max_c, edges[eid].cost );
+			}
+			potential[u] = max_c >> 1;
+		}
+	}
+
+	void find_maximal_matching () {
+		// Find a maximal matching naively.
+		for ( int u = 1; u <= N; ++u ) if ( !mate[u] ) {
+			for ( int eid = ofs[u]; eid < ofs[u + 1]; ++eid ) {
+				auto& e = edges[eid]; int v = e.to;
+				if ( mate[v] > 0 || reduced_cost ( u, v, e ) > 0 ) continue;
+				mate[u] = v; mate[v] = u;
+				break;
+			}
+		}
+	}
+
+private:
+	int N, B, S; // N = |V|, B = (|V| - 1) / 2, S = N + B + 1
+	vector<int> ofs;
+	vector<Edge> edges;
+
+	Queue<int> que;
+	vector<int> mate, surface, base;
+	vector<Link> link;
+	vector<Label> label;
+	vector<cost_t> potential;
+
+	vector<int> unused_bid; int unused_bid_idx_;
+	vector<Node> node;
+
+	// for O(nm log n) implementation
+	vector<int> heavy, group;
+	vector<cost_t> time_created, lazy, slack;
+	vector<int> best_from;
+
+	cost_t time_current_;
+	Event event1;
+	ModifiableHeap<EdgeEvent> heap2;
+	ModifiableHeaps<EdgeEvent> heap2s;
+	FastHeap<EdgeEvent> heap3;
+	ModifiableHeap<cost_t> heap4;
+};
+
+using MWM = MaximumWeightedMatching<int>;
+
+pair < double, vector < int > > solve_cpp ( Graph < double > g, vector < int > in_set );

compile.sh

+#!/bin/bash
+
+g++ main.cpp consts.cpp blossom.cpp simulate_annealing.cpp -o main -O2 -std=c++11 -lpthread

consts.cpp

+#include "consts.h"
+
+const double EPS = 1e-9;
+
+const double WEIGHT_ROAD = 1.5;
+const double WEIGHT_TRAIL = 1;
+
+const double WEIGHT_PARK = 0.03;
+const double WEIGHT_PICNIC = 0.06;
+const double WEIGHT_RESTROOM = 0.05;
+const double WEIGHT_CAMP = 0.06;

consts.h

+#pragma once
+
+extern const double EPS;
+
+extern const double WEIGHT_ROAD;
+extern const double WEIGHT_TRAIL;
+
+extern const double WEIGHT_PARK;
+extern const double WEIGHT_PICNIC;
+extern const double WEIGHT_RESTROOM;
+extern const double WEIGHT_CAMP;

data/1/desc.txt

+图片宽度为1813像素，长度为1312像素
+左上角为原点，先列后行（如右下角坐标为(1813,1312))
+图上298像素代表1/4 miles
+*说明：P-Parking, A-Picnic Area, R-Restrooms, B-Boat Ramp, V-Pavilion
+
+P 422 639
+P 319 871
+P 400 971
+P 551 917
+P 808 650
+P 756 445
+P 527 461
+P 582 416
+P 609 561
+P 1153 589
+P 1177 607
+P 1112 716
+P 1035 783
+P 884 854
+P 1127 1084
+A 383 906
+A 726 898
+A 744 852
+A 811 677
+A 704 559
+A 620 418
+A 807 438
+A 1129 576
+A 1204 665
+A 1111 745
+A 1037 845
+A 1328 945
+A 1265 1033
+A 1290 1104
+A 1147 1166
+A 1056 1180
+A 401 688
+R 455 975
+R 683 871
+R 635 580
+R 773 465
+R 1154 566
+R 1069 1047
+B 561 944
+V 674 549
+G 1193 1230

data/1/point.txt

+P 422 639
+P 319 871
+P 400 971
+P 551 917
+P 808 650
+P 756 445
+P 527 461
+P 582 416
+P 609 561
+P 1153 589
+P 1177 607
+P 1112 716
+P 1035 783
+P 884 854
+P 1127 1084
+A 383 906
+A 726 898
+A 744 852
+A 811 677
+A 704 559
+A 620 418
+A 807 438
+A 1129 576
+A 1204 665
+A 1111 745
+A 1037 845
+A 1328 945
+A 1265 1033
+A 1290 1104
+A 1147 1166
+A 1056 1180
+A 401 688
+R 455 975
+R 683 871
+R 635 580
+R 773 465
+R 1154 566
+R 1069 1047
+B 561 944
+V 674 549
+G 1193 1230

data/2/desc.txt

+图片宽度为2542像素，长度为3274像素
+左上角为原点，先列后行（如右下角坐标为(2542,3274))
+图上424像素代表1 miles
+*说明：P-Parking, F-Fire Station, C-Campground, B-Bunkers
+
+P 2031 917
+P 808 2571
+P 2257 1167
+P 1171 284
+P 1717 2195
+F 2210 1445
+F 703 2835
+C 719 2553
+B 1208 710
+G 2259 1171

data/2/point.txt

+P 2031 917
+P 808 2571
+P 2257 1167
+P 1171 284
+P 1717 2195
+F 2210 1445
+F 703 2835
+C 719 2553
+B 1208 710
+G 2259 1171

getpath.cpp

+#include <bits/stdc++.h>
+using namespace std;
+
+#define LL long long
+#define LD long double
+#define SC(t,x) static_cast<t>(x)
+#define AR(t) vector < t >
+#define PII pair < int, int >
+#define PLL pair < LL, LL >
+#define PIL pair < int, LL >
+#define PLI pair < LL, int >
+#define PQ priority_queue
+#define GT(type) greater < type >
+#define MP make_pair
+#define PB push_back
+#define PF push_front
+#define POB pop_back
+#define POF pop_front
+#define PRF first
+#define PRS second
+#define INIT(ar,val) memset ( ar, val, sizeof ( ar ) )
+#define lp(loop,start,end) for ( int loop = start; loop < end; ++loop )
+#define lpd(loop,start,end) for ( int loop = start; loop > end; --loop )
+#define lpi(loop,start,end) for ( int loop = start; loop <= end; ++loop )
+#define lpdi(loop,start,end) for ( int loop = start; loop >= end; --loop )
+#define qmax(a,b) (((a)>(b))?(a):(b))
+#define qmin(a,b) (((a)<(b))?(a):(b))
+#define qabs(a) (((a)>=0)?(a):(-(a)))
+
+const int MAXN = 100007;
+const int MAXM = MAXN * 2;
+
+struct eT
+{
+	void setd ( int _u, int _v, int _l )
+	{
+		u = _u, v = _v, last = _l;
+	}
+	
+	int u, v, last;
+}edge[MAXM*2];
+
+int n, m;
+int ke, la[MAXN];
+int id[MAXM];
+
+int deg[MAXN];
+int ans[MAXM], ka;
+bool vis[MAXM];
+
+int getdist ( pair < int, int > x, pair < int, int > y )
+{
+	return ( x.first - y.first ) * ( x.first - y.first ) + ( x.second - y.second ) * ( x.second - y.second );
+}
+
+void dfs ( int now )
+{
+	for ( int i = la[now]; ~i; i = la[now] ){
+		la[now] = edge[i].last;
+		if ( !vis[i>>1] ){
+			vis[i>>1] = true;
+			dfs ( edge[i].v );
+			ans[++ka] = ( ( i & 1 ) ? -1 : 1 ) * id[i];
+		}
+	}
+}
+
+void euler ()
+{
+	lp ( i, 0, ke ) deg[edge[i].u] ^= 1;
+	
+	lpi ( i, 1, n ){
+		if ( deg[i] ){
+			printf ( "NO\n" );
+			return;
+		}
+	}
+	
+	lpi ( i, 1, n ){
+		if ( ~la[i] ){
+			dfs ( i );
+			break;
+		}
+	}
+}
+
+int main (int argc, char **argv)
+{
+	if (argc != 3){
+		cerr << "3 arguments required!" << endl;
+		return 1;
+	}
+
+	ifstream INFO (argv[1]);
+	INFO >> n >> m;
+	vector < pair < int, int > > point ( n + 1 );
+	for ( int i = 1; i <= n; ++i ){
+		int x, y;
+		INFO >> x >> y;
+		point[i] = make_pair ( x, y );
+	}
+	vector < pair < int, int > > es ( m + 1 );
+	for ( int i = 1; i <= m; ++i ){
+		int u, v;
+		INFO >> u >> v;
+		es[i] = make_pair ( u, v );
+	}
+	INFO.close ();
+
+	ifstream RES (argv[2]);
+	int S;
+	RES >> S;
+	double t1, t2, t3;
+	RES >> t1 >> t2 >> t3;
+	vector < int > count ( m + 1 );
+	for ( int i = 1; i <= m; ++i ){
+//		cerr << count[i] << " ";
+		RES >> count[i];
+	}
+	cerr << endl;
+	RES.close ();
+
+	int cm = 0;
+	INIT ( la, -1 );
+	for ( int i = 1; i <= m; ++i ){
+		cm += count[i];
+		for ( int c = 0; c < count[i]; ++c ){
+			edge[ke].setd ( es[i].first, es[i].second, la[es[i].first] );
+			id[ke] = i;
+			la[es[i].first] = ke++;
+			edge[ke].setd ( es[i].second, es[i].first, la[es[i].second] );
+			id[ke] = i;
+			la[es[i].second] = ke++;
+		}
+	}
+
+//	for ( int i = 0; i < ke; i += 2 ) cerr << edge[i].u << " " << edge[i].v << endl;
+
+	euler ();
+
+	if (ka == cm){
+		vector < int > vec;
+		for ( int i = cm; i > 0; --i ){
+			if ( ans[i] > 0 ) vec.push_back (es[ans[i]].first);
+			else vec.push_back (es[-ans[i]].second);
+		}
+		for ( int i = cm; i > 1; --i) {
+			if ( ans[i] > 0 ) assert ( es[ans[i]].second == vec[cm - i + 1] );
+			else assert ( es[-ans[i]].first == vec[cm - i + 1] );
+		}
+//		for ( auto x : vec ) cerr << x << " "; cerr << endl;
+		int beg = -1;
+		for ( int i = 0; i < cm; ++i ){
+			if ( vec[i] == S ){
+				beg = i;
+				break;
+			}
+		}
+		assert ( beg != -1 );
+		for ( int i = 0; i < cm; ++i ){
+			if ( !i || getdist ( point[vec[(beg+i-1)%cm]], point[vec[(beg+i)%cm]] ) > 2 ){
+				cout << point[vec[(beg+i)%cm]].first << " " << point[vec[(beg+i)%cm]].second << endl;
+			}
+		}
+	}else{
+		cerr << "ERROR, no euler cycle exists!" << endl;
+	}
+}
\ No newline at end of file

graph.h

+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <utility>
+#include <algorithm>
+
+template < typename T >
+struct Edge
+{
+	Edge ( int u, int v, int id, T w );
+
+	int u, v, id;
+	T w;
+};
+
+struct Floyd;
+
+template < typename T >
+class Graph
+{
+public:
+	Graph ( int n );
+
+	void addedge ( int u, int v, T w );
+	std::vector < Edge < T > > & getedges ( int u );
+
+	int get_n () const;
+	int get_m () const;
+
+	Edge < T > get_edge_with_id ( int id ) const;
+
+private:
+	friend struct Floyd;
+
+	int n, m;
+	std::vector < std::vector < Edge < T > > > edge;
+	std::vector < Edge < T > > alle;
+};
+
+struct Floyd
+{
+	Floyd ( Graph < double > g )
+		: original_g ( g ), graph ( g.get_n () ), dist ( g.get_n () + 1, std::vector < double > ( g.get_n () + 1, INFINITY ) ), pre ( g.get_n () + 1, std::vector < int > ( g.get_n () + 1 ) )
+	{
+		int n = g.get_n ();
+		for ( int i = 1; i <= n; ++i ) {
+			dist[i][i] = 0;
+			pre[i][i] = -1;
+			for ( auto v : g.edge[i] ) {
+				dist[i][v.v] = v.w;
+				pre[i][v.v] = v.id;
+			}
+		}
+		for ( int k = 1; k <= n; ++k ) {
+			for ( int i = 1; i <= n; ++i ) {
+				for ( int j = 1; j <= n; ++j ) {
+					if ( dist[i][k] + dist[k][j] < dist[i][j] ) {
+						dist[i][j] = dist[i][k] + dist[k][j];
+						pre[i][j] = pre[k][j];
+					}
+				}
+			}
+		}
+		for ( int i = 1; i <= n; ++i ) {
+			for ( int j = i + 1; j <= n; ++j ) {
+				graph.addedge ( i, j, dist[i][j] );
+			}
+		}
+	}
+
+	std::vector < int > get_path ( int u, int v )
+	{
+		std::vector < int > ans;
+		for ( int t = v; t != u; t = ( ( original_g.get_edge_with_id ( pre[u][t] ).u == t ) ? original_g.get_edge_with_id ( pre[u][t] ).v : original_g.get_edge_with_id ( pre[u][t] ).u ) ) {
+			ans.push_back ( pre[u][t] );
+		}
+		reverse ( ans.begin (), ans.end () );
+		return ans;
+	}
+
+	std::vector < std::vector < double > > dist;
+	std::vector < std::vector < int > > pre;
+	Graph < double > original_g;
+	Graph < double > graph;
+};
+
+
+
+template < typename T >
+Graph<T>::Graph ( int n )
+	: n ( n ), m ( 0 ), edge ( n + 2 )
+{
+}
+
+template<typename T>
+void Graph<T>::addedge ( int u, int v, T w )
+{
+	++m;
+	edge[u].push_back ( Edge < T > ( u, v, m, w ) );
+	edge[v].push_back ( Edge < T > ( v, u, m, w ) );
+	alle.push_back ( Edge < T > ( u, v, m, w ) );
+}
+
+template<typename T>
+std::vector<Edge<T>>& Graph<T>::getedges ( int u )
+{
+	return edge[u];
+}
+
+template<typename T>
+int Graph<T>::get_n () const
+{
+	return n;
+}
+
+template<typename T>
+int Graph<T>::get_m () const
+{
+	return m;
+}
+
+template<typename T>
+Edge<T> Graph<T>::get_edge_with_id ( int id ) const
+{
+	return alle[id - 1];
+}
+
+struct Weight
+{
+	Weight ( double length = 0, double weight = 0 )
+		: length ( length ), weight ( weight )
+	{}
+
+	double length;
+	double weight;
+};
+
+typedef Graph < Weight > WGraph;
+
+template<typename T>
+Edge<T>::Edge ( int u, int v, int id, T w )
+	: u ( u ), v ( v ), id ( id ), w ( w )
+{
+}

main.cpp

+#include <bits/stdc++.h>
+
+#include "blossom.h"
+#include "simulate_annealing.h"
+#include "consts.h"
+
+using namespace std;
+
+int main ( int argc, char **argv )
+{
+	if ( argc != 4 ) {
+		cerr << "4 arguments required!" << endl;
+		return 1;
+	}
+
+	double max_length;
+	{
+		stringstream str_buf;
+		str_buf << argv[3];
+		str_buf >> max_length;
+	}
+
+	int S;
+
+	int n, m, spc;
+	ifstream IN ( argv[1] );
+	IN >> n >> m >> spc;
+	vector < int > node_type ( n + 2 );
+	for ( int i = 1; i <= n; ++i ) {
+		IN >> node_type[i];
+	}
+	vector < int > deg ( n + 2, 0 );
+	Graph < double > g ( n );
+	double sum_w = 0;
+	for ( int i = 1; i <= m; ++i ) {
+		int u, v;
+		double w;
+		IN >> u >> v >> w;
+		g.addedge ( u, v, w );
+		++deg[u], ++deg[v];
+		sum_w += w;
+	}
+	vector < double > node_weight ( n + 2 );
+	for ( int i = 1; i <= spc; ++i ) {
+		string sp_type;
+		int node;
+		IN >> sp_type >> node;
+		if ( sp_type[0] == 'P' ) node_weight[node] += WEIGHT_PARK;
+		else if ( sp_type[0] == 'A' ) node_weight[node] += WEIGHT_PICNIC;
+		else if ( sp_type[0] == 'R' ) node_weight[node] += WEIGHT_RESTROOM;
+		else if ( sp_type[0] == 'C' ) node_weight[node] += WEIGHT_CAMP;
+		else if ( sp_type[0] == 'G' ) S = node;
+	}
+	IN.close ();
+
+	WGraph wg ( n );
+	double all_weight_sum = 0;
+	for ( int i = 1; i <= n; ++i ) {
+		all_weight_sum += node_weight[i];
+	}
+	double c[3] = { 0, 0, 0 };
+	for ( int i = 1; i <= m; ++i ) {
+		auto ge = g.get_edge_with_id ( i );
+		int tu = node_type[ge.u], tv = node_type[ge.v];
+		double coef;
+		if ( tu == 1 && tv == 1 ) {
+			coef = WEIGHT_TRAIL;
+			c[0] += ge.w;
+		} else if ( tu == 2 && tv == 2 ) {
+			coef = WEIGHT_ROAD;
+			c[1] += ge.w;
+		} else {
+			coef = ( WEIGHT_TRAIL + WEIGHT_ROAD ) / 2;
+			c[2] += ge.w;
+		}
+		wg.addedge ( ge.u, ge.v, Weight ( ge.w, coef ) );
+		all_weight_sum += ge.w * coef;
+	}
+//	cerr << c[0] << " " << c[1] << " " << c[2] << endl;
+
+	Floyd floyd ( g );
+
+	vector < int > all_set;
+	for ( int i = 1; i <= n; ++i ) {
+		if ( deg[i] & 1 ) {
+			all_set.push_back ( i );
+		}
+	}
+	cout << "Length of Chinese Postman Problem: " << solve_cpp ( floyd.graph, all_set ).first + sum_w << "/" << sum_w << endl;
+
+	cout << endl << "Running simulated annealing with parameters:" << endl;
+	cout << "Node count: " << n << endl;
+	cout << "Edge count: " << m << endl;
+	cout << "Starting point: " << S << endl;
+	cout << "Maximum length: " << max_length << endl;
+	cout << "Sum of weight of all nodes and edges: " << all_weight_sum << endl;
+
+	auto sa_result = solve_sa ( wg, floyd, node_weight, S, max_length );
+
+	cout << "Result of simulated annealing: " << sa_result.weight_sum << "/" << all_weight_sum << ", path length: " << sa_result.path_length << endl;
+
+	ofstream OUT ( argv[2] );
+	OUT << S << endl;
+	OUT << sa_result.weight_sum << " " << all_weight_sum << " " << sa_result.path_length << endl;
+	for ( int i = 1; i <= m; ++i ) {
+		OUT << sa_result.edge_count[i] << " ";
+	}
+	OUT << endl;
+	OUT.close ();
+}
\ No newline at end of file

png_to_graph.py

+#!/usr/bin/env python3
+
+import sys
+import os
+import png
+import math
+import argparse
+
+
+
+PIXEL_LENGTH = 0
+
+BLACK = (0, 0, 0)
+RED = (255, 0, 0)
+WHITE = (255, 255, 255)
+
+
+
+def main ():
+	parser = argparse.ArgumentParser (prog='png_to_graph.py', description='Convert PNG to graph')
+	parser.add_argument ('-i', '--input', dest='input', metavar='input_file', help='Input file name', required=True)
+	parser.add_argument ('-o', '--output', dest='output', metavar='output_file', help='Output file name', required=True)
+	parser.add_argument ('-d', '--desc', dest='desc', metavar='desc_file', help='Description file name', required=True)
+	parser.add_argument ('-s', '--shrink', dest='shrink', help='Shrink the graph', action='store_true', default=False)
+	parser.add_argument ('-l', '--length', dest='length', metavar='length', help='Number of pixels of a mile', type=int, required=True)
+	parser.add_argument ('-v', '--vertex', dest='vertex', metavar='vertex', help='Output the vertex information to this file', default=None)
+
+	args = parser.parse_args ()
+
+	PIXEL_LENGTH = 1 / args.length
+
+	png_reader = png.Reader (file=open (args.input, 'rb'))
+	png_info = png_reader.asRGB ()
+	
+	W, H = png_info[0], png_info[1]
+
+	pixels_raw = list (png_info[2])
+
+	pixels = [[WHITE] * (W + 2)]
+	for i in range (H):
+		pixels.append ([WHITE])
+		for j in range (W):
+			pixels[i+1].append ((pixels_raw[i][3*j], pixels_raw[i][3*j+1], pixels_raw[i][3*j+2] ))
+		pixels[i+1].append (WHITE)
+	pixels.append ([WHITE] * (W + 2))
+
+	vertex = {}
+	color = [0]
+	edge = []
+
+	def getcolor (p):
+		if pixels[p[0]][p[1]] == BLACK:
+			return 1
+		else:
+			return 2
+
+	def getdist (u, v):
+		def sqr(x):
+			return x*x
+		return math.sqrt (sqr(u[0]-v[0]) + sqr(u[1]-v[1]))
+
+	def addedge (u, v):
+		def addnode (u):
+			if u not in vertex:
+				vertex[u] = len(vertex) + 1
+				color.append (getcolor (u))
+		addnode (u)
+		addnode (v)
+		if vertex[u] < vertex[v]:
+			edge.append ((vertex[u], vertex[v], getdist (u, v)))
+
+	for i in range (1, H+1):
+		for j in range (1, W+1):
+			if pixels[i][j] != WHITE:
+				for dx in range (-1, 2):
+					for dy in range (-1, 2):
+						if (dx or dy) and (pixels[i+dx][j+dy] != WHITE):
+							addedge ((i, j), (i+dx, j+dy))
+
+	sp_points = []
+	spp_set = set ()
+	with open (args.desc) as DESC:
+		for line in DESC.readlines ():
+			content = line.strip ().split ()
+			if len(content) == 0:
+				continue
+			tp = content[0]
+			y = int (content[1])
+			x = int (content[2])
+			min_dist = None
+			min_p = None
+			for i in vertex.keys ():
+				cur_dist = getdist (i, (x, y))
+				if min_dist is None or cur_dist < min_dist:
+					min_dist = cur_dist
+					min_p = vertex[i]
+			sp_points.append ((tp, min_p))
+			spp_set.add (min_p)
+
+	if args.shrink:
+		n = len(vertex)
+		g = [{} for i in range (n+1)]
+
+		def insedge (u, v, w):
+			g[u][v] = w
+			g[v][u] = w
+
+		def deledge (u, v):
+			g[u].pop (v)
+			g[v].pop (u)
+
+		for (u, v, w) in edge:
+			insedge (u, v, w)
+		
+		rem_vertex = {}
+		for i in range (1, n+1):
+			rem_vertex[i] = 0
+
+		for i in range (1, n+1):
+			if len(g[i]) == 2 and i not in spp_set:
+				gl = list(g[i].keys ())
+				gw = sum (g[i].values ())
+				if color[i] == color[gl[0]] and color[i] == color[gl[1]]:
+					deledge (i, gl[0])
+					deledge (i, gl[1])
+					rem_vertex.pop (i)
+					insedge (gl[0], gl[1], gw)
+
+		rem_count = 0
+		for i in rem_vertex:
+			rem_count += 1
+			rem_vertex[i] = rem_count
+		new_vertex = {}
+		new_color = [0] * (rem_count + 1)
+		new_edge = []
+		for i in vertex:
+			if vertex[i] in rem_vertex:
+				new_vertex[i] = rem_vertex[vertex[i]]
+				new_color[new_vertex[i]] = color[vertex[i]]
+		for i in rem_vertex:
+			ci = rem_vertex[i]
+			for j in g[i]:
+				cj = rem_vertex[j]
+				if ci < cj:
+					new_edge.append ((ci, cj, g[i][j]))
+		vertex = new_vertex
+		edge = new_edge
+		color = new_color
+		for i in range (len(sp_points)):
+			sp_points[i] = (sp_points[i][0], rem_vertex[sp_points[i][1]])
+		new_spp_set = set ()
+		for i in spp_set:
+			new_spp_set.add (rem_vertex[i])
+		spp_set = new_spp_set
+
+	print ('Total length of edges:', sum (map (lambda x: x[2], edge)) * PIXEL_LENGTH)
+
+	with open (args.output, 'w') as OUT:
+		print (len(vertex), len(edge), len(sp_points), file=OUT)
+		print (' '.join (map(str, color[1:])), file=OUT)
+		for i in edge:
+			print (i[0], i[1], i[2] * PIXEL_LENGTH, file=OUT)
+		for i in sp_points:
+			print (i[0], i[1], file=OUT)
+
+	if args.vertex is not None:
+		rev_ver = [None] * (rem_count + 1)
+		for i in vertex:
+			rev_ver[vertex[i]] = i
+		with open (args.vertex, 'w') as VERT:
+			print (len(vertex), len(edge), file=VERT)
+			for i in rev_ver[1:]:
+				print (i[0], i[1], file=VERT)
+			for i in edge:
+				print (i[0], i[1], file=VERT)
+	
+
+
+if __name__ == '__main__':
+	main ()

simulate_annealing.cpp

+#include <bits/stdc++.h>
+
+#include "blossom.h"
+#include "simulate_annealing.h"
+
+#include <random>
+#include <chrono>
+
+using namespace std;
+
+static const double DIST_DROP_RATE = 0.999;
+static const double P_DROP_RATE = 0.9998;
+static const double DIST_START = 5;
+static const double P_START = 1;
+static const int ROUND_COUNT = 50000;
+
+class SA_Solver
+{
+public:
+	struct State
+	{
+		State ( SA_Solver *solver = nullptr )
+			: inv ( ( solver == nullptr ? 0 : solver->n + 2 ), 0 ), ine ( ( solver == nullptr ? 0 : solver->m + 2 ), false ), solver ( solver )
+		{}
+
+		void addedge ( int ei )
+		{
+			ine[ei] = true;
+			auto edge = solver->g.get_edge_with_id ( ei );
+			++inv[edge.u];
+			++inv[edge.v];
+		}
+
+		void deledge ( int ei )
+		{
+			ine[ei] = false;
+			auto edge = solver->g.get_edge_with_id ( ei );
+			--inv[edge.u];
+			--inv[edge.v];
+		}
+
+		State get_next ( double p ) const
+		{
+			State res = *this;
+			int ec = 0;
+			for ( int i = 1; i <= solver->m; ++i ) {
+				if ( res.ine[i] ) {
+					++ec;
+				}
+			}
+			for ( int i = 1; i <= solver->m; ++i ) {
+				if ( res.ine[i] && solver->urd ( solver->random_engine ) < p * solver->m / ec ) {
+//					cerr << "DEL : " << i << " " << "(" << solver->g.get_edge_with_id ( i ).u << "," << solver->g.get_edge_with_id ( i ).v << ") " << endl;
+					res.deledge ( i );
+				}
+			}
+			int nc = 0;
+			for ( int i = 1; i <= solver->n; ++i ) {
+				if ( res.inv[i] ) {
+					for ( auto edge : solver->g.getedges ( i ) ) {
+						if ( !res.ine[edge.id] ) {
+							++nc;
+						}
+					}
+				}
+			}
+			for ( int i = 1; i <= solver->n; ++i ) {
+				if ( res.inv[i] ) {
+					for ( auto edge : solver->g.getedges ( i ) ) {
+						if ( !res.ine[edge.id] && solver->urd ( solver->random_engine ) < p * solver->m / nc ) {
+//							cerr << "ADD : " << edge.id << " " << "(" << solver->g.get_edge_with_id ( edge.id ).u << "," << solver->g.get_edge_with_id ( edge.id ).v << ") " << endl;
+							res.addedge ( edge.id );
+						}
+					}
+				}
+			}
+			//for ( int i = 1; i <= solver->n; ++i ) {
+			//	if ( res.inv[i] ) cerr << i << " ";
+			//}
+			//cerr << endl;
+			//for ( int i = 1; i <= solver->m; ++i ) {
+			//	if ( res.ine[i] ) cerr << "(" << solver->g.get_edge_with_id ( i ).u << "," << solver->g.get_edge_with_id ( i ).v << ") ";
+			//}
+			//cerr << endl;
+			//cerr << ec << " " << nc << endl;
+			return res;
+		}
+
+		pair < pair < double, double >, vector < int > > get_length () const
+		{
+			//for ( int i = 1; i <= solver->n; ++i ) {
+			//	if ( inv[i] ) cerr << i << " ";
+			//}
+			//cerr << endl;
+			//for ( int i = 1; i <= solver->m; ++i ) {
+			//	if ( ine[i] ) cerr << "(" << solver->g.get_edge_with_id ( i ).u << "," << solver->g.get_edge_with_id ( i ).v << ") ";
+			//}
+			//cerr << endl;
+			vector < bool > visit ( solver->n + 2 );
+			for ( int i = 1; i <= solver->n; ++i ) {
+				if ( inv[i] ) {
+					dfs ( i, visit );
+					break;
+				}
+			}
+			for ( int i = 1; i <= solver->n; ++i ) {
+				if ( inv[i] && !visit[i] ) {
+					return make_pair ( make_pair ( -INFINITY, -INFINITY ), vector < int > () );
+				}
+			}
+
+			double sum_e = 0, sum_w = 0;
+			for ( int i = 1; i <= solver->n; ++i ) {
+				if ( inv[i] ) {
+					sum_w += solver->node_weight[i];
+				}
+			}
+			//for ( int i = 1; i <= solver->m; ++i ) {
+			//	if ( ine[i] ) {
+			//		sum_e += solver->g.get_edge_with_id ( i ).w.length;
+			//	}
+			//}
+			vector < int > odd;
+			for ( int i = 1; i <= solver->n; ++i ) {
+				if ( inv[i] & 1 ) {
+					odd.push_back ( i );
+				}
+			}
+			auto result = solve_cpp ( solver->md.graph, odd );
+			sum_e += result.first;
+			vector < int > es ( solver->m + 1 );
+			for ( int i = 1; i <= solver->m; ++i ) {
+				if ( ine[i] ) {
+					if ( !es[i] ) {
+						sum_e += solver->g.get_edge_with_id ( i ).w.length;
+						sum_w += solver->g.get_edge_with_id ( i ).w.weight * solver->g.get_edge_with_id ( i ).w.length;
+					}
+					++es[i];
+				}
+			}
+			for ( auto x : odd ) {
+				int y = result.second[x];
+				if ( x < y ) {
+					vector < int > ces = solver->md.get_path ( x, y );
+					for ( auto e : ces ) {
+						if ( !es[e] ) {
+							sum_e += solver->g.get_edge_with_id ( e ).w.length;
+							sum_w += solver->g.get_edge_with_id ( e ).w.weight * solver->g.get_edge_with_id ( e ).w.length;
+						}
+						++es[e];
+					}
+				}
+			}
+			return make_pair ( make_pair ( sum_e, sum_w ), es );
+		}
+
+		double get_value ( double max_length ) const
+		{
+			pair < double, double > res = get_length ().first;
+			if ( res.first > max_length ) return max_length - res.first;
+			else return res.second;
+		}
+
+		void dfs ( int now, vector < bool > &visit ) const
+		{
+			if ( !inv[now] || visit[now] ) {
+				return;
+			}
+			visit[now] = true;
+			for ( auto edge : solver->g.getedges ( now ) ) {
+				if ( ine[edge.id] ) {
+					dfs ( edge.v, visit );
+				}
+			}
+		}
+
+		vector < int > inv;
+		vector < bool > ine;
+		SA_Solver *solver;
+	};
+
+public:
+	SA_Solver ( WGraph g, Floyd md, vector < double > node_weight )
+		: g ( g ), md ( md ), node_weight ( node_weight ), n ( g.get_n () ), m ( g.get_m () )
+		, random_engine ( static_cast < unsigned > ( chrono::system_clock::now ().time_since_epoch ().count () ) )
+		, urd ( 0.0, 1.0 )
+	{}
+
+	State solve ( int S, double max_length )
+	{
+		double disturb = DIST_START / m, p = P_START;
+		State cur = State ( this ), now;
+		cur.inv[S] = 2;
+		State mp = cur;
+		double cv = cur.get_value ( max_length ), maxv = cv, nv;
+		int output_freq = 2000 / n + 1;
+		for ( int rd = 1; rd <= ROUND_COUNT; ++rd ) {
+			now = cur.get_next ( disturb );
+			nv = now.get_value ( max_length );
+			if ( rd % output_freq == 0 ) {
+				cerr << "\r";
+				for ( int i = 1; i <= 100; ++i ) {
+					cerr << " ";
+				}
+				cerr << "\r";
+				cerr << "round " << rd << "/" << ROUND_COUNT << " : " << setw ( 10 ) << fixed << nv << " " << setw ( 10 ) << fixed << cv << " " << setw ( 10 ) << fixed << maxv << " : with param " << setw ( 10 ) << fixed << disturb << " " << setw ( 10 ) << fixed << p;
+			}
+			if ( nv > cv || urd ( random_engine ) < exp ( ( nv - cv ) / p ) ) {
+				cur = now;
+				cv = nv;
+			}
+			if ( cv > maxv ) {
+				mp = cur;
+				maxv = cv;
+			}
+			disturb = max ( disturb * DIST_DROP_RATE, 2.0 / m );
+			p *= P_DROP_RATE;
+		}
+		cerr << "\r";
+		for ( int i = 1; i <= 100; ++i ) {
+			cerr << " ";
+		}
+		cerr << "\r";
+		return mp;
+	}
+
+private:
+	friend struct State;
+
+	int n, m;
+	WGraph g;
+	Floyd md;
+	vector < double > node_weight;
+	mt19937 random_engine;
+	uniform_real_distribution < double > urd;
+};
+
+sa_result_t solve_sa ( WGraph g, Floyd md, vector<double> node_weight, int S, double max_length )
+{
+	SA_Solver solver ( g, md, node_weight );
+	SA_Solver::State result = solver.solve ( S, max_length );
+	auto eval = result.get_length ();
+	if ( eval.first.first > max_length ) return sa_result_t ( 0, 0, vector < int > ( g.get_n () + 2 ) );
+	return sa_result_t ( eval.first.first, eval.first.second, eval.second );
+}
+
+sa_result_t::sa_result_t ( double path_length, double weight_sum, std::vector<int> edge_count )
+	: path_length ( path_length ), weight_sum ( weight_sum ), edge_count ( edge_count)
+{
+}

simulate_annealing.h

+#pragma once
+
+#include <vector>
+#include "graph.h"
+
+struct sa_result_t
+{
+	sa_result_t ( double path_length, double weight_sum, std::vector < int > edge_count );
+
+	double path_length;
+	double weight_sum;
+	std::vector < int > edge_count;
+};
+
+sa_result_t solve_sa ( WGraph g, Floyd md, std::vector < double > node_weight, int S, double max_length );