text:
// Copyright (c) 2011 Niels Martin Hansen <nielsm@aegisub.org>
//
// Permission to use, copy, modify, and distribute this software for any
// purpose with or without fee is hereby granted, provided that the above
// copyright notice and this permission notice appear in all copies.
//
// THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
// WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
// MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
// ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
// WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
// ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
// OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
//
// Aegisub Project http://www.aegisub.org/

/// @file subtitle_format_ebu3264.cpp
/// @see subtitle_format_ebu3264.h
/// @ingroup subtitle_io

// This implements support for the EBU tech 3264 (1991) subtitling data exchange format.
// Work on support for this format was sponsored by Bandai.

#include "config.h"

#include "subtitle_format_ebu3264.h"

#include "ass_dialogue.h"
#include "ass_file.h"
#include "ass_style.h"
#include "compat.h"
#include "dialog_export_ebu3264.h"
#include "options.h"
#include "text_file_writer.h"

#include <libaegisub/charset_conv.h>
#include <libaegisub/exception.h>
#include <libaegisub/io.h>
#include <libaegisub/line_wrap.h>
#include <libaegisub/of_type_adaptor.h>

#include <boost/algorithm/string/replace.hpp>
#include <fstream>

namespace
{
#pragma pack(push, 1)
	/// General Subtitle Information block as it appears in the file
	struct BlockGSI
	{
		char cpn[3];    ///< code page number
		char dfc[8];    ///< disk format code
		char dsc;       ///< display standard code
		char cct[2];    ///< character code table number
		char lc[2];     ///< language code
		char opt[32];   ///< original programme title
		char oet[32];   ///< original episode title
		char tpt[32];   ///< translated programme title
		char tet[32];   ///< translated episode title
		char tn[32];    ///< translator name
		char tcd[32];   ///< translator contact details
		char slr[16];   ///< subtitle list reference code
		char cd[6];     ///< creation date
		char rd[6];     ///< revision date
		char rn[2];     ///< revision number
		char tnb[5];    ///< total number of TTI blocks
		char tns[5];    ///< total number of subtitles
		char tng[3];    ///< total number of subtitle groups
		char mnc[2];    ///< maximum number of displayable characters in a row
		char mnr[2];    ///< maximum number of displayable rows
		char tcs;       ///< time code: status
		char tcp[8];    ///< time code: start of programme
		char tcf[8];    ///< time code: first in-cue
		char tnd;       ///< total number of disks
		char dsn;       ///< disk sequence number
		char co[3];     ///< country of origin
		char pub[32];   ///< publisher
		char en[32];    ///< editor's name
		char ecd[32];   ///< editor's contact details
		char unused[75];
		char uda[576];  ///< user defined area
	};

	/// Text and Timing Information block as it appears in the file
	struct BlockTTI
	{
		uint8_t     sgn; ///< subtitle group number
		uint16_t    sn;  ///< subtitle number
		uint8_t     ebn; ///< extension block number
		uint8_t     cs;  ///< cumulative status
		EbuTimecode tci; ///< time code in
		EbuTimecode tco; ///< time code out
		uint8_t     vp;  ///< vertical position
		uint8_t     jc;  ///< justification code
		uint8_t     cf;  ///< comment flag
		char    tf[112]; ///< text field
	};
#pragma pack(pop)

	/// A block of text with basic formatting information
	struct EbuFormattedText
	{
		std::string text; ///< Text in this block
		bool underline;   ///< Is this block underlined?
		bool italic;      ///< Is this block italic?
		bool word_start;  ///< Is it safe to line-wrap between this block and the previous one?
		EbuFormattedText(std::string t, bool u = false, bool i = false, bool ws = true) : text(std::move(t)), underline(u), italic(i), word_start(ws) { }
	};
	typedef std::vector<EbuFormattedText> EbuTextRow;

	/// Formatting character constants
	const unsigned char EBU_FORMAT_ITALIC[]     = "\x81\x80";
	const unsigned char EBU_FORMAT_UNDERLINE[]  = "\x83\x82";
	const unsigned char EBU_FORMAT_BOXING[]     = "\x85\x84";
	const unsigned char EBU_FORMAT_LINEBREAK    = '\x8a';
	const unsigned char EBU_FORMAT_UNUSED_SPACE = '\x8f';

	/// intermediate format
	class EbuSubtitle
	{
		void ProcessOverrides(AssDialogueBlockOverride *ob, bool &underline, bool &italic, int &align, bool style_underline, bool style_italic)
		{
			for (auto const& t : ob->Tags)
			{
				if (t.Name == "\\u")
					underline = t.Params[0].Get<bool>(style_underline);
				else if (t.Name == "\\i")
					italic = t.Params[0].Get<bool>(style_italic);
				else if (t.Name == "\\an")
					align = t.Params[0].Get<int>(align);
				else if (t.Name == "\\a" && !t.Params[0].omitted)
					align = AssStyle::SsaToAss(t.Params[0].Get<int>());
			}
		}

		void SetAlignment(int ass_alignment)
		{
			if (ass_alignment < 1 || ass_alignment > 9)
				ass_alignment = 2;

			vertical_position = static_cast<VerticalPosition>(ass_alignment / 3);
			justification_code = static_cast<JustificationCode>((ass_alignment - 1) % 3 + 1);
		}

	public:
		enum CumulativeStatus
		{
			NotCumulative    = 0,
			CumulativeStart  = 1,
			CulumativeMiddle = 2,
			CumulativeEnd    = 3
		};

		enum JustificationCode
		{
			UnchangedPresentation = 0,
			JustifyLeft           = 1,
			JustifyCentre         = 2,
			JustifyRight          = 3
		};

		// note: not set to constants from spec
		enum VerticalPosition
		{
			PositionTop    = 2,
			PositionMiddle = 1,
			PositionBottom = 0
		};

		int group_number; ///< always 0 for compat
		/// subtitle number is assigned when generating blocks
		CumulativeStatus cumulative_status; ///< always NotCumulative for compat
		int time_in;       ///< frame number
		int time_out;      ///< frame number
		bool comment_flag; ///< always false for compat
		JustificationCode justification_code; ///< never Unchanged presentation for compat
		VerticalPosition vertical_position;   ///< translated to row on tti conversion
		std::vector<EbuTextRow> text_rows;    ///< text split into rows, still unicode

		EbuSubtitle()
		: group_number(0)
		, cumulative_status(NotCumulative)
		, time_in(0)
		, time_out(0)
		, comment_flag(false)
		, justification_code(JustifyCentre)
		, vertical_position(PositionBottom)
		, text_rows()
		{
		}

		void SplitLines(int max_width, int split_type)
		{
			// split_type is an SSA wrap style number
			if (split_type == 2) return; // no wrapping here!
			if (split_type < 0) return;
			if (split_type > 4) return;

			std::vector<EbuTextRow> new_text;
			new_text.reserve(text_rows.size());

			for (auto const& row : text_rows)
			{
				// Get lengths of each word
				std::vector<size_t> word_lengths;
				for (auto const& cur_block : row)
				{
					if (cur_block.word_start)
						word_lengths.push_back(0);
					word_lengths.back() += cur_block.text.size();
				}

				std::vector<size_t> split_points = agi::get_wrap_points(word_lengths, (size_t)max_width, (agi::WrapMode)split_type);

				if (split_points.empty())
				{
					// Line doesn't need splitting, so copy straight over
					new_text.push_back(row);
					continue;
				}

				// Apply the splits
				new_text.emplace_back();
				size_t cur_word = 0;
				size_t split_point = 0;
				for (auto const& cur_block : row)
				{
					if (cur_block.word_start && split_point < split_points.size())
					{
						if (split_points[split_point] == cur_word)
						{
							new_text.emplace_back();
							++split_point;
						}
						++cur_word;
					}

					new_text.back().push_back(cur_block);
				}
			}

			// replace old text
			swap(text_rows, new_text);
		}

		bool CheckLineLengths(int max_width) const
		{
			for (auto const& row : text_rows)
			{
				int line_length = 0;
				for (auto const& block : row)
					line_length += block.text.size();

				if (line_length > max_width)
					// early return as soon as any line is over length
					return false;
			}
			// no lines failed
			return true;
		}

		void SetTextFromAss(AssDialogue *line, bool style_underline, bool style_italic, int align, int wrap_mode)
		{
			boost::ptr_vector<AssDialogueBlock> blocks(line->ParseTags());

			text_rows.clear();
			text_rows.emplace_back();

			// current row being worked on
			EbuTextRow *cur_row = &text_rows.back();

			// create initial text part
			cur_row->emplace_back("", style_underline, style_italic, true);

			bool underline = style_underline, italic = style_italic;

			for (auto& b : blocks)
			{
				switch (b.GetType())
				{
					case AssBlockType::PLAIN:
					// find special characters and convert them
					{
						std::string text = b.GetText();

						boost::replace_all(text, "\\t", " ");

						size_t start = 0;
						for (size_t i = 0; i < text.size(); ++i)
						{
							if (text[i] != ' ' && (i + 1 >= text.size() || text[i] != '\\' || (text[i + 1] != 'N' && text[i + 1] != 'n')))
								continue;

							// add first part of text to current part
							cur_row->back().text.append(begin(text) + start, begin(text) + i);

							// process special character
							if (text[i] == '\\' && (text[i + 1] == 'N' || wrap_mode == 1))
							{
								// create a new row with current style
								text_rows.emplace_back();
								cur_row = &text_rows.back();
								cur_row->emplace_back("", underline, italic, true);
							}
							else // if (substr == " " || substr == "\\n")
							{
								cur_row->back().text.append(" ");
								cur_row->emplace_back("", underline, italic, true);
							}

							start = i + (text[i] == '\\');
						}

						// add the remaining text
						cur_row->back().text.append(begin(text) + start, end(text));

						// convert \h to regular spaces
						// done after parsing so that words aren't split on \h
						boost::replace_all(cur_row->back().text, "\\h", " ");
					}
					break;

					case AssBlockType::OVERRIDE:
					// find relevant tags and process them
					{
						AssDialogueBlockOverride *ob = static_cast<AssDialogueBlockOverride*>(&b);
						ob->ParseTags();
						ProcessOverrides(ob, underline, italic, align, style_underline, style_italic);

						// apply any changes
						if (underline != cur_row->back().underline || italic != cur_row->back().italic)
						{
							if (cur_row->back().text.empty())
							{
								// current part is empty, we can safely change formatting on it
								cur_row->back().underline = underline;
								cur_row->back().italic = italic;
							}
							else
							{
								// create a new empty part with new style
								cur_row->emplace_back("", underline, italic, false);
							}
						}
					}
					break;

				default:
					// ignore block, we don't want to output it (drawing or comment)
					break;
				}
			}

			SetAlignment(align);
		}
	};

	std::vector<EbuSubtitle> convert_subtitles(AssFile &copy, EbuExportSettings const& export_settings)
	{
		SubtitleFormat::StripComments(copy);
		copy.Sort();
		SubtitleFormat::RecombineOverlaps(copy);
		SubtitleFormat::MergeIdentical(copy);

		int line_wrap_type = copy.GetScriptInfoAsInt("WrapStyle");

		agi::vfr::Framerate fps = export_settings.GetFramerate();
		EbuTimecode tcofs = export_settings.timecode_offset;
		int timecode_bias = fps.FrameAtSmpte(tcofs.h, tcofs.m, tcofs.s, tcofs.s);

		AssStyle default_style;
		std::vector<EbuSubtitle> subs_list;
		subs_list.reserve(copy.Line.size());

		// convert to intermediate format
		for (auto line : copy.Line | agi::of_type<AssDialogue>())
		{
			// add a new subtitle and work on it
			subs_list.emplace_back();
			EbuSubtitle &imline = subs_list.back();

			// some defaults for compatibility
			imline.group_number = 0;
			imline.comment_flag = false;
			imline.cumulative_status = EbuSubtitle::NotCumulative;

			// convert times
			imline.time_in = fps.FrameAtTime(line->Start) + timecode_bias;
			imline.time_out = fps.FrameAtTime(line->End) + timecode_bias;
			if (export_settings.inclusive_end_times)
				// cheap and possibly wrong way to ensure exclusive times, subtract one frame from end time
				imline.time_out -= 1;

			// convert alignment from style
			AssStyle *style = copy.GetStyle(line->Style);
			if (!style)
				style = &default_style;

			// add text, translate formatting
			imline.SetTextFromAss(line, style->underline, style->italic, style->alignment, line_wrap_type);

			// line breaking handling
			if (export_settings.line_wrapping_mode == EbuExportSettings::AutoWrap)
				imline.SplitLines(export_settings.max_line_length, line_wrap_type);
			else if (export_settings.line_wrapping_mode == EbuExportSettings::AutoWrapBalance)
				imline.SplitLines(export_settings.max_line_length, agi::Wrap_Balanced);
			else if (!imline.CheckLineLengths(export_settings.max_line_length))
			{
				if (export_settings.line_wrapping_mode == EbuExportSettings::AbortOverLength)
					throw Ebu3264SubtitleFormat::ConversionFailed(from_wx(wxString::Format(_("Line over maximum length: %s"), line->Text.get())), nullptr);
				else // skip over-long lines
					subs_list.pop_back();
			}
		}

		// produce an empty line if there are none
		// (it still has to contain a space to not get ignored)
		if (subs_list.empty())
		{
			subs_list.emplace_back();
			subs_list.back().text_rows.emplace_back();
			subs_list.back().text_rows.back().emplace_back(" ");
		}

		return subs_list;
	}

	std::string convert_subtitle_line(EbuSubtitle const& sub, agi::charset::IconvWrapper *encoder, bool enable_formatting)
	{
		std::string fullstring;
		for (auto const& row : sub.text_rows)
		{
			if (!fullstring.empty())
				fullstring += EBU_FORMAT_LINEBREAK;

			// formatting is reset at the start of every row, so keep track per row
			bool underline = false, italic = false;
			for (auto const& block : row)
			{
				if (enable_formatting)
				{
					// insert codes for changed formatting
					if (underline != block.underline)
						fullstring += EBU_FORMAT_UNDERLINE[block.underline];
					if (italic != block.italic)
						fullstring += EBU_FORMAT_ITALIC[block.italic];

					underline = block.underline;
					italic = block.italic;
				}

				// convert text to specified encoding
				fullstring += encoder->Convert(block.text);
			}
		}
		return fullstring;
	}

	void smpte_at_frame(agi::vfr::Framerate const& fps, int frame, EbuTimecode &tc)
	{
		int h=0, m=0, s=0, f=0;
		fps.SmpteAtFrame(frame, &h, &m, &s, &f);
		tc.h = h;
		tc.m = m;
		tc.s = s;
		tc.f = f;
	}

	std::vector<BlockTTI> create_blocks(std::vector<EbuSubtitle> const& subs_list, EbuExportSettings const& export_settings)
	{
		auto encoder = export_settings.GetTextEncoder();
		auto fps = export_settings.GetFramerate();

		// Teletext captions are 1-23; Open subtitles are 0-99
		uint8_t min_row = 0;
		uint8_t max_row = 100;
		if (export_settings.display_standard != EbuExportSettings::DSC_Open) {
			min_row = 1;
			max_row = 24;
		}

		uint16_t subtitle_number = 0;

		std::vector<BlockTTI> tti;
		tti.reserve(subs_list.size());
		for (auto const& sub : subs_list)
		{
			std::string fullstring = convert_subtitle_line(sub, encoder.get(),
				export_settings.display_standard == EbuExportSettings::DSC_Open);

			// construct a base block that can be copied and filled
			BlockTTI base;
			base.sgn = sub.group_number;
			base.sn = subtitle_number++;
			base.ebn = 255;
			base.cf = sub.comment_flag;
			memset(base.tf, EBU_FORMAT_UNUSED_SPACE, sizeof(base.tf));
			smpte_at_frame(fps, sub.time_in, base.tci);
			smpte_at_frame(fps, sub.time_out, base.tco);
			base.cs = sub.cumulative_status;

			if (export_settings.translate_alignments)
			{
				// vertical position
				if (sub.vertical_position == EbuSubtitle::PositionTop)
					base.vp = min_row;
				else if (sub.vertical_position == EbuSubtitle::PositionMiddle)
					base.vp = std::min<uint8_t>(min_row, max_row / 2 - (max_row / 5 * sub.text_rows.size()));
				else //if (sub.vertical_position == EbuSubtitle::PositionBottom)
					base.vp = max_row - 1;

				base.jc = sub.justification_code;
			}
			else
			{
				base.vp = max_row - 1;
				base.jc = EbuSubtitle::JustifyCentre;
			}

			// produce blocks from string
			static const size_t block_size = sizeof(((BlockTTI*)nullptr)->tf);
			uint8_t num_blocks = 0;
			for (size_t pos = 0; pos < fullstring.size(); pos += block_size)
			{
				size_t bytes_remaining = fullstring.size() - pos;

				tti.push_back(base);
				// write an extension block number if the remaining text doesn't fit in the block
				tti.back().ebn = bytes_remaining >= block_size ? num_blocks++ : 255;

				std::copy(&fullstring[pos], &fullstring[pos + std::min(block_size, bytes_remaining)], tti.back().tf);

				// Write another block for the terminator if we exactly used up
				// the last block
				if (bytes_remaining == block_size)
					tti.push_back(base);
			}
		}

		return tti;
	}

#ifdef _MSC_VER
#define snprintf _snprintf
#endif

	BlockGSI create_header(AssFile const& copy, EbuExportSettings const& export_settings)
	{
		std::string scriptinfo_title = copy.GetScriptInfo("Title");
		std::string scriptinfo_translation = copy.GetScriptInfo("Original Translation");
		std::string scriptinfo_editing = copy.GetScriptInfo("Original Editing");

		agi::charset::IconvWrapper gsi_encoder("UTF-8", "CP850");

		BlockGSI gsi;
		memset(&gsi, 0x20, sizeof(gsi)); // fill with spaces
		memcpy(gsi.cpn, "850", 3);
		switch (export_settings.tv_standard)
		{
			case EbuExportSettings::STL23:
			case EbuExportSettings::STL24:
				memcpy(gsi.dfc, "STL24.01", 8);
				break;
			case EbuExportSettings::STL29:
			case EbuExportSettings::STL29drop:
			case EbuExportSettings::STL30:
				memcpy(gsi.dfc, "STL30.01", 8);
				break;
			case EbuExportSettings::STL25:
			default:
				memcpy(gsi.dfc, "STL25.01", 8);
				break;
		}
		gsi.dsc = '0' + (int)export_settings.display_standard;
		gsi.cct[0] = '0';
		gsi.cct[1] = '0' + (int)export_settings.text_encoding;
		if (export_settings.text_encoding == EbuExportSettings::utf8)
			memcpy(gsi.cct, "U8", 2);
		memcpy(gsi.lc, "00", 2);
		gsi_encoder.Convert(scriptinfo_title.c_str(), scriptinfo_title.size(), gsi.opt, 32);
		gsi_encoder.Convert(scriptinfo_translation.c_str(), scriptinfo_translation.size(), gsi.tn, 32);
		{
			char buf[20];
			time_t now;
			time(&now);
			tm *thetime = localtime(&now);
			strftime(buf, 20, "AGI-%y%m%d%H%M%S", thetime);
			memcpy(gsi.slr, buf, 16);
			strftime(buf, 20, "%y%m%d", thetime);
			memcpy(gsi.cd, buf, 6);
			memcpy(gsi.rd, buf, 6);
			memcpy(gsi.rn, "00", 2);
			memcpy(gsi.tng, "001", 3);
			snprintf(gsi.mnc, 2, "%02u", (unsigned int)export_settings.max_line_length);
			memcpy(gsi.mnr, "99", 2);
			gsi.tcs = '1';
			EbuTimecode tcofs = export_settings.timecode_offset;
			snprintf(gsi.tcp, 8, "%02u%02u%02u%02u", (unsigned int)tcofs.h, (unsigned int)tcofs.m, (unsigned int)tcofs.s, (unsigned int)tcofs.s);
		}
		gsi.tnd = '1';
		gsi.dsn = '1';
		memcpy(gsi.co, "NTZ", 3); // neutral zone!
		gsi_encoder.Convert(scriptinfo_editing.c_str(), scriptinfo_editing.size(), gsi.en, 32);
		if (export_settings.text_encoding == EbuExportSettings::utf8)
			strncpy(gsi.uda, "This file was exported by Aegisub using non-standard UTF-8 encoding for the subtitle blocks. The TTI.TF field contains UTF-8-encoded text interspersed with the standard formatting codes, which are not encoded. GSI.CCT is set to 'U8' to signify this.", sizeof(gsi.uda));

		return gsi;
	}

	EbuExportSettings get_export_config(wxWindow *parent)
	{
		EbuExportSettings s("Subtitle Format/EBU STL");

		// Disable the busy cursor set by the exporter while the dialog is visible
		wxEndBusyCursor();
		int res = EbuExportConfigurationDialog(parent, s).ShowModal();
		wxBeginBusyCursor();

		if (res != wxID_OK)
			throw agi::UserCancelException("EBU/STL export");

		s.Save();
		return s;
	}

} // namespace {

Ebu3264SubtitleFormat::Ebu3264SubtitleFormat()
: SubtitleFormat("EBU subtitling data exchange format (EBU tech 3264, 1991)")
{
}

std::vector<std::string> Ebu3264SubtitleFormat::GetWriteWildcards() const
{
	std::vector<std::string> formats;
	formats.push_back("stl");
	return formats;
}

void Ebu3264SubtitleFormat::WriteFile(const AssFile *src, agi::fs::path const& filename, std::string const&) const
{
	// collect data from user
	EbuExportSettings export_settings = get_export_config(nullptr);
	AssFile copy(*src);

	std::vector<EbuSubtitle> subs_list = convert_subtitles(copy, export_settings);
	std::vector<BlockTTI> tti = create_blocks(subs_list, export_settings);
	BlockGSI gsi = create_header(copy, export_settings);

	BlockTTI &block0 = tti.front();
	snprintf(gsi.tcf, 8, "%02u%02u%02u%02u", (unsigned int)block0.tci.h, (unsigned int)block0.tci.m, (unsigned int)block0.tci.s, (unsigned int)block0.tci.f);
	snprintf(gsi.tnb, 5, "%5u", (unsigned int)tti.size());
	snprintf(gsi.tns, 5, "%5u", (unsigned int)subs_list.size());

	// write file
	agi::io::Save f(filename, true);
	f.Get().write((const char *)&gsi, sizeof(gsi));
	for (auto const& block : tti)
		f.Get().write((const char *)&block, sizeof(block));
}

max_stars_repo_name: rcombs/Aegisub

max_issues_repo_name: rcombs/Aegisub

max_forks_repo_name: rcombs/Aegisub


text:
#pragma once

#include <boost/interprocess/managed_mapped_file.hpp>
#include <boost/interprocess/containers/map.hpp>
#include <boost/interprocess/containers/set.hpp>
#include <boost/interprocess/containers/flat_map.hpp>
#include <boost/interprocess/containers/deque.hpp>
#include <boost/interprocess/containers/string.hpp>
#include <boost/interprocess/allocators/allocator.hpp>
#include <boost/interprocess/sync/interprocess_sharable_mutex.hpp>
#include <boost/interprocess/sync/sharable_lock.hpp>
#include <boost/interprocess/sync/file_lock.hpp>

#include <boost/multi_index_container.hpp>

#include <boost/chrono.hpp>
#include <boost/config.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/thread.hpp>
#include <boost/throw_exception.hpp>

#include <array>
#include <atomic>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <typeindex>
#include <typeinfo>

#ifndef CHAINBASE_NUM_RW_LOCKS
   #define CHAINBASE_NUM_RW_LOCKS 10
#endif

#ifdef CHAINBASE_CHECK_LOCKING
   #define CHAINBASE_REQUIRE_READ_LOCK(m, t) require_read_lock(m, typeid(t).name())
   #define CHAINBASE_REQUIRE_WRITE_LOCK(m, t) require_write_lock(m, typeid(t).name())
#else
   #define CHAINBASE_REQUIRE_READ_LOCK(m, t)
   #define CHAINBASE_REQUIRE_WRITE_LOCK(m, t)
#endif

namespace chainbase {

   namespace bip = boost::interprocess;
   namespace bfs = boost::filesystem;
   using std::unique_ptr;
   using std::vector;

   template<typename T>
   using allocator = bip::allocator<T, bip::managed_mapped_file::segment_manager>;

   typedef bip::basic_string< char, std::char_traits< char >, allocator< char > > shared_string;

   template<typename T>
   using shared_vector = std::vector<T, allocator<T> >;

   constexpr char _db_dirty_flag_string[] = "db_dirty_flag";

   struct strcmp_less
   {
      bool operator()( const shared_string& a, const shared_string& b )const
      {
         return less( a.c_str(), b.c_str() );
      }

      bool operator()( const shared_string& a, const std::string& b )const
      {
         return less( a.c_str(), b.c_str() );
      }

      bool operator()( const std::string& a, const shared_string& b )const
      {
         return less( a.c_str(), b.c_str() );
      }
      private:
         inline bool less( const char* a, const char* b )const
         {
            return std::strcmp( a, b ) < 0;
         }
   };

   typedef boost::interprocess::interprocess_sharable_mutex read_write_mutex;
   typedef boost::interprocess::sharable_lock< read_write_mutex > read_lock;
   typedef boost::unique_lock< read_write_mutex > write_lock;

   /**
    *  Object ID type that includes the type of the object it references
    */
   template<typename T>
   class oid {
      public:
         oid( int64_t i = 0 ):_id(i){}

         oid& operator++() { ++_id; return *this; }

         friend bool operator < ( const oid& a, const oid& b ) { return a._id < b._id; }
         friend bool operator > ( const oid& a, const oid& b ) { return a._id > b._id; }
         friend bool operator == ( const oid& a, const oid& b ) { return a._id == b._id; }
         friend bool operator != ( const oid& a, const oid& b ) { return a._id != b._id; }
         friend std::ostream& operator<<(std::ostream& s, const oid& id) {
            s << boost::core::demangle(typeid(oid<T>).name()) << '(' << id._id << ')'; return s;
         }

         int64_t _id = 0;
   };

   template<uint16_t TypeNumber, typename Derived>
   struct object
   {
      typedef oid<Derived> id_type;
      static const uint16_t type_id = TypeNumber;
   };

   /** this class is ment to be specified to enable lookup of index type by object type using
    * the SET_INDEX_TYPE macro.
    **/
   template<typename T>
   struct get_index_type {};

   /**
    *  This macro must be used at global scope and OBJECT_TYPE and INDEX_TYPE must be fully qualified
    */
   #define CHAINBASE_SET_INDEX_TYPE( OBJECT_TYPE, INDEX_TYPE )  \
   namespace chainbase { template<> struct get_index_type<OBJECT_TYPE> { typedef INDEX_TYPE type; }; }

   #define CHAINBASE_DEFAULT_CONSTRUCTOR( OBJECT_TYPE ) \
   template<typename Constructor, typename Allocator> \
   OBJECT_TYPE( Constructor&& c, Allocator&&  ) { c(*this); }

   template< typename value_type >
   class undo_state
   {
      public:
         typedef typename value_type::id_type                      id_type;
         typedef allocator< std::pair<const id_type, value_type> > id_value_allocator_type;
         typedef allocator< id_type >                              id_allocator_type;

         template<typename T>
         undo_state( allocator<T> al )
         :old_values( id_value_allocator_type( al.get_segment_manager() ) ),
          removed_values( id_value_allocator_type( al.get_segment_manager() ) ),
          new_ids( id_allocator_type( al.get_segment_manager() ) ){}

         typedef boost::interprocess::map< id_type, value_type, std::less<id_type>, id_value_allocator_type >  id_value_type_map;
         typedef boost::interprocess::set< id_type, std::less<id_type>, id_allocator_type >                    id_type_set;

         id_value_type_map            old_values;
         id_value_type_map            removed_values;
         id_type_set                  new_ids;
         id_type                      old_next_id = 0;
         int64_t                      revision = 0;
   };

   /**
    * The code we want to implement is this:
    *
    * ++target; try { ... } finally { --target }
    *
    * In C++ the only way to implement finally is to create a class
    * with a destructor, so that's what we do here.
    */
   class int_incrementer
   {
      public:
         int_incrementer( int32_t& target ) : _target(target)
         { ++_target; }
         ~int_incrementer()
         { --_target; }

         int32_t get()const
         { return _target; }

      private:
         int32_t& _target;
   };

   /**
    *  The value_type stored in the multiindex container must have a integer field with the name 'id'.  This will
    *  be the primary key and it will be assigned and managed by generic_index.
    *
    *  Additionally, the constructor for value_type must take an allocator
    */
   template<typename MultiIndexType>
   class generic_index
   {
      public:
         typedef bip::managed_mapped_file::segment_manager             segment_manager_type;
         typedef MultiIndexType                                        index_type;
         typedef typename index_type::value_type                       value_type;
         typedef bip::allocator< generic_index, segment_manager_type > allocator_type;
         typedef undo_state< value_type >                              undo_state_type;

         generic_index( allocator<value_type> a )
         :_stack(a),_indices( a ),_size_of_value_type( sizeof(typename MultiIndexType::node_type) ),_size_of_this(sizeof(*this)){}

         void validate()const {
            if( sizeof(typename MultiIndexType::node_type) != _size_of_value_type || sizeof(*this) != _size_of_this )
               BOOST_THROW_EXCEPTION( std::runtime_error("content of memory does not match data expected by executable") );
         }

         /**
          * Construct a new element in the multi_index_container.
          * Set the ID to the next available ID, then increment _next_id and fire off on_create().
          */
         template<typename Constructor>
         const value_type& emplace( Constructor&& c ) {
            auto new_id = _next_id;

            auto constructor = [&]( value_type& v ) {
               v.id = new_id;
               c( v );
            };

            auto insert_result = _indices.emplace( constructor, _indices.get_allocator() );

            if( !insert_result.second ) {
               BOOST_THROW_EXCEPTION( std::logic_error("could not insert object, most likely a uniqueness constraint was violated") );
            }

            ++_next_id;
            on_create( *insert_result.first );
            return *insert_result.first;
         }

         template<typename Modifier>
         void modify( const value_type& obj, Modifier&& m ) {
            on_modify( obj );
            auto ok = _indices.modify( _indices.iterator_to( obj ), m );
            if( !ok ) BOOST_THROW_EXCEPTION( std::logic_error( "Could not modify object, most likely a uniqueness constraint was violated" ) );
         }

         void remove( const value_type& obj ) {
            on_remove( obj );
            _indices.erase( _indices.iterator_to( obj ) );
         }

         template<typename CompatibleKey>
         const value_type* find( CompatibleKey&& key )const {
            auto itr = _indices.find( std::forward<CompatibleKey>(key) );
            if( itr != _indices.end() ) return &*itr;
            return nullptr;
         }

         template<typename CompatibleKey>
         const value_type& get( CompatibleKey&& key )const {
            auto ptr = find( key );
            if( !ptr ) {
               std::stringstream ss;
               ss << "key not found (" << boost::core::demangle( typeid( key ).name() ) << "): " << key;
               BOOST_THROW_EXCEPTION( std::out_of_range( ss.str().c_str() ) );
            }
            return *ptr;
         }

         const index_type& indices()const { return _indices; }

         class session {
            public:
               session( session&& mv )
               :_index(mv._index),_apply(mv._apply){ mv._apply = false; }

               ~session() {
                  if( _apply ) {
                     _index.undo();
                  }
               }

               /** leaves the UNDO state on the stack when session goes out of scope */
               void push()   { _apply = false; }
               /** combines this session with the prior session */
               void squash() { if( _apply ) _index.squash(); _apply = false; }
               void undo()   { if( _apply ) _index.undo();  _apply = false; }

               session& operator = ( session&& mv ) {
                  if( this == &mv ) return *this;
                  if( _apply ) _index.undo();
                  _apply = mv._apply;
                  mv._apply = false;
                  return *this;
               }

               int64_t revision()const { return _revision; }

            private:
               friend class generic_index;

               session( generic_index& idx, int64_t revision )
               :_index(idx),_revision(revision) {
                  if( revision == -1 )
                     _apply = false;
               }

               generic_index& _index;
               bool           _apply = true;
               int64_t        _revision = 0;
         };

         session start_undo_session( bool enabled ) {
            if( enabled ) {
               _stack.emplace_back( _indices.get_allocator() );
               _stack.back().old_next_id = _next_id;
               _stack.back().revision = ++_revision;
               return session( *this, _revision );
            } else {
               return session( *this, -1 );
            }
         }

         int64_t revision()const { return _revision; }

         /**
          *  Restores the state to how it was prior to the current session discarding all changes
          *  made between the last revision and the current revision.
          */
         void undo() {
            if( !enabled() ) return;

            const auto& head = _stack.back();

            // 还原改变的数据
            for( auto& item : head.old_values ) {
               auto ok = _indices.modify( _indices.find( item.second.id ), [&]( value_type& v ) {
                  v = std::move( item.second );
               });
               if( !ok ) BOOST_THROW_EXCEPTION( std::logic_error( "Could not modify object, most likely a uniqueness constraint was violated" ) );
            }

            // 删除新增的数据
            for( auto id : head.new_ids )
            {
               _indices.erase( _indices.find( id ) );
            }
            _next_id = head.old_next_id;

            // 恢复删除的数据
            for( auto& item : head.removed_values ) {
               bool ok = _indices.emplace( std::move( item.second ) ).second;
               if( !ok ) BOOST_THROW_EXCEPTION( std::logic_error( "Could not restore object, most likely a uniqueness constraint was violated" ) );
            }

            _stack.pop_back();
            --_revision;
         }

         /**
          *  This method works similar to git squash, it merges the change set from the two most
          *  recent revision numbers into one revision number (reducing the head revision number)
          *
          *  This method does not change the state of the index, only the state of the undo buffer.
          *
          *  合并_stack顶层的两个undo_state
          */
         void squash()
         {
            if( !enabled() ) return;
            if( _stack.size() == 1 ) {
               _stack.pop_front();
               --_revision;
               return;
            }

            auto& state = _stack.back();
            auto& prev_state = _stack[_stack.size()-2];

            // An object's relationship to a state can be:
            // in new_ids            : new
            // in old_values (was=X) : upd(was=X)
            // in removed (was=X)    : del(was=X)
            // not in any of above   : nop
            //
            // When merging A=prev_state and B=state we have a 4x4 matrix of all possibilities:
            //
            //                   |--------------------- B ----------------------|
            //
            //                +------------+------------+------------+------------+
            //                | new        | upd(was=Y) | del(was=Y) | nop        |
            //   +------------+------------+------------+------------+------------+
            // / | new        | N/A        | new       A| nop       C| new       A|
            // | +------------+------------+------------+------------+------------+
            // | | upd(was=X) | N/A        | upd(was=X)A| del(was=X)C| upd(was=X)A|
            // A +------------+------------+------------+------------+------------+
            // | | del(was=X) | N/A        | N/A        | N/A        | del(was=X)A|
            // | +------------+------------+------------+------------+------------+
            // \ | nop        | new       B| upd(was=Y)B| del(was=Y)B| nop      AB|
            //   +------------+------------+------------+------------+------------+
            //
            // Each entry was composed by labelling what should occur in the given case.
            //
            // Type A means the composition of states contains the same entry as the first of the two merged states for that object.
            // Type B means the composition of states contains the same entry as the second of the two merged states for that object.
            // Type C means the composition of states contains an entry different from either of the merged states for that object.
            // Type N/A means the composition of states violates causal timing.
            // Type AB means both type A and type B simultaneously.
            //
            // The merge() operation is defined as modifying prev_state in-place to be the state object which represents the composition of
            // state A and B.
            //
            // Type A (and AB) can be implemented as a no-op; prev_state already contains the correct value for the merged state.
            // Type B (and AB) can be implemented by copying from state to prev_state.
            // Type C needs special case-by-case logic.
            // Type N/A can be ignored or assert(false) as it can only occur if prev_state and state have illegal values
            // (a serious logic error which should never happen).
            //

            // We can only be outside type A/AB (the nop path) if B is not nop, so it suffices to iterate through B's three containers.

            for( const auto& item : state.old_values )
            {
               if( prev_state.new_ids.find( item.second.id ) != prev_state.new_ids.end() )
               {
                  // new+upd -> new, type A
                  continue;
               }
               if( prev_state.old_values.find( item.second.id ) != prev_state.old_values.end() )
               {
                  // upd(was=X) + upd(was=Y) -> upd(was=X), type A
                  continue;
               }
               // del+upd -> N/A
               assert( prev_state.removed_values.find(item.second.id) == prev_state.removed_values.end() );
               // nop+upd(was=Y) -> upd(was=Y), type B
               prev_state.old_values.emplace( std::move(item) );
            }

            // *+new, but we assume the N/A cases don't happen, leaving type B nop+new -> new
            for( auto id : state.new_ids )
               prev_state.new_ids.insert(id);

            // *+del
            for( auto& obj : state.removed_values )
            {
               if( prev_state.new_ids.find(obj.second.id) != prev_state.new_ids.end() )
               {
                  // new + del -> nop (type C)
                  prev_state.new_ids.erase(obj.second.id);
                  continue;
               }
               auto it = prev_state.old_values.find(obj.second.id);
               if( it != prev_state.old_values.end() )
               {
                  // upd(was=X) + del(was=Y) -> del(was=X)
                  prev_state.removed_values.emplace( std::move(*it) );
                  prev_state.old_values.erase(obj.second.id);
                  continue;
               }
               // del + del -> N/A
               assert( prev_state.removed_values.find( obj.second.id ) == prev_state.removed_values.end() );
               // nop + del(was=Y) -> del(was=Y)
               prev_state.removed_values.emplace( std::move(obj) ); //[obj.second->id] = std::move(obj.second);
            }

            _stack.pop_back();
            --_revision;
         }

         /**
          * Discards all undo history prior to revision
          */
         void commit( int64_t revision )
         {
            while( _stack.size() && _stack[0].revision <= revision )
            {
               _stack.pop_front();
            }
         }

         /**
          * Unwinds all undo states
          */
         void undo_all()
         {
            while( enabled() )
               undo();
         }

         void set_revision( uint64_t revision )
         {
            if( _stack.size() != 0 )
               BOOST_THROW_EXCEPTION( std::logic_error("cannot set revision while there is an existing undo stack") );

            if( revision > std::numeric_limits<int64_t>::max() )
               BOOST_THROW_EXCEPTION( std::logic_error("revision to set is too high") );

            _revision = static_cast<int64_t>(revision);
         }

         void remove_object( int64_t id )
         {
            const value_type* val = find( typename value_type::id_type(id) );
            if( !val ) BOOST_THROW_EXCEPTION( std::out_of_range( boost::lexical_cast<std::string>(id) ) );
            remove( *val );
         }

         std::pair<int64_t, int64_t> undo_stack_revision_range()const {
            int64_t begin = _revision;
            int64_t end   = _revision;

            if( _stack.size() > 0 ) {
               begin = _stack.front().revision - 1;
               end   = _stack.back().revision;
            }

            return {begin, end};
         }

         const auto& stack()const { return _stack; }

      private:
         bool enabled()const { return _stack.size(); }

         void on_modify( const value_type& v ) {
            if( !enabled() ) return;

            auto& head = _stack.back();

            if( head.new_ids.find( v.id ) != head.new_ids.end() )
               return;

            auto itr = head.old_values.find( v.id );
            if( itr != head.old_values.end() )
               return;

            head.old_values.emplace( std::pair< typename value_type::id_type, const value_type& >( v.id, v ) );
         }

         void on_remove( const value_type& v ) {
            if( !enabled() ) return;

            auto& head = _stack.back();
            if( head.new_ids.count(v.id) ) {
               head.new_ids.erase( v.id );
               return;
            }

            auto itr = head.old_values.find( v.id );
            if( itr != head.old_values.end() ) {
               head.removed_values.emplace( std::move( *itr ) );
               head.old_values.erase( v.id );
               return;
            }

            if( head.removed_values.count( v.id ) )
               return;

            head.removed_values.emplace( std::pair< typename value_type::id_type, const value_type& >( v.id, v ) );
         }

         void on_create( const value_type& v ) {
            if( !enabled() ) return;
            auto& head = _stack.back();

            head.new_ids.insert( v.id );
         }

         boost::interprocess::deque< undo_state_type, allocator<undo_state_type> > _stack;

         /**
          *  Each new session increments the revision, a squash will decrement the revision by combining
          *  the two most recent revisions into one revision.
          *
          *  Commit will discard all revisions prior to the committed revision.
          */
         int64_t                         _revision = 0;
         typename value_type::id_type    _next_id = 0;
         index_type                      _indices;
         uint32_t                        _size_of_value_type = 0;
         uint32_t                        _size_of_this = 0;
   };

   class abstract_session {
      public:
         virtual ~abstract_session(){};
         virtual void push()             = 0;
         virtual void squash()           = 0;
         virtual void undo()             = 0;
         virtual int64_t revision()const  = 0;
   };

   template<typename SessionType>
   class session_impl : public abstract_session
   {
      public:
         session_impl( SessionType&& s ):_session( std::move( s ) ){}

         virtual void push() override  { _session.push();  }
         virtual void squash() override{ _session.squash(); }
         virtual void undo() override  { _session.undo();  }
         virtual int64_t revision()const override  { return _session.revision();  }
      private:
         SessionType _session;
   };

   class abstract_index
   {
      public:
         abstract_index( void* i ):_idx_ptr(i){}
         virtual ~abstract_index(){}
         virtual void     set_revision( uint64_t revision ) = 0;
         virtual unique_ptr<abstract_session> start_undo_session( bool enabled ) = 0;

         virtual int64_t revision()const = 0;
         virtual void    undo()const = 0;
         virtual void    squash()const = 0;
         virtual void    commit( int64_t revision )const = 0;
         virtual void    undo_all()const = 0;
         virtual uint32_t type_id()const  = 0;
         virtual uint64_t row_count()const = 0;
         virtual const std::string& type_name()const = 0;
         virtual std::pair<int64_t, int64_t> undo_stack_revision_range()const = 0;

         virtual void remove_object( int64_t id ) = 0;

         void* get()const { return _idx_ptr; }
      private:
         void* _idx_ptr;
   };

   template<typename BaseIndex>
   class index_impl : public abstract_index {
      public:
         index_impl( BaseIndex& base ):abstract_index( &base ),_base(base){}

         virtual unique_ptr<abstract_session> start_undo_session( bool enabled ) override {
            return unique_ptr<abstract_session>(new session_impl<typename BaseIndex::session>( _base.start_undo_session( enabled ) ) );
         }

         virtual void     set_revision( uint64_t revision ) override { _base.set_revision( revision ); }
         virtual int64_t  revision()const  override { return _base.revision(); }
         virtual void     undo()const  override { _base.undo(); }
         virtual void     squash()const  override { _base.squash(); }
         virtual void     commit( int64_t revision )const  override { _base.commit(revision); }
         virtual void     undo_all() const override {_base.undo_all(); }
         virtual uint32_t type_id()const override { return BaseIndex::value_type::type_id; }
         virtual uint64_t row_count()const override { return _base.indices().size(); }
         virtual const std::string& type_name() const override { return BaseIndex_name; }
         virtual std::pair<int64_t, int64_t> undo_stack_revision_range()const override { return _base.undo_stack_revision_range(); }

         virtual void     remove_object( int64_t id ) override { return _base.remove_object( id ); }
      private:
         BaseIndex& _base;
         std::string BaseIndex_name = boost::core::demangle( typeid( typename BaseIndex::value_type ).name() );
   };

   template<typename IndexType>
   class index : public index_impl<IndexType> {
      public:
         index( IndexType& i ):index_impl<IndexType>( i ){}
   };

   class read_write_mutex_manager
   {
      public:
         read_write_mutex_manager()
         {
            _current_lock = 0;
         }

         ~read_write_mutex_manager(){}

         void next_lock()
         {
            _current_lock++;
            new( &_locks[ _current_lock % CHAINBASE_NUM_RW_LOCKS ] ) read_write_mutex();
         }

         read_write_mutex& current_lock()
         {
            return _locks[ _current_lock % CHAINBASE_NUM_RW_LOCKS ];
         }

         uint32_t current_lock_num()
         {
            return _current_lock;
         }

      private:
         std::array< read_write_mutex, CHAINBASE_NUM_RW_LOCKS >     _locks;
         std::atomic< uint32_t >                                    _current_lock;
   };

   /**
    *  This class
    */
   class database
   {
      public:
         enum open_flags {
            read_only     = 0, // 只读权限
            read_write    = 1 // 读写权限
         };

         using database_index_row_count_multiset = std::multiset<std::pair<unsigned, std::string>>;

         database(const bfs::path& dir, open_flags write = read_only, uint64_t shared_file_size = 0, bool allow_dirty = false);
         ~database();
         database(database&&) = default;
         database& operator=(database&&) = default;
         bool is_read_only() const { return _read_only; }
         void flush();
         void set_require_locking( bool enable_require_locking );

#ifdef CHAINBASE_CHECK_LOCKING
         void require_lock_fail( const char* method, const char* lock_type, const char* tname )const;

         void require_read_lock( const char* method, const char* tname )const
         {
            if( BOOST_UNLIKELY( _enable_require_locking & _read_only & (_read_lock_count <= 0) ) )
               require_lock_fail(method, "read", tname);
         }

         void require_write_lock( const char* method, const char* tname )
         {
            if( BOOST_UNLIKELY( _enable_require_locking & (_write_lock_count <= 0) ) )
               require_lock_fail(method, "write", tname);
         }
#endif

         struct session {
            public:
               session( session&& s ):_index_sessions( std::move(s._index_sessions) ),_revision( s._revision ){}
               session( vector<std::unique_ptr<abstract_session>>&& s ):_index_sessions( std::move(s) )
               {
                  if( _index_sessions.size() )
                     _revision = _index_sessions[0]->revision();
               }

               ~session() {
                  undo();
               }

               void push()
               {
                  for( auto& i : _index_sessions ) i->push();
                  _index_sessions.clear();
               }

               void squash()
               {
                  for( auto& i : _index_sessions ) i->squash();
                  _index_sessions.clear();
               }

               void undo()
               {
                  for( auto& i : _index_sessions ) i->undo();
                  _index_sessions.clear();
               }

               int64_t revision()const { return _revision; }

            private:
               friend class database;
               session(){}

               vector< std::unique_ptr<abstract_session> > _index_sessions;
               int64_t _revision = -1;
         };

         session start_undo_session( bool enabled );

         int64_t revision()const {
             if( _index_list.size() == 0 ) return -1;
             return _index_list[0]->revision();
         }

         void undo();
         void squash();
         void commit( int64_t revision );
         void undo_all();

         void set_revision( uint64_t revision )
         {
             CHAINBASE_REQUIRE_WRITE_LOCK( "set_revision", uint64_t );
             for( auto i : _index_list ) i->set_revision( revision );
         }

         template<typename MultiIndexType>
         void add_index() {
            const uint16_t type_id = generic_index<MultiIndexType>::value_type::type_id;
            typedef generic_index<MultiIndexType>          index_type;
            typedef typename index_type::allocator_type    index_alloc;

            std::string type_name = boost::core::demangle( typeid( typename index_type::value_type ).name() );

            if( !( _index_map.size() <= type_id || _index_map[ type_id ] == nullptr ) ) {
               BOOST_THROW_EXCEPTION( std::logic_error( type_name + "::type_id is already in use" ) );
            }

            index_type* idx_ptr = _segment->find< index_type >( type_name.c_str() ).first;
            bool first_time_adding = false;
            if( !idx_ptr ) {
               if( _read_only ) {
                  BOOST_THROW_EXCEPTION( std::runtime_error( "unable to find index for " + type_name + " in read only database" ) );
               }
               first_time_adding = true;
               idx_ptr = _segment->construct< index_type >( type_name.c_str() )( index_alloc( _segment->get_segment_manager() ) );
             }

            idx_ptr->validate();

            // Ensure the undo stack of added index is consistent with the other indices in the database
            if( _index_list.size() > 0 ) {
               auto expected_revision_range = _index_list.front()->undo_stack_revision_range();
               auto added_index_revision_range = idx_ptr->undo_stack_revision_range();

               if( added_index_revision_range.first != expected_revision_range.first ||
                   added_index_revision_range.second != expected_revision_range.second ) {

                  if( !first_time_adding ) {
                     BOOST_THROW_EXCEPTION( std::logic_error(
                        "existing index for " + type_name + " has an undo stack (revision range [" +
                        std::to_string(added_index_revision_range.first) + ", " + std::to_string(added_index_revision_range.second) +
                        "]) that is inconsistent with other indices in the database (revision range [" +
                        std::to_string(expected_revision_range.first) + ", " + std::to_string(expected_revision_range.second) +
                        "]); corrupted database?"
                     ) );
                  }

                  if( _read_only ) {
                     BOOST_THROW_EXCEPTION( std::logic_error(
                        "new index for " + type_name +
                        " requires an undo stack that is consistent with other indices in the database; cannot fix in read-only mode"
                     ) );
                  }

                  idx_ptr->set_revision( static_cast<uint64_t>(expected_revision_range.first) );
                  while( idx_ptr->revision() < expected_revision_range.second ) {
                     idx_ptr->start_undo_session(true).push();
                  }
               }
            }

            if( type_id >= _index_map.size() )
               _index_map.resize( type_id + 1 );

            auto new_index = new index<index_type>( *idx_ptr );
            _index_map[ type_id ].reset( new_index );
            _index_list.push_back( new_index );
         }

         auto get_segment_manager() -> decltype( ((bip::managed_mapped_file*)nullptr)->get_segment_manager()) {
            return _segment->get_segment_manager();
         }

         auto get_segment_manager()const -> std::add_const_t< decltype( ((bip::managed_mapped_file*)nullptr)->get_segment_manager() ) > {
            return _segment->get_segment_manager();
         }

         size_t get_free_memory()const
         {
            return _segment->get_segment_manager()->get_free_memory();
         }

         template<typename MultiIndexType>
         const generic_index<MultiIndexType>& get_index()const
         {
            CHAINBASE_REQUIRE_READ_LOCK("get_index", typename MultiIndexType::value_type);
            typedef generic_index<MultiIndexType> index_type;
            typedef index_type*                   index_type_ptr;
            assert( _index_map.size() > index_type::value_type::type_id );
            assert( _index_map[index_type::value_type::type_id] );
            return *index_type_ptr( _index_map[index_type::value_type::type_id]->get() );
         }

         template<typename MultiIndexType, typename ByIndex>
         auto get_index()const -> decltype( ((generic_index<MultiIndexType>*)( nullptr ))->indices().template get<ByIndex>() )
         {
            CHAINBASE_REQUIRE_READ_LOCK("get_index", typename MultiIndexType::value_type);
            typedef generic_index<MultiIndexType> index_type;
            typedef index_type*                   index_type_ptr;
            assert( _index_map.size() > index_type::value_type::type_id );
            assert( _index_map[index_type::value_type::type_id] );
            return index_type_ptr( _index_map[index_type::value_type::type_id]->get() )->indices().template get<ByIndex>();
         }

         template<typename MultiIndexType>
         generic_index<MultiIndexType>& get_mutable_index()
         {
            CHAINBASE_REQUIRE_WRITE_LOCK("get_mutable_index", typename MultiIndexType::value_type);
            typedef generic_index<MultiIndexType> index_type;
            typedef index_type*                   index_type_ptr;
            assert( _index_map.size() > index_type::value_type::type_id );
            assert( _index_map[index_type::value_type::type_id] );
            return *index_type_ptr( _index_map[index_type::value_type::type_id]->get() );
         }

         template< typename ObjectType, typename IndexedByType, typename CompatibleKey >
         const ObjectType* find( CompatibleKey&& key )const
         {
             CHAINBASE_REQUIRE_READ_LOCK("find", ObjectType);
             typedef typename get_index_type< ObjectType >::type index_type;
             const auto& idx = get_index< index_type >().indices().template get< IndexedByType >();
             auto itr = idx.find( std::forward< CompatibleKey >( key ) );
             if( itr == idx.end() ) return nullptr;
             return &*itr;
         }

         template< typename ObjectType >
         const ObjectType* find( oid< ObjectType > key = oid< ObjectType >() ) const
         {
             CHAINBASE_REQUIRE_READ_LOCK("find", ObjectType);
             typedef typename get_index_type< ObjectType >::type index_type;
             const auto& idx = get_index< index_type >().indices();
             auto itr = idx.find( key );
             if( itr == idx.end() ) return nullptr;
             return &*itr;
         }

         template< typename ObjectType, typename IndexedByType, typename CompatibleKey >
         const ObjectType& get( CompatibleKey&& key )const
         {
             CHAINBASE_REQUIRE_READ_LOCK("get", ObjectType);
             auto obj = find< ObjectType, IndexedByType >( std::forward< CompatibleKey >( key ) );
             if( !obj ) {
                std::stringstream ss;
                ss << "unknown key (" << boost::core::demangle( typeid( key ).name() ) << "): " << key;
                BOOST_THROW_EXCEPTION( std::out_of_range( ss.str().c_str() ) );
             }
             return *obj;
         }

         template< typename ObjectType >
         const ObjectType& get( const oid< ObjectType >& key = oid< ObjectType >() )const
         {
             CHAINBASE_REQUIRE_READ_LOCK("get", ObjectType);
             auto obj = find< ObjectType >( key );
             if( !obj ) {
                std::stringstream ss;
                ss << "unknown key (" << boost::core::demangle( typeid( key ).name() ) << "): " << key._id;
                BOOST_THROW_EXCEPTION( std::out_of_range( ss.str().c_str() ) );
             }
             return *obj;
         }

         template<typename ObjectType, typename Modifier>
         void modify( const ObjectType& obj, Modifier&& m )
         {
             CHAINBASE_REQUIRE_WRITE_LOCK("modify", ObjectType);
             typedef typename get_index_type<ObjectType>::type index_type;
             get_mutable_index<index_type>().modify( obj, m );
         }

         template<typename ObjectType>
         void remove( const ObjectType& obj )
         {
             CHAINBASE_REQUIRE_WRITE_LOCK("remove", ObjectType);
             typedef typename get_index_type<ObjectType>::type index_type;
             return get_mutable_index<index_type>().remove( obj );
         }

         template<typename ObjectType, typename Constructor>
         const ObjectType& create( Constructor&& con )
         {
             CHAINBASE_REQUIRE_WRITE_LOCK("create", ObjectType);
             typedef typename get_index_type<ObjectType>::type index_type;
             return get_mutable_index<index_type>().emplace( std::forward<Constructor>(con) );
         }

         template< typename Lambda >
         auto with_read_lock( Lambda&& callback, uint64_t wait_micro = 1000000 ) const -> decltype( (*(Lambda*)nullptr)() )
         {
            read_lock lock( _rw_manager->current_lock(), bip::defer_lock_type() );
#ifdef CHAINBASE_CHECK_LOCKING
            BOOST_ATTRIBUTE_UNUSED
            int_incrementer ii( _read_lock_count );
#endif

            if( !wait_micro )
            {
               lock.lock();
            }
            else
            {

               if( !lock.timed_lock( boost::posix_time::microsec_clock::local_time() + boost::posix_time::microseconds( wait_micro ) ) )
                  BOOST_THROW_EXCEPTION( std::runtime_error( "unable to acquire lock" ) );
            }

            return callback();
         }

         template< typename Lambda >
         auto with_write_lock( Lambda&& callback, uint64_t wait_micro = 1000000 ) -> decltype( (*(Lambda*)nullptr)() )
         {
            if( _read_only )
               BOOST_THROW_EXCEPTION( std::logic_error( "cannot acquire write lock on read-only process" ) );

            write_lock lock( _rw_manager->current_lock(), boost::defer_lock_t() );
#ifdef CHAINBASE_CHECK_LOCKING
            BOOST_ATTRIBUTE_UNUSED
            int_incrementer ii( _write_lock_count );
#endif

            if( !wait_micro )
            {
               lock.lock();
            }
            else
            {
               while( !lock.timed_lock( boost::posix_time::microsec_clock::local_time() + boost::posix_time::microseconds( wait_micro ) ) )
               {
                  _rw_manager->next_lock();
                  std::cerr << "Lock timeout, moving to lock " << _rw_manager->current_lock_num() << std::endl;
                  lock = write_lock( _rw_manager->current_lock(), boost::defer_lock_t() );
               }
            }

            return callback();
         }

         database_index_row_count_multiset row_count_per_index()const {
            database_index_row_count_multiset ret;
            for(const auto& ai_ptr : _index_map) {
               if(!ai_ptr)
                  continue;
               ret.emplace(make_pair(ai_ptr->row_count(), ai_ptr->type_name()));
            }
            return ret;
         }

      private:
         unique_ptr<bip::managed_mapped_file>                        _segment;
         unique_ptr<bip::managed_mapped_file>                        _meta;
         read_write_mutex_manager*                                   _rw_manager = nullptr;
         bool                                                        _read_only = false;
         bip::file_lock                                              _flock;

         /**
          * This is a sparse list of known indices kept to accelerate creation of undo sessions
          */
         vector<abstract_index*>                                     _index_list;

         /**
          * This is a full map (size 2^16) of all possible index designed for constant time lookup
          */
         vector<unique_ptr<abstract_index>>                          _index_map;

         bfs::path                                                   _data_dir;

#ifdef CHAINBASE_CHECK_LOCKING
         int32_t                                                     _read_lock_count = 0;
         int32_t                                                     _write_lock_count = 0;
         bool                                                        _enable_require_locking = false;
#endif

         void                                                        _msync_database();
   };

   template<typename Object, typename... Args>
   using shared_multi_index_container = boost::multi_index_container<Object,Args..., chainbase::allocator<Object> >;
}  // namepsace chainbase

ext: hpp

max_stars_repo_name: happyscene/learneos

max_issues_repo_name: happyscene/learneos

max_forks_repo_name: happyscene/learneos


text:
#include <iostream>
#include <algorithm>
#include <cmath>
#include <vector>
#include <map>
#include <set>
#include <tuple>
#include <stdbool.h>
#include <bitset>
#include <string>
#include <boost/multiprecision/cpp_int.hpp>

namespace mp = boost::multiprecision;
using namespace std;

int mizuyari(vector<int> h) {
    // 0が含まれない区間で分割統治する
    // for(const auto& x : h) {
    //     cout << x << ", ";
    // }
    // cout << endl;

    int count = 0;
    for(int i = 0 ; i < h.size() ; ) {
        if (!h[i]) {
            i++;
            continue;
        }
        vector<int> sub; // 分割統治するvector
        int j;
        for(j = i ; j < h.size() ; ++j ) {
            if(!h[j]) {
                break;
            } else {
                sub.push_back(h[j] - 1); // 水やりする
            }
        }
        // (i, j]区間でいける
        count++; // 今回の水やり分
        count += mizuyari(sub); // 水をやった後の残りはまかした
        // 後処理
        i = j;
    }
    return count;
}

int main(void) {
    int n;
    cin >> n;
    vector<int> h(n);
    for(int i = 0 ; i < n ; ++i) {
        cin >> h[i];
    }
    auto ans = mizuyari(h);
    cout << ans << endl;
    return 0;
}

ext: cpp

max_stars_repo_name: kamiyaowl/atcoder

max_issues_repo_name: kamiyaowl/atcoder

max_forks_repo_name: kamiyaowl/atcoder


text:
#define BOOST_TEST_MODULE matrix

#include <boost/test/unit_test.hpp>
#include <boost/test/test_case_template.hpp>
#include <boost/mpl/list.hpp>

#include <mla/matrix/all.h++>

typedef boost::mpl::list<
	mla::matrix::DenseRowMajor<float>,
	mla::matrix::DenseRowMajor<double>,
	mla::matrix::Diagonal<float>,
	mla::matrix::Diagonal<double>,
	mla::matrix::SparseDOK<float>,
	mla::matrix::SparseDOK<double>,
	mla::matrix::SparseCRS<float>,
	mla::matrix::SparseCRS<double>
> matrix_type_list;

BOOST_AUTO_TEST_SUITE(test_matrix)

BOOST_AUTO_TEST_CASE_TEMPLATE( Declaration, MatrixType, matrix_type_list )
{
	size_t row_size = 3;
	size_t column_size = 3;

	MatrixType m(row_size,column_size);

	BOOST_CHECK_EQUAL(m.rows(), row_size);
	BOOST_CHECK_EQUAL(m.columns(), column_size);
}

BOOST_AUTO_TEST_CASE_TEMPLATE( Clearing, MatrixType, matrix_type_list )
{
	size_t row_size = 3;
	size_t column_size = 3;

	MatrixType m(row_size,column_size);

	m.setZero();

	for(unsigned int i = 0; i < row_size; i++)
	{
		for(unsigned int j = 0; j < column_size; j++)
		{
			BOOST_CHECK_EQUAL(m.getValue(i, j), (typename MatrixType::scalar_type)0);
		}
	}

}

BOOST_AUTO_TEST_CASE_TEMPLATE( set_eye_square, MatrixType, matrix_type_list )
{
	size_t row_size = 3;
	size_t column_size = 3;

	MatrixType m(row_size,column_size);

	m.setEye();

	for(unsigned int i = 0; i < row_size; i++)
	{
		for(unsigned int j = 0; j < column_size; j++)
		{
			BOOST_CHECK_EQUAL(m.getValue(i, j), (typename MatrixType::scalar_type)(i==j?1:0) );
		}
	}

}

BOOST_AUTO_TEST_CASE_TEMPLATE( set_eye_rectangular_row, MatrixType, matrix_type_list )
{
	size_t row_size = 3;
	size_t column_size = 5;

	MatrixType m(row_size,column_size);

	m.setEye();

	for(unsigned int i = 0; i < row_size; i++)
	{
		for(unsigned int j = 0; j < column_size; j++)
		{
			BOOST_CHECK_EQUAL(m.getValue(i, j), (typename MatrixType::scalar_type)(i==j?1:0) );
		}
	}

}

BOOST_AUTO_TEST_CASE_TEMPLATE( set_eye_rectangular_column, MatrixType, matrix_type_list )
{
	size_t row_size = 5;
	size_t column_size = 3;

	MatrixType m(row_size, column_size);

	m.setEye();

	for(unsigned int i = 0; i < row_size; i++)
	{
		for(unsigned int j = 0; j < column_size; j++)
		{
			BOOST_CHECK_EQUAL(m.getValue(i, j), (typename MatrixType::scalar_type)(i==j?1:0) );
		}
	}

}

BOOST_AUTO_TEST_CASE_TEMPLATE( reference_assigning_one, MatrixType, matrix_type_list )
{
	if( mla::matrix::Traits<MatrixType>::is_writeable() )
	{
		MatrixType m(3,3);
		m.setZero();

		m(1,1) = 1.0;

		BOOST_CHECK_EQUAL(m.getValue(0, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(0, 1), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(0, 2), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 1), 1.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 2), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 1), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 2), 0.0);
	}
}

BOOST_AUTO_TEST_CASE_TEMPLATE( reference_assigning_same_row_ordered, MatrixType, matrix_type_list )
{
	if( mla::matrix::Traits<MatrixType>::is_writeable() )
	{
		MatrixType m(3,3);
		m.setZero();

		m(1,0) = 1.0;
		m(1,1) = 1.0;

		BOOST_CHECK_EQUAL(m.getValue(0, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(0, 1), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(0, 2), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 0), 1.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 1), 1.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 2), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 1), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 2), 0.0);
	}
}

BOOST_AUTO_TEST_CASE_TEMPLATE( reference_assigning_same_row_unordered, MatrixType, matrix_type_list )
{
	if( mla::matrix::Traits<MatrixType>::is_writeable() )
	{
		MatrixType m(3,3);
		m.setZero();

		m(1,1) = 1.0;
		m(1,0) = 1.0;

		BOOST_CHECK_EQUAL(m.getValue(0, 0), 0);
		BOOST_CHECK_EQUAL(m.getValue(0, 1), 0);
		BOOST_CHECK_EQUAL(m.getValue(0, 2), 0);
		BOOST_CHECK_EQUAL(m.getValue(1, 0), 1);
		BOOST_CHECK_EQUAL(m.getValue(1, 1), 1);
		BOOST_CHECK_EQUAL(m.getValue(1, 2), 0);
		BOOST_CHECK_EQUAL(m.getValue(2, 0), 0);
		BOOST_CHECK_EQUAL(m.getValue(2, 1), 0);
		BOOST_CHECK_EQUAL(m.getValue(2, 2), 0);
	}
}

BOOST_AUTO_TEST_CASE_TEMPLATE( reference_assigning_different_row_ordered, MatrixType, matrix_type_list )
{
	if( mla::matrix::Traits<MatrixType>::is_writeable() )
	{
		MatrixType m(3,3);
		m.setZero();

		m(1,1) = 1.0;
		m(2,1) = 1.0;

		BOOST_CHECK_EQUAL(m.getValue(0, 0), 0);
		BOOST_CHECK_EQUAL(m.getValue(0, 1), 0);
		BOOST_CHECK_EQUAL(m.getValue(0, 2), 0);
		BOOST_CHECK_EQUAL(m.getValue(1, 0), 0);
		BOOST_CHECK_EQUAL(m.getValue(1, 1), 1);
		BOOST_CHECK_EQUAL(m.getValue(1, 2), 0);
		BOOST_CHECK_EQUAL(m.getValue(2, 0), 0);
		BOOST_CHECK_EQUAL(m.getValue(2, 1), 1);
		BOOST_CHECK_EQUAL(m.getValue(2, 2), 0);
	}
}

BOOST_AUTO_TEST_CASE_TEMPLATE( value_setting, MatrixType, matrix_type_list )
{
	if( mla::matrix::Traits<MatrixType>::is_writeable() )
	{
		MatrixType m(3,3);
		m.setZero();

		m.setValue( 1, 1, (typename MatrixType::scalar_type)1);

		BOOST_CHECK_EQUAL(m.getValue(1, 1), 1.0);
		BOOST_CHECK_EQUAL(m.getValue(0, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(0, 1), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(0, 2), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 1), 1.0);
		BOOST_CHECK_EQUAL(m.getValue(1, 2), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 0), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 1), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 2), 0.0);
	}
}

BOOST_AUTO_TEST_CASE_TEMPLATE( matrix_assigning, MatrixType, matrix_type_list )
{
	if( mla::matrix::Traits<MatrixType>::is_writeable() )
	{
		MatrixType m(3,3), n;
		m.setZero();

		m.setValue( 1, 1, 1.0);

		n = m;

		BOOST_CHECK_EQUAL(n.getValue(1, 1), 1.0);
		BOOST_CHECK_EQUAL(n.getValue(0, 0), 0.0);
		BOOST_CHECK_EQUAL(n.getValue(0, 1), 0.0);
		BOOST_CHECK_EQUAL(n.getValue(0, 2), 0.0);
		BOOST_CHECK_EQUAL(n.getValue(1, 0), 0.0);
		BOOST_CHECK_EQUAL(n.getValue(1, 1), 1.0);
		BOOST_CHECK_EQUAL(n.getValue(1, 2), 0.0);
		BOOST_CHECK_EQUAL(n.getValue(2, 0), 0.0);
		BOOST_CHECK_EQUAL(n.getValue(2, 1), 0.0);
		BOOST_CHECK_EQUAL(m.getValue(2, 2), 0.0);
	}
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: ruimaciel/mla

max_issues_repo_name: ruimaciel/mla

max_forks_repo_name: ruimaciel/mla


text: #include <boost/spirit/include/qi_match_attr.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include <iostream>
#include <fstream>
#include <cstring>

#include "Common.hpp"
#include "BUSData.h"

#include "bustools_predict.h"
#include <iostream>
#include <fstream>
#include <cmath>

#include <Eigen/Core>
#include <unsupported/Eigen/SpecialFunctions>
#include <LBFGSB.h>
#include <cppoptlib/meta.h>
#include <cppoptlib/boundedproblem.h>
#include <cppoptlib/solver/lbfgsbsolver.h>
#include <thread>
#include <memory>
#include <mutex>

//based on the boost implementation
//r is the size parameter
inline double DensityNegBin(double k, double r, double mu) {
	//so, the standard parameterization is r and p, but we use mean (mu) and size
	//mean = mu = r(1-p)/p <=> p = r/(mu+r), size (i.e. r) in this parameterization and r in the standard one is the same.
	double p = r / (mu + r);
	return exp(lgamma(r + k) - lgamma(r) - lgamma(k + 1)) * pow(p, r) * pow((1 - p), k);
}

//So, this one is different from the negative binomial below in that the zero is not included.
//We here want to look at the log likelihood for observed molecules (i.e. not including the zeros) and how they would fit the negative
//binomial, which means that the density function needs to be modified. The probability for zero needs to be set to zero, and the other
//values scaled up to get to a sum of 1.
inline double ZTNBLogLikelihood(const double* hist, size_t histLen, double size, double mu) {
	//first, get the probability of zero, probZero - all other values should be divided by 1-probZero to normalize the sum of non-zeros to 1
	double probZero = DensityNegBin(0, size, mu);
	double scale = 1 - probZero;

	//now calculate the log likelihood
	double ll = 0;
	for (size_t i = 0; i < histLen; ++i) {
		ll += log(DensityNegBin(double(i + 1), size, mu) / scale) * (*(hist + i));
	}

	return ll;
}

inline double NegBinLogLikelihood(const double* hist, size_t histLen, double size, double mu, double zeroCount) {
	//loglikelihood for nonzero counts
	double ll = 0;
	for (size_t i = 0; i < histLen; ++i) {
		ll += log(DensityNegBin(double(i + 1), size, mu)) * (*(hist + i));
	}
	//...and for the zero counts
	ll += log(DensityNegBin(0, size, mu)) * zeroCount;

	return ll;
}

//using LBFGSpp
class Optimizer
{
public:
	//make the variables public to be able to update them without a lot of extra code
	double mean = 0;
	double zeroCounts = 0;
	double totMol = 0;
	const double* hist = nullptr;
	size_t histLen = 0;
	Optimizer() {}

	//returns objective function value. The grad vector should be modified (i.e. set gradients)
	//the vector x is the in parameters, in this case the size (i.e. r) parameter in the distribution
	//The math here is the same as in preseqR
	double operator()(const Eigen::VectorXd& x, Eigen::VectorXd& grad) {
		double firstTerm = Eigen::numext::digamma(x[0]) * zeroCounts;
		for (size_t i = 0; i < histLen; ++i) {
			firstTerm += *(hist+i) * Eigen::numext::digamma(double(i + 1 + x[0]));
		}
		firstTerm /= totMol;

		double secondTerm = Eigen::numext::digamma(x[0]);

		double thirdTerm = log(x[0]) - log(x[0] + mean);

		grad[0] = -firstTerm + secondTerm - thirdTerm;

		return -NegBinLogLikelihood(hist, histLen, x[0], mean, zeroCounts) / totMol;
	}
};

//using CppOptimizationLibrary
class Optimizer2 : public cppoptlib::BoundedProblem<double, 1> {
public:
	//make the variables public to be able to update them without a lot of extra code
	double mean = 0;
	double zeroCounts = 0;
	double totMol = 0;
	const double* hist = nullptr;
	size_t histLen = 0;

	//The math here is the same as in preseqR

	double value(const TVector& x) {
		return -NegBinLogLikelihood(hist, histLen, x[0], mean, zeroCounts) / totMol;
	}

	void gradient(const TVector& x, TVector& grad)
	{
		double firstTerm = Eigen::numext::digamma(x[0]) * zeroCounts;
		for (size_t i = 0; i < histLen; ++i) {
			firstTerm += *(hist + i) * Eigen::numext::digamma(double(i + 1 + x[0]));
		}
		firstTerm /= totMol;

		double secondTerm = Eigen::numext::digamma(x[0]);

		double thirdTerm = log(x[0]) - log(x[0] + mean);

		grad[0] = -firstTerm + secondTerm - thirdTerm;
	}
};

//size and mu are both in and out parameters - the in parameters are the inital values in the expectation maximization algorithm
//the hist is just a pointer into the large hist vector to avoid unnecessary copying
//returns the log likelihood
//This function uses LBFGSpp
double PredictZTNBEmAlg1(const double* hist, size_t histLen, double& size, double& mu) {
	//estimate start values
	double histSum = 0; //S in the R code
	double histWeights = 0;
	for (size_t i = 0; i < histLen; ++i) {
		histSum += *(hist + i);
		histWeights += (*(hist + i)) * (i+1);
	}
	mu = (histWeights - histSum)/histSum;
	size = 1;

	double zeroProbability = DensityNegBin(0, size, mu);

	//estimate the total number of molecules
	double totMol = histSum / (1 - zeroProbability); //L in the R code

	//the estimated number of zero counts
	double zeroCounts = totMol * zeroProbability;

	//estimate the mean and variance
	double mean = 0, variance = 0;
	double precalcMeanSum = 0; //reuse this sum in the EM loop below since it is constant throughout the algorithm
	for (size_t i = 0; i < histLen; ++i) {
		precalcMeanSum += *(hist + i) * double(i + 1);//remember that the hist index is off by one, i.e. that the hist[0] corresponds to molecules with 1 copy and so forth
	}
	mean = precalcMeanSum / totMol;
	//..and the variance
	for (size_t i = 0; i < histLen; ++i) {
		variance += *(hist + i) * pow(double(i + 1) - mean, 2);
	}
	variance = (variance + mean * mean * zeroCounts) / (totMol - 1);

	Optimizer op;

	// set all params
	op.mean = mean;
	op.zeroCounts = zeroCounts;
	op.totMol = totMol;
	op.hist = hist;
	op.histLen = histLen;

	//using LBFGSpp
	LBFGSpp::LBFGSBParam<double> param;
	//So, the epsilon and max iterations here are for the LBFGS algorithm, not the EM, that comes later!
	param.epsilon = 1e-8; //It seems that the optim function in R uses about 1e-8
	param.max_iterations = 100; //defaults to 100 in the optim function in R
	param.max_linesearch = 500; //don't want it to fail...

	// Create solver and function object
	LBFGSpp::LBFGSBSolver<double> solver(param);

	// Initial guess
	Eigen::VectorXd x = Eigen::VectorXd::Zero(1);

	if (variance > mean) {
		x[0] = mean * mean / (variance - mean);
	}
	else {
		x[0] = size;
	}

	// Bounds
	Eigen::VectorXd lb = Eigen::VectorXd::Constant(1, 0.0001);
	Eigen::VectorXd ub = Eigen::VectorXd::Constant(1, 10000);

	// x will be overwritten to be the best point found
	double fx = 0;
	int numIter = solver.minimize(op, x, fx, lb, ub);

	size_t iter = 0;
	double lastNegLL = 10000000000000;
	double currNegLL = -ZTNBLogLikelihood(hist, histLen, x[0], mean);

	//termination criteria for the EM algorithm
	const double MAX_ERROR = 1e-8;
	const double MAX_ERROR_FAST = 1e-5;
	const size_t MAX_ITER = 100000;
	const size_t ITER_FAST_LIMIT = 200; //if the number of iterations goes over this number, start using the lower error threshold, MAX_ERROR_FAST, to quicken things up

	//The EM algorithm starts here
	while (fabs(lastNegLL - currNegLL) / histSum > MAX_ERROR&&
		iter < MAX_ITER &&
		!(iter >= ITER_FAST_LIMIT && fabs(lastNegLL - currNegLL) / histSum <= MAX_ERROR_FAST))
	{
		lastNegLL = currNegLL;

		//update distribution params
		size = x[0];
		mu = mean;

		// E step: estimate the number of unobserved species

		zeroProbability = DensityNegBin(0, size, mu);
		totMol = histSum / (1 - zeroProbability);
		zeroCounts = totMol * zeroProbability;

		//mean and variance
		mean = precalcMeanSum / totMol;
		//...and variance
		for (size_t i = 0; i < histLen; ++i) {
			variance += *(hist + i) * pow(double(i + 1) - mean, 2);
		}
		variance = (variance + mean * mean * zeroCounts) / (totMol - 1);

		// M step: estimate the parameters size and mu
		//rerun the LBFGSB

		// set params
		op.mean = mean;
		op.zeroCounts = zeroCounts;
		op.totMol = totMol;

		if (variance > mean) {
			x[0] = mean * mean / (variance - mean);
		}
		else {
			x[0] = size;
		}

		//avoid printing of warning text in console
		if (x[0] > 10000.0) {
			x[0] = 10000.0;
		}
		if (x[0] < 0.0001) {
			x[0] = 0.0001;
		}

		int numIter = solver.minimize(op, x, fx, lb, ub);

		currNegLL = -ZTNBLogLikelihood(hist, histLen, x[0], mean);

		//std::cout << "error: " << (lastNegLL - currNegLL) / histSum << "\n";

		++iter;

		//if (iter > 300) break;//tmp
	}

	//std::cout << "iterations: " << iter << " error: " << (lastNegLL - currNegLL) / histSum << "\n";

	return -currNegLL;
}

//Similar to above, but uses CppOptimizationLibrary
double PredictZTNBEmAlg2(const double* hist, size_t histLen, double& size, double& mu) {

	//estimate start values
	double histSum = 0; //S in the R code
	double histWeights = 0;
	for (size_t i = 0; i < histLen; ++i) {
		histSum += *(hist + i);
		histWeights += (*(hist + i)) * (i+1);
	}
	mu = (histWeights - histSum)/histSum;
	size = 1;

	double zeroProbability = DensityNegBin(0, size, mu);

	//estimate the total number of molecules
	double totMol = histSum / (1 - zeroProbability); //L in the R code

	//the estimated number of zero counts
	double zeroCounts = totMol * zeroProbability;

	//estimate the mean and variance
	double mean = 0, variance = 0;
	double precalcMeanSum = 0; //reuse this sum in the EM loop below since it is constant throughout the algorithm
	for (size_t i = 0; i < histLen; ++i) {
		precalcMeanSum += *(hist + i) * double(i + 1);//remember that the hist index is off by one, i.e. that the hist[0] corresponds to molecules with 1 copy and so forth
	}
	mean = precalcMeanSum / totMol;
	//..and the variance
	for (size_t i = 0; i < histLen; ++i) {
		variance += *(hist + i) * pow(double(i + 1) - mean, 2);
	}
	variance = (variance + mean * mean * zeroCounts) / (totMol - 1);

	Optimizer2 op;
	// set all params
	op.mean = mean;
	op.zeroCounts = zeroCounts;
	op.totMol = totMol;
	op.hist = hist;
	op.histLen = histLen;

	//using CppOptimizationLibrary
	op.setLowerBound(Optimizer2::TVector::Ones(1) * 0.0001);
	op.setUpperBound(Optimizer2::TVector::Ones(1) * 10000);
	cppoptlib::LbfgsbSolver<Optimizer2> solver;
	//change the stop criteria to speed things up
	cppoptlib::LbfgsbSolver<Optimizer2>::TCriteria crit = cppoptlib::LbfgsbSolver<Optimizer2>::TCriteria::defaults();
	crit.iterations = 1000; //default is 10000
	solver.setStopCriteria(crit);
	Optimizer2::TVector x = Optimizer2::TVector::Zero(1);

	if (variance > mean) {
		x[0] = mean * mean / (variance - mean);
	}
	else {
		x[0] = size;
	}

//	solver.setDebug(cppoptlib::DebugLevel::Low);
	// x will be overwritten to be the best point found
	solver.minimize(op, x);

	// Bounds
	Eigen::VectorXd lb = Eigen::VectorXd::Constant(1, 0.0001);
	Eigen::VectorXd ub = Eigen::VectorXd::Constant(1, 10000);

	size_t iter = 0;
	double lastNegLL = 10000000000000;
	double currNegLL = -ZTNBLogLikelihood(hist, histLen, x[0], mean);

	//termination criteria for the EM algorithm
	const double MAX_ERROR = 1e-8;
	const double MAX_ERROR_FAST = 1e-5;
	const size_t MAX_ITER = 400; //Since this algorithm is pretty slow, don't iterate too many times. The algorithm may get stuck here, on a few values only, don't let them dictate the total execution time too much.
	const size_t ITER_FAST_LIMIT = 200; //if the number of iterations goes over this number, start using the lower error threshold, MAX_ERROR_FAST, to quicken things up
	const double LARGE_NEG_LL = 100000000000000000.0;

	double bestNegLL = LARGE_NEG_LL;
	double bestMu = -1;
	double bestSize = -1;

	//The EM algorithm starts here
	while (fabs(lastNegLL - currNegLL) / histSum > MAX_ERROR&&
		iter < MAX_ITER &&
		!(iter >= ITER_FAST_LIMIT && fabs(lastNegLL - currNegLL) / histSum <= MAX_ERROR_FAST))
	{
		//temp, remove
		//double err = fabs(lastNegLL - currNegLL);

		lastNegLL = currNegLL;

		//update distribution params
		size = x[0];
		mu = mean;

		// E step: estimate the number of unobserved species

		zeroProbability = DensityNegBin(0, size, mu);
		totMol = histSum / (1 - zeroProbability);
		zeroCounts = totMol * zeroProbability;

		//mean and variance
		mean = precalcMeanSum / totMol;
		//...and variance
		for (size_t i = 0; i < histLen; ++i) {
			variance += *(hist + i) * pow(double(i + 1) - mean, 2);
		}
		variance = (variance + mean * mean * zeroCounts) / (totMol - 1);

		// M step: estimate the parameters size and mu
		//rerun the LBFGSB

		// set params
		op.mean = mean;
		op.zeroCounts = zeroCounts;
		op.totMol = totMol;

		if (variance > mean) {
			x[0] = mean * mean / (variance - mean);
		}
		else {
			x[0] = size;
		}

		//avoid printing of warning text in console
		if (x[0] > 10000.0) {
			x[0] = 10000.0;
		}
		if (x[0] < 0.0001) {
			x[0] = 0.0001;
		}

		solver.minimize(op, x);
		//std::cout << "val: " << x[0] << "\n";

		currNegLL = -ZTNBLogLikelihood(hist, histLen, x[0], mean);

		//std::cout << "Iteration: " << iter << " x val: " << x[0] << " ll " << currNegLL << " lldiff: " << err << "\n";

		//keep track of the best values we had in case it doesn't converge - better to return those
		if (currNegLL < bestNegLL) {
			bestNegLL = currNegLL;
			bestMu = mu;
			bestSize = size;
		}

		//std::cout << "error: " << (lastNegLL - currNegLL) / histSum << "\n";

		++iter;

		//if (iter > 300) break;//tmp
	}

	//std::cout << "iterations: " << iter << " error: " << (lastNegLL - currNegLL) / histSum << "\n";

	if (bestNegLL < LARGE_NEG_LL) {
		currNegLL = bestNegLL;
		mu = bestMu;
		size = bestSize;
	}

	//std::cout << "Iteration: " << iter << " mu: " << mu << " size: " << size << " ll: " << currNegLL << "\n";

	return -currNegLL;
}

//size and mu are out parameters describing the negative binomial
double PredictZTNBForGene(const double* hist, size_t histLen, double t, double& size, double& mu, int index) {
	double histSum = 0; //S in the R code
	for (size_t i = 0; i < histLen; ++i) {
		histSum += *(hist + i);
	}
	if (histSum == 0.0) {
		return 0.0;//nothing to do...
	}
	//initial values of negative binomial, taken from R code
	size = 1.0;
	mu = 0.5;
	//fit the ZTNB (will update size and mu)
	//So, the trick here is to first use Alg1 - it is faster, but fails sometimes. If it fails,
	//use Alg2
	try {
		//std::cout << "Alg1: " << index << "\n";
		PredictZTNBEmAlg1(hist, histLen, size, mu);
		//std::cout << "Alg1 done: " << index << "\n";
	}
	catch (std::exception&)
	{
		//std::cout << "Alg2: " << index << "\n";
		PredictZTNBEmAlg2(hist, histLen, size, mu);
		//std::cout << "Alg2 done: " << index << "\n";
	}

	//std::cout << "Mu: " << mu << " Size: " << size << "\n";

	//estimate the total number of molecules
	double zeroProbability = DensityNegBin(0, size, mu);
	double totMol = histSum / (1 - zeroProbability); //L in the R code

	//The prediction is based on the assumption that the size parameter remains the same, while
	//the mean is scaled up with t. We then look at how many non-zero molecules we would get
	double zeroProbScaled = DensityNegBin(0, size, mu*t);
	double result = totMol * (1 - zeroProbScaled);
	if (result < histSum) { //safety check, we should never get fewer molecules than what we started with!
		result = histSum;
	}

	return result;
}

//This is a thread scheduler
class PredictionExecuter
{
public:
	PredictionExecuter(const std::vector<double>& hists, const std::vector<size_t>& histLengths, double t, std::vector<double>& predVals, std::vector<double>& sizeVals, std::vector<double>& muVals, const uint32_t histmax)
		: m_hists(hists)
		, m_histLengths(histLengths)
		, m_t(t)
		, m_predVals(predVals)
		, m_sizeVals(sizeVals)
		, m_muVals(muVals)
		, m_histmax(histmax)
	{}
	void Execute(int numThreads) {
		//Some test code:
		//std::cout << "Predicting 1345 \n";
		//double size = 0;
		//double mu = 0;
		//PredictZTNBForGene(&m_hists[1345 * m_histmax], m_histLengths[1345], m_t, size, mu, 1345);
		//std::cout << "End Predicting 1345 \n";

		//create the threads
		for (int i = 0; i < numThreads; ++i) {
			m_threads.push_back(std::shared_ptr<std::thread> (new std::thread(&PredictionExecuter::ThreadFunc, this)));
		}
		//wait for them to finish
		for (int i = 0; i < numThreads; ++i) {
			m_threads[i]->join();
		}
		std::cerr << "\n";
	}
private:
	void ThreadFunc() {
		size_t i = 0;
		while (GetNextIndex(i)) {
			double size = 0;
			double mu = 0;
			m_predVals[i] = PredictZTNBForGene(&m_hists[i * m_histmax], m_histLengths[i], m_t, size, mu, i);
			m_sizeVals[i] = size;
			m_muVals[i] = mu;
		}
	}
	bool GetNextIndex(size_t& index) {
		std::lock_guard<std::mutex> lg(m_mutex);
		if (m_nextIndex == m_histLengths.size()) {
			return false;
		} else {
			std::cerr << "\rProcessing gene: " << m_nextIndex + 1 << " of " << m_histLengths.size();//so, we use indexes starting at 1 in the printout
			index = m_nextIndex++;
			return true;
		}
	}
	std::mutex m_mutex;
	size_t m_nextIndex = 0;
	const std::vector<double>& m_hists;
	const std::vector<size_t>& m_histLengths;
	double m_t;
	std::vector<double>& m_predVals;
	std::vector<double>& m_sizeVals;
	std::vector<double>& m_muVals;
	const uint32_t m_histmax;
	std::vector<std::shared_ptr<std::thread>> m_threads;
};

void bustools_predict(Bustools_opt &opt) {

	//Prepare and load histograms
	//////////////////

	//read the hist file
	std::string hist_ifn = opt.predict_input + ".hist.txt";
	std::string counts_ifn = opt.predict_input + ".mtx";
	std::string gene_ifn = opt.predict_input + ".genes.txt";
	std::string barcode_ifn = opt.predict_input + ".barcodes.txt";

	std::string corr_counts_ofn = opt.output + ".mtx";
	std::string nb_params_ofn = opt.output + ".nb_params.txt";
	std::string corr_gene_ofn = opt.output + ".genes.txt";
	std::string corr_barcode_ofn = opt.output + ".barcodes.txt";

	//just copy the genes and barcodes files, they don't change
	copy_file(gene_ifn, corr_gene_ofn);
	copy_file(barcode_ifn, corr_barcode_ofn);

	//Get gene list associated with count matrix and histograms
	std::vector<std::string> genes;
	parseGenesList(gene_ifn, genes);

	//Allocate histograms
	//Indexed as gene*histmax + histIndex
	size_t n_genes = genes.size();
	const uint32_t histmax = 100;//set molecules with more than histmax copies to histmax
	std::vector<double> histograms = std::vector<double>(n_genes * histmax, 0);
	std::vector<size_t> histogramLengths = std::vector<size_t>(n_genes, 0);
	std::string line;

	//load histograms file
	{
		std::ifstream inf(hist_ifn);
		size_t geneIndex = 0;
		while (std::getline(inf, line)) {
			std::stringstream ss(line);
			double num = 0;
			size_t histIndex = 0;
			while (ss >> num) {
				histograms[geneIndex * histmax + histIndex++] = num;
			}
			histogramLengths[geneIndex] = histIndex;
			++geneIndex;
		}
	}

	//fix the histograms if they have only a single non-zero value (i.e. for example looks like this: 1 0 0 0)
	//however, do not fix completely empty histograms
	//the strategy is to always add a one after the last item in the histogram
	//(this will be next to the only non-zero number, which will be last)
	//the strategy may lead to that a gene gets more counts. This, however, doesn't matter since we only calculate
	//a scaling factor for the gene, so that will be normalized out.
	for (size_t i = 0; i < histogramLengths.size(); ++i) {
		size_t nonZeros = 0;
		for (size_t j = 0; j < histogramLengths[i]; ++j) {
			if (histograms[i * histmax + j] != 0) {
				++nonZeros;
				if (nonZeros > 1) {
					break;
				}
			}
		}
		if (nonZeros == 1) {
			if (histogramLengths[i] == histmax) { //rare case, add a count before instead of after
				histograms[i * histmax + histmax - 2] = 1;
			} else {
				histograms[i * histmax + histogramLengths[i]] = 1;
				histogramLengths[i]++;//also extend the histogram one step;
			}
			//so, if we can, remove one count from the current value
			if (histograms[i * histmax + histogramLengths[i] - 1] > 1) {
				histograms[i * histmax + histogramLengths[i] - 1] -= 1.0;
			}
		}
	}

	//Predict
	//////////////////

	std::vector<double> predVals(n_genes, 0);
	std::vector<double> sizeVals(n_genes, 0);
	std::vector<double> muVals(n_genes, 0);
	int numThreads = std::thread::hardware_concurrency();
	std::cerr << "Using " << numThreads << " threads\n";
	PredictionExecuter pe(histograms, histogramLengths, opt.predict_t, predVals, sizeVals, muVals, histmax);
	pe.Execute(numThreads);

	//calculate the sum of all histograms and all predvals:
	double histSum = 0;
	std::vector<double> umisPerGene(n_genes, 0);
	std::vector<double> countsPerGene(n_genes, 0);
	double predSum = 0;

	for (size_t i = 0; i < predVals.size(); ++i) {
		for (size_t j = 0; j < histogramLengths[i]; ++j) {
			umisPerGene[i] += histograms[i * histmax + j];
			countsPerGene[i] += histograms[i * histmax + j]*(j+1);
		}
		histSum += umisPerGene[i];
		predSum += predVals[i];
	}

	//so the genes should be scaled according to the following:
	//geneScaling = histSum/predSum * predVal/umisPerGene;
	std::vector<double> geneScaling(n_genes, 1.0);//1.0 for all empty genes
	double globScale = histSum / predSum;
	for (size_t i = 0; i < predVals.size(); ++i) {
		if (umisPerGene[i] > 0) {
			geneScaling[i] = globScale * predVals[i] / umisPerGene[i];
		}
	}

	//Modify counts
	//////////////////

	std::cerr << "Creating corrected counts matrix...\n";

	//read and write the counts matrix
	{
		std::ifstream inf(counts_ifn);
		std::string test;
		std::ofstream of(corr_counts_ofn);
		//first number is cells, second genes
		//first handle header, we just leave that untouched
		while (std::getline(inf, line)) {
			of << line << '\n';
			if ((!line.empty()) && line[0] != '%') {
				break;
			}
		}
		//now all data with scaling
		size_t cell = 0, gene = 0;
		double count = 0;
		while (std::getline(inf, line)) {
			if ((!line.empty()) && line[0] == '%') {
				of << line << '\n';
				continue;
			}
			std::stringstream ss(line);
			if (ss >> cell >> gene >> count) {
				of << cell << " " << gene << " " << count * geneScaling[gene-1] << '\n';
			}
		}
	}

	//write negative binomial params (file with header)
	std::cerr << "Writing negative binomial params...\n";
	{
		std::ofstream of(nb_params_ofn);

		//header
		of << "gene\tmu\tsize\tUMIs\tcounts\n";

		for (size_t i = 0; i < genes.size(); ++i) {
			of << genes[i] << '\t' << muVals[i] << '\t' << sizeVals[i] << '\t' << umisPerGene[i] << '\t' << countsPerGene[i] << '\n';
		}
	}

}

max_stars_repo_name: Yenaled/bustools

max_issues_repo_name: Yenaled/bustools

max_forks_repo_name: Yenaled/bustools


text:
/*
 * GridTools
 *
 * Copyright (c) 2014-2019, ETH Zurich
 * All rights reserved.
 *
 * Please, refer to the LICENSE file in the root directory.
 * SPDX-License-Identifier: BSD-3-Clause
 */

#pragma once

#include <type_traits>

#include <boost/preprocessor.hpp>

#include "defs.hpp"

// internal
#define GT_META_INTERNAL_APPLY(fun, args) BOOST_PP_REMOVE_PARENS(fun)<BOOST_PP_REMOVE_PARENS(args)>

#if GT_BROKEN_TEMPLATE_ALIASES

/**
 * backward compatible way to call function
 */
#define GT_META_CALL(fun, args) typename GT_META_INTERNAL_APPLY(fun, args)::type

/**
 * backward compatible way to define an alias to the function composition
 */
#define GT_META_DEFINE_ALIAS(name, fun, args) \
    struct name : GT_META_INTERNAL_APPLY(fun, args) {}

#define GT_META_LAZY_NAMESPACE inline namespace lazy

#define GT_META_DELEGATE_TO_LAZY(fun, signature, args) static_assert(1, "")

#else

/**
 * backward compatible way to call function
 */
#define GT_META_CALL(fun, args) GT_META_INTERNAL_APPLY(fun, args)
/**
 * backward compatible way to define an alias to the function composition
 */
#define GT_META_DEFINE_ALIAS(name, fun, args) using name = GT_META_INTERNAL_APPLY(fun, args)

#define GT_META_LAZY_NAMESPACE namespace lazy

#define GT_META_DELEGATE_TO_LAZY(fun, signature, args) \
    template <BOOST_PP_REMOVE_PARENS(signature)>       \
    using fun = typename lazy::fun<BOOST_PP_REMOVE_PARENS(args)>::type

#endif

/**
 *  NVCC bug workaround: sizeof... works incorrectly within template alias context.
 */
#ifdef __NVCC__

namespace gridtools {
    namespace meta {
        template <class... Ts>
        struct sizeof_3_dots : std::integral_constant<std::size_t, sizeof...(Ts)> {};
    } // namespace meta
} // namespace gridtools

#define GT_SIZEOF_3_DOTS(Ts) ::gridtools::meta::sizeof_3_dots<Ts...>::value
#else
#define GT_SIZEOF_3_DOTS(Ts) sizeof...(Ts)
#endif

max_stars_repo_name: mbianco/gridtools

max_issues_repo_name: mbianco/gridtools

max_forks_repo_name: mbianco/gridtools