Working on item_value

CMakeLists.txt

@@ -32,7 +32,7 @@ endif()
 # set the project name
 project(
 	libcifpp
-	VERSION 9.0.1
+	VERSION 9.0.5
 	LANGUAGES CXX C)
 
 list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
@@ -43,6 +43,7 @@ include(GenerateExportHeader)
 include(CTest)
 include(ExternalProject)
 include(FetchContent)
+include(VersionString)
 
 # When building with ninja-multiconfig, build both debug and release by default
 if(CMAKE_GENERATOR STREQUAL "Ninja Multi-Config")
@@ -168,14 +169,19 @@ if(MSVC)
 	endforeach()
 endif()
 
-# First check if <format> is available
-find_file(FMT NAME format)
+set(CMAKE_CXX_STANDARD 20)
 
-if(FMT EQUAL "FMT-NOTFOUND")
-	if(NOT (fmt_FOUND OR TARGET fmt))
-		find_package(fmt REQUIRED)
-		message(FATAL_ERROR "fmt not found, compiler too old, you're out of luck")
-	endif()
+# Using fast_float for float parsing, but only if needed
+try_compile(STD_CHARCONV_COMPILING
+	SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/cmake/test-charconv.cpp)
+
+if(NOT STD_CHARCONV_COMPILING)
+	message(NOTICE "libcifpp: Using fast_float for std::from_chars")
+	FetchContent_Declare(fastfloat
+		GIT_REPOSITORY "https://github.com/fastfloat/fast_float"
+		GIT_TAG v8.0.2
+		EXCLUDE_FROM_ALL)
+	FetchContent_MakeAvailable(fastfloat)
 endif()
 
 find_package(Threads)
@@ -189,10 +195,6 @@ include(FindPkgConfig)
 
 if(PKG_CONFIG_FOUND)
 	pkg_check_modules(PCRE2 IMPORTED_TARGET libpcre2-8)
-
-	if(PCRE2_FOUND)
-		message(STATUS "Using pcre2 found using pkg-config")
-	endif()
 endif()
 
 if(NOT PCRE2_FOUND)
@@ -220,10 +222,6 @@ else()
 	set(EIGEN_INCLUDE_DIR ${SOURCE_DIR})
 endif()
 
-# Create a revision file, containing the current git version info
-include(VersionString)
-write_version_header(${CMAKE_CURRENT_SOURCE_DIR}/src/ LIB_NAME "LibCIFPP")
-
 # SymOp data table
 if(CIFPP_RECREATE_SYMOP_DATA)
 	# The tool to create the table
@@ -245,6 +243,9 @@ if(CIFPP_RECREATE_SYMOP_DATA)
 		"$ENV{CLIBD}/symop.lib")
 endif()
 
+# Create a revision file, containing the current git version info
+write_version_header("${CMAKE_CURRENT_SOURCE_DIR}/src/" LIB_NAME "LibCIFPP")
+
 # Sources
 set(project_sources
 	${CMAKE_CURRENT_SOURCE_DIR}/src/category.cpp
@@ -258,18 +259,18 @@ set(project_sources
 	${CMAKE_CURRENT_SOURCE_DIR}/src/validate.cpp
 	${CMAKE_CURRENT_SOURCE_DIR}/src/text.cpp
 	${CMAKE_CURRENT_SOURCE_DIR}/src/utilities.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/atom_type.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/compound.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/point.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/symmetry.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/model.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/cif2pdb.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb_record.hpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.hpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/reconstruct.cpp
-	${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/validate-pdbx.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/atom_type.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/compound.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/point.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/symmetry.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/model.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/cif2pdb.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb_record.hpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.hpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/pdb2cif_remark_3.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/reconstruct.cpp
+	# ${CMAKE_CURRENT_SOURCE_DIR}/src/pdb/validate-pdbx.cpp
 )
 
 set(project_headers
@@ -347,8 +348,8 @@ else()
 	target_link_libraries(cifpp PRIVATE $<BUILD_INTERFACE:pcre2s>)
 endif()
 
-if(fmt_FOUND)
-	target_link_libraries(cifpp PUBLIC fmt)
+if(NOT STD_CHARCONV_COMPILING)
+	target_link_libraries(cifpp PRIVATE FastFloat::fast_float)
 endif()
 
 if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
@@ -487,7 +488,6 @@ if(CIFPP_DATA_DIR AND CIFPP_DOWNLOAD_CCD)
 endif()
 
 set(CONFIG_TEMPLATE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/cmake/cifpp-config.cmake.in)
-set(REQUIRE_FMT ${fmt_FOUND})
 
 configure_package_config_file(
 	${CONFIG_TEMPLATE_FILE} ${CMAKE_CURRENT_BINARY_DIR}/cifpp/cifpp-config.cmake

changelog

@@ -1,3 +1,19 @@
+Version 9.0.5
+- Added exists to compound_factory
+- Added sub_matrix, fix and extend determinant calculation
+- Added yet another structure::create_non_poly
+- Remove revision.hpp file in make clean (new VersionString.cmake)
+
+Version 9.0.4
+- Fix various stopping and reconstruction errors
+
+Version 9.0.3
+- Reconstruction fixed when some entity ids are missing
+
+Version 9.0.2
+- Fix code that reconstructs sequences, could throw a map::at
+- Many optimisations in validation and reconstruction code.
+
 Version 9.0.1
 - Use pcre2 from pkg-config if available, if not
   build a version from the original code.

cmake/VersionString.cmake

@@ -238,7 +238,7 @@ function(write_version_header dir)
 		if(res EQUAL 0)
 			set(REVISION_STRING "${out}")
 		else()
-			message(STATUS "Git hash not found, does this project has a 'build' tag?")
+			message(STATUS "Git hash not found, does this project have a 'build' tag?")
 		endif()
 	else()
 		message(STATUS "Git hash not found")

cmake/cifpp-config.cmake.in

@@ -8,8 +8,5 @@ include(CMakeFindDependencyMacro)
 
 find_dependency(Threads)
 find_dependency(ZLIB REQUIRED)
-if(@REQUIRE_FMT@)
-    find_dependency(fmt REQUIRED)
-endif()
 
 check_required_components(cifpp)

cmake/test-charconv.cpp

@@ -0,0 +1,17 @@
+#include <charconv>
+#include <cassert>
+#include <cstring>
+
+int main()
+{
+	float v;
+	char s[] = "1.0";
+
+	auto r = std::from_chars(s, s + strlen(s), v);
+
+	assert(r.ec == std::errc{});
+	assert(r.ptr = s + strlen(s));
+	assert(v == 1.0f);
+
+	return 0;
+}

include/cif++/category.hpp

@@ -28,12 +28,14 @@
 
 #include "cif++/forward_decl.hpp"
 
-#include "cif++/condition.hpp"
+// #include "cif++/condition.hpp"
+#include "cif++/item.hpp"
 #include "cif++/iterator.hpp"
 #include "cif++/row.hpp"
 #include "cif++/text.hpp"
 
-#include <array>
+// #include <array>
+#include <functional>
 
 /** \file category.hpp
  * Documentation for the cif::category class
@@ -108,16 +110,6 @@ class multiple_results_error : public std::runtime_error
 	}
 };
 
-// --------------------------------------------------------------------
-// These should be moved elsewhere, one day.
-
-/// \cond
-template <typename _Tp>
-inline constexpr bool is_optional_v = false;
-template <typename _Tp>
-inline constexpr bool is_optional_v<std::optional<_Tp>> = true;
-/// \endcond
-
 // --------------------------------------------------------------------
 
 /// The class category is a sequence container for rows of data values.
@@ -181,12 +173,8 @@ class category
 
 	const std::string &name() const { return m_name; } ///< Returns the name of the category
 
-	[[deprecated("use key_items instead")]] iset key_fields() const; ///< Returns the cif::iset of key item names. Retrieved from the @ref category_validator for this category
-
 	iset key_items() const; ///< Returns the cif::iset of key item names. Retrieved from the @ref category_validator for this category
 
-	[[deprecated("use key_item_indices instead")]] std::set<uint16_t> key_field_indices() const; ///< Returns a set of indices for the key items.
-
 	std::set<uint16_t> key_item_indices() const; ///< Returns a set of indices for the key items.
 
 	/// @brief Set the validator for this category to @a v
@@ -336,7 +324,7 @@ class category
 	struct key_element_type
 	{
 		std::string name;         ///< Name of the item
-		std::string value;        ///< Value to be found
+		item_value value;         ///< Value to be found
 		bool may_be_null = false; ///< If true, value should be same or empty
 	};
 
@@ -961,7 +949,7 @@ class category
 			for (auto i = b; i != e; ++i)
 			{
 				// item_value *new_item = this->create_item(*i);
-				r->append(add_item(i->name()), { i->value() });
+				r->append(add_item(i->name()), i->value());
 			}
 		}
 		catch (...)
@@ -1005,7 +993,7 @@ class category
 
 	// --------------------------------------------------------------------
 
-	using value_provider_type = std::function<std::string_view(std::string_view)>;
+	using value_provider_type = std::function<item_value(const item_value &)>;
 
 	/// \brief Update a single item named @a item_name in the rows that match
 	/// \a cond to values provided by a callback function \a value_provider
@@ -1036,7 +1024,7 @@ class category
 	/// That means, child categories are updated if the links are absolute
 	/// and unique. If they are not, the child category rows are split.
 
-	void update_value(condition &&cond, std::string_view item_name, std::string_view value)
+	void update_value(condition &&cond, std::string_view item_name, item_value value)
 	{
 		auto rs = find(std::move(cond));
 		std::vector<row_handle> rows;
@@ -1049,66 +1037,12 @@ class category
 	/// That means, child categories are updated if the links are absolute
 	/// and unique. If they are not, the child category rows are split.
 
-	void update_value(const std::vector<row_handle> &rows, std::string_view item_name, std::string_view value)
+	void update_value(const std::vector<row_handle> &rows, std::string_view item_name, item_value value)
 	{
-		update_value(rows, item_name, [value](std::string_view)
+		update_value(rows, item_name, [value](const item_value &v)
 			{ return value; });
 	}
 
-	// --------------------------------------------------------------------
-	// Naming used to be very inconsistent. For backward compatibility,
-	// the old function names are here as deprecated variants.
-
-	/// \brief Return the index number for \a column_name
-	[[deprecated("Use get_item_ix instead")]] uint16_t get_column_ix(std::string_view column_name) const
-	{
-		return get_item_ix(column_name);
-	}
-
-	/// @brief Return the name for column with index @a ix
-	/// @param ix The index number
-	/// @return The name of the column
-	[[deprecated("use get_item_name instead")]] std::string_view get_column_name(uint16_t ix) const
-	{
-		return get_item_name(ix);
-	}
-
-	/// @brief Make sure a item with name @a item_name is known and return its index number
-	/// @param item_name The name of the item
-	/// @return The index number of the item
-	[[deprecated("use add_item instead")]] uint16_t add_column(std::string_view item_name)
-	{
-		return add_item(item_name);
-	}
-
-	/** @brief Remove column name @a colum_name
-	 * @param column_name The column to be removed
-	 */
-	[[deprecated("use remove_item instead")]] void remove_column(std::string_view column_name)
-	{
-		remove_item(column_name);
-	}
-
-	/** @brief Rename column @a from_name to @a to_name */
-	[[deprecated("use rename_item instead")]] void rename_column(std::string_view from_name, std::string_view to_name)
-	{
-		rename_item(from_name, to_name);
-	}
-
-	/// @brief Return whether a column with name @a name exists in this category
-	/// @param name The name of the column
-	/// @return True if the column exists
-	[[deprecated("use has_item instead")]] bool has_column(std::string_view name) const
-	{
-		return has_item(name);
-	}
-
-	/// @brief Return the cif::iset of columns in this category
-	[[deprecated("use get_items instead")]] iset get_columns() const
-	{
-		return get_items();
-	}
-
 	// --------------------------------------------------------------------
 	/// \brief Return the index number for \a item_name
 
@@ -1117,7 +1051,7 @@ class category
 	/// @brief Return the name for item with index @a ix
 	/// @param ix The index number
 	/// @return The name of the item
-	std::string_view get_item_name(uint16_t ix) const
+	const std::string &get_item_name(uint16_t ix) const
 	{
 		if (ix >= m_items.size())
 			throw std::out_of_range("item index is out of range");
@@ -1200,7 +1134,7 @@ class category
 	}
 
   private:
-	void update_value(row *row, uint16_t item, std::string_view value, bool updateLinked, bool validate = true);
+	void update_value(row *row, uint16_t item, item_value value, bool updateLinked, bool validate = true);
 
 	void erase_orphans(condition &&cond, category &parent);
 

include/cif++/compound.hpp

@@ -179,7 +179,7 @@ class compound
 	friend class compound_factory_impl;
 	friend class local_compound_factory_impl;
 
-	compound(cif::datablock &db);
+	compound(datablock &db);
 
 	std::string m_id;
 	std::string m_name;
@@ -270,11 +270,15 @@ class compound_factory
 		return is_std_base(res_name) or is_std_peptide(res_name);
 	}
 
+	/// Return whether @a res_name is water
 	bool is_water(std::string_view res_name) const
 	{
 		return res_name == "HOH" or res_name == "H2O" or res_name == "WAT";
 	}
 
+	/// Return whether @a res_name already exists, without creating it.
+	bool exists(std::string_view res_name) const;
+
 	/// \brief Create the compound object for \a id
 	///
 	/// This will create the compound instance for \a id if it doesn't exist already.
@@ -331,14 +335,14 @@ class compound_factory
 class compound_source
 {
   public:
-	compound_source(const cif::file &file)
+	compound_source(const file &file)
 	{
-		cif::compound_factory::instance().push_dictionary(file);
+		compound_factory::instance().push_dictionary(file);
 	}
 
 	~compound_source()
 	{
-		cif::compound_factory::instance().pop_dictionary();
+		compound_factory::instance().pop_dictionary();
 	}
 };
 

include/cif++/condition.hpp

@@ -371,7 +371,7 @@ namespace detail
 	{
 		key_equals_condition_impl(item &&i)
 			: m_item_name(i.name())
-			, m_value(std::forward<item>(i).value())
+			, m_value(std::forward<item_value>(i.value()))
 		{
 		}
 
@@ -409,7 +409,7 @@ namespace detail
 		std::string m_item_name;
 		uint16_t m_item_ix = 0;
 		bool m_icase = false;
-		std::string m_value;
+		item_value m_value;
 		std::optional<row_handle> m_single_hit;
 	};
 
@@ -466,111 +466,11 @@ namespace detail
 
 		std::string m_item_name;
 		uint16_t m_item_ix = 0;
-		std::string m_value;
+		item_value &m_value;
 		bool m_icase = false;
 		std::optional<row_handle> m_single_hit;
 	};
 
-	struct key_equals_number_condition_impl : public condition_impl
-	{
-		key_equals_number_condition_impl(const std::string &name, double v)
-			: m_item_name(name)
-			, m_value(v)
-		{
-		}
-
-		condition_impl *prepare(const category &c) override;
-
-		bool test(row_handle r) const override
-		{
-			return m_single_hit.has_value() ? *m_single_hit == r : r[m_item_ix].compare(m_value) == 0;
-		}
-
-		void str(std::ostream &os) const override
-		{
-			os << m_item_name << " == " << m_value;
-		}
-
-		virtual std::optional<row_handle> single() const override
-		{
-			return m_single_hit;
-		}
-
-		virtual bool equals(const condition_impl *rhs) const override
-		{
-			if (typeid(*rhs) == typeid(key_equals_number_condition_impl))
-			{
-				auto ri = static_cast<const key_equals_number_condition_impl *>(rhs);
-				if (m_single_hit.has_value() or ri->m_single_hit.has_value())
-					return m_single_hit == ri->m_single_hit;
-				else
-					// watch out, both m_item_ix might be the same while item_names might be diffent (in case they both do not exist in the category)
-					return m_item_ix == ri->m_item_ix and m_value == ri->m_value and m_item_name == ri->m_item_name;
-			}
-			return this == rhs;
-		}
-
-		std::string m_item_name;
-		uint16_t m_item_ix = 0;
-		double m_value;
-		std::optional<row_handle> m_single_hit;
-	};
-
-	struct key_equals_number_or_empty_condition_impl : public condition_impl
-	{
-		key_equals_number_or_empty_condition_impl(key_equals_number_condition_impl *equals)
-			: m_item_name(equals->m_item_name)
-			, m_value(equals->m_value)
-			, m_single_hit(equals->m_single_hit)
-		{
-		}
-
-		condition_impl *prepare(const category &c) override
-		{
-			m_item_ix = get_item_ix(c, m_item_name);
-			return this;
-		}
-
-		bool test(row_handle r) const override
-		{
-			bool result = false;
-			if (m_single_hit.has_value())
-				result = *m_single_hit == r;
-			else
-				result = r[m_item_ix].empty() or r[m_item_ix].compare(m_value) == 0;
-			return result;
-		}
-
-		void str(std::ostream &os) const override
-		{
-			os << '(' << m_item_name << " == " << m_value << " OR " << m_item_name << " IS NULL)";
-		}
-
-		virtual std::optional<row_handle> single() const override
-		{
-			return m_single_hit;
-		}
-
-		virtual bool equals(const condition_impl *rhs) const override
-		{
-			if (typeid(*rhs) == typeid(key_equals_number_or_empty_condition_impl))
-			{
-				auto ri = static_cast<const key_equals_number_or_empty_condition_impl *>(rhs);
-				if (m_single_hit.has_value() or ri->m_single_hit.has_value())
-					return m_single_hit == ri->m_single_hit;
-				else
-					// watch out, both m_item_ix might be the same while item_names might be diffent (in case they both do not exist in the category)
-					return m_item_ix == ri->m_item_ix and m_value == ri->m_value and m_item_name == ri->m_item_name;
-			}
-			return this == rhs;
-		}
-
-		std::string m_item_name;
-		uint16_t m_item_ix = 0;
-		double m_value;
-		std::optional<row_handle> m_single_hit;
-	};
-
 	struct key_compare_condition_impl : public condition_impl
 	{
 		template <typename COMP>
@@ -622,7 +522,7 @@ namespace detail
 
 		bool test(row_handle r) const override
 		{
-			std::string_view txt = r[m_item_ix].text();
+			auto txt = r[m_item_ix].get<std::string>();
 			return std::regex_match(txt.begin(), txt.end(), mRx);
 		}
 
@@ -693,7 +593,7 @@ namespace detail
 			{
 				try
 				{
-					std::string_view txt = r[f].text();
+					auto txt = r[f].get<std::string>();
 					if (std::regex_match(txt.begin(), txt.end(), mRx))
 					{
 						result = true;
@@ -970,26 +870,6 @@ inline condition operator or(condition &&a, condition &&b)
 				return condition(new detail::key_equals_or_empty_condition_impl(ci));
 		}
 
-		if (typeid(*a.m_impl) == typeid(detail::key_equals_number_condition_impl) and
-			typeid(*b.m_impl) == typeid(detail::key_is_empty_condition_impl))
-		{
-			auto ci = static_cast<detail::key_equals_number_condition_impl *>(a.m_impl);
-			auto ce = static_cast<detail::key_is_empty_condition_impl *>(b.m_impl);
-
-			if (ci->m_item_name == ce->m_item_name)
-				return condition(new detail::key_equals_number_or_empty_condition_impl(ci));
-		}
-
-		if (typeid(*b.m_impl) == typeid(detail::key_equals_number_condition_impl) and
-			typeid(*a.m_impl) == typeid(detail::key_is_empty_condition_impl))
-		{
-			auto ci = static_cast<detail::key_equals_number_condition_impl *>(b.m_impl);
-			auto ce = static_cast<detail::key_is_empty_condition_impl *>(a.m_impl);
-
-			if (ci->m_item_name == ce->m_item_name)
-				return condition(new detail::key_equals_number_or_empty_condition_impl(ci));
-		}
-
 		return condition(new detail::or_condition_impl(std::move(a), std::move(b)));
 	}
 
@@ -1066,20 +946,10 @@ struct key
 template <typename T>
 concept Numeric = ((std::is_floating_point_v<T> or std::is_integral_v<T>) and not std::is_same_v<T, bool>);
 
-/**
- * @brief Operator to create an equals condition based on a key @a key and a numeric value @a v
- */
-template <Numeric T>
-condition operator==(const key &key, const T &v)
-{
-	// TODO: change key_equals_etc... to use std::variant<double,int64_t> or something
-	return condition(new detail::key_equals_number_condition_impl(key.m_item_name, static_cast<double>(v)));
-}
-
 /**
  * @brief Operator to create an equals condition based on a key @a key and a value @a value
  */
-inline condition operator==(const key &key, std::string_view value)
+inline condition operator==(const key &key, item_value value)
 {
 	if (not value.empty())
 		return condition(new detail::key_equals_condition_impl({ key.m_item_name, value }));
@@ -1087,16 +957,6 @@ inline condition operator==(const key &key, std::string_view value)
 		return condition(new detail::key_is_empty_condition_impl(key.m_item_name));
 }
 
-/**
- * @brief Operator to create an equals condition based on a key @a key and a value @a value
- */
-template <typename T>
-	requires std::is_same_v<T, bool>
-inline condition operator==(const key &key, T value)
-{
-	return condition(new detail::key_equals_condition_impl({ key.m_item_name, value ? "y" : "n" }));
-}
-
 /**
  * @brief Operator to create a not equals condition based on a key @a key and a value @a v
  */

include/cif++/datablock.hpp

@@ -43,7 +43,7 @@ namespace cif
 
 /**
  * @brief A datablock is a list of category objects with some additional features
- * 
+ *
  */
 
 class datablock : public std::list<category>
@@ -53,7 +53,7 @@ class datablock : public std::list<category>
 
 	/**
 	 * @brief Construct a new datablock object with name @a name
-	 * 
+	 *
 	 * @param name The name for the new datablock
 	 */
 	datablock(std::string_view name)
@@ -80,7 +80,7 @@ class datablock : public std::list<category>
 	{
 		std::swap(a.m_name, b.m_name);
 		std::swap(a.m_validator, b.m_validator);
-		std::swap(static_cast<std::list<category>&>(a), static_cast<std::list<category>&>(b));
+		std::swap(static_cast<std::list<category> &>(a), static_cast<std::list<category> &>(b));
 	}
 
 	// --------------------------------------------------------------------
@@ -92,7 +92,7 @@ class datablock : public std::list<category>
 
 	/**
 	 * @brief Set the name of this datablock to @a name
-	 * 
+	 *
 	 * @param name The new name
 	 */
 	void set_name(std::string_view name)
@@ -102,27 +102,27 @@ class datablock : public std::list<category>
 
 	/**
 	 * @brief Attempt to load the dictionary specified in audit_conform category
-	 * 
+	 *
 	 */
 	void load_dictionary();
 
 	/**
 	 * @brief Set the validator object to @a v
-	 * 
+	 *
 	 * @param v The new validator object, may be null
 	 */
 	void set_validator(const validator *v);
 
 	/**
 	 * @brief Get the validator object
-	 * 
+	 *
 	 * @return const validator* The validator or nullptr if there is none
 	 */
 	const validator *get_validator() const;
 
 	/**
 	 * @brief Validates the content of this datablock and all its content
-	 * 
+	 *
 	 * @return true If the content is valid
 	 * @return false If the content is not valid
 	 */
@@ -131,7 +131,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Validates all contained data for valid links between parents and children
 	 * as defined in the validator
-	 * 
+	 *
 	 * @return true If all links are valid
 	 * @return false If all links are not valid
 	 */
@@ -140,7 +140,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Strip removes all categories and items that are invalid according
 	 * to the assigned validator. Will also add a valid audit_conform block.
-	 * 
+	 *
 	 * @return true if the remaining datablock is valid
 	 */
 	bool strip();
@@ -150,7 +150,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return the category named @a name, will create a new and empty
 	 * category named @a name if it does not exist.
-	 * 
+	 *
 	 * @param name The name of the category to return
 	 * @return category& Reference to the named category
 	 */
@@ -159,7 +159,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return the const category named @a name, will return a reference
 	 * to a static empty category if it was not found.
-	 * 
+	 *
 	 * @param name The name of the category to return
 	 * @return category& Reference to the named category
 	 */
@@ -168,7 +168,7 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return a pointer to the category named @a name or nullptr if
 	 * it does not exist.
-	 * 
+	 *
 	 * @param name The name of the category
 	 * @return category* Pointer to the category found or nullptr
 	 */
@@ -177,13 +177,12 @@ class datablock : public std::list<category>
 	/**
 	 * @brief Return a pointer to the category named @a name or nullptr if
 	 * it does not exist.
-	 * 
+	 *
 	 * @param name The name of the category
 	 * @return category* Pointer to the category found or nullptr
 	 */
 	const category *get(std::string_view name) const;
 
-
 	/**
 	 * @brief Return true if this datablock contains a non-empty category
 	 */
@@ -197,7 +196,7 @@ class datablock : public std::list<category>
 	 * new one if it is not found. The result is a tuple of an iterator
 	 * pointing to the category and a boolean indicating whether the category
 	 * was created or not.
-	 * 
+	 *
 	 * @param name The name for the category
 	 * @return std::tuple<iterator, bool> A tuple containing an iterator pointing
 	 * at the category and a boolean indicating whether the category was newly

include/cif++/iterator.hpp

@@ -26,6 +26,7 @@
 
 #pragma once
 
+#include "cif++/condition.hpp"
 #include "cif++/row.hpp"
 
 #include <array>
@@ -179,7 +180,7 @@ class iterator_impl
 	template <std::size_t... Is>
 	tuple_type get(std::index_sequence<Is...>) const
 	{
-		return m_current ? tuple_type{ m_current[m_item_ix[Is]].template as<Ts>()... } : tuple_type{};
+		return m_current ? tuple_type{ m_current[m_item_ix[Is]].template get<Ts>()... } : tuple_type{};
 	}
 
 	row_handle m_current;
@@ -422,7 +423,7 @@ class iterator_impl<Category, T>
   private:
 	value_type get() const
 	{
-		return m_current ?  m_current[m_item_ix].template as<value_type>() : value_type{};
+		return m_current ?  m_current[m_item_ix].template get<value_type>() : value_type{};
 	}
 
 	row_handle m_current;

include/cif++/matrix.hpp

@@ -124,6 +124,23 @@ class matrix_expression
 
 		return os;
 	}
+
+	template <typename M2>
+	constexpr bool operator==(const matrix_expression<M2> &m) const
+	{
+		bool same = false;
+		if (dim_m() == m.dim_m() and dim_n() == m.dim_n())
+		{
+			same = true;
+			for (std::size_t i = 0; same and i < m.dim_m(); ++i)
+			{
+				for (std::size_t j = 0; same and j < m.dim_n(); ++j)
+					same = operator()(i, j) == m(i, j);
+			}
+		}
+
+		return same;
+	}
 };
 
 // --------------------------------------------------------------------
@@ -594,6 +611,35 @@ auto operator*(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
 
 // --------------------------------------------------------------------
 
+template <typename M2>
+class sub_matrix : public matrix_expression<sub_matrix<M2>>
+{
+  public:
+	sub_matrix(const M2 &m, int i, int j)
+		: m_m(m)
+		, m_i(i)
+		, m_j(j)
+	{
+	}
+
+	constexpr std::size_t dim_m() const { return m_m.dim_m() - 1; } ///< Return dimension m
+	constexpr std::size_t dim_n() const { return m_m.dim_n() - 1; } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return m_m(
+			i >= m_i ? i + 1 : i,
+			j >= m_j ? j + 1 : j);
+	}
+
+  private:
+	const M2 &m_m;
+	std::size_t m_i, m_j;
+};
+
+// --------------------------------------------------------------------
+
 /** Generic routine to calculate the determinant of a matrix
  * 
  * @note This is currently only implemented for fixed matrices of size 3x3
@@ -605,11 +651,23 @@ auto determinant(const M &m);
 template <typename F = float>
 auto determinant(const matrix3x3<F> &m)
 {
-	return (m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) +
-			m(0, 1) * (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) +
-			m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)));
+	return (m(0, 0) * ((m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))) +
+			m(0, 1) * ((m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2))) +
+			m(0, 2) * ((m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0))));
 }
 
+/** Implementation of the determinant function for fixed size matrices of size 4x4 */
+template <typename F = float>
+F determinant(const matrix4x4<F> &m)
+{
+	return m(0, 0) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 0))) -
+	       m(0, 1) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 1))) +
+	       m(0, 2) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 2))) -
+	       m(0, 3) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 3)));
+}
+
+// --------------------------------------------------------------------
+
 /** Generic routine to calculate the inverse of a matrix
  * 
  * @note This is currently only implemented for fixed matrices of size 3x3

include/cif++/model.hpp

@@ -29,6 +29,7 @@
 #include "cif++/atom_type.hpp"
 #include "cif++/datablock.hpp"
 #include "cif++/point.hpp"
+#include "cif++/row.hpp"
 
 #include <memory>
 #include <numeric>
@@ -134,14 +135,20 @@ class atom
 
 		row_handle row_aniso()
 		{
+			row_handle result{};
 			auto cat = m_db.get("atom_site_anisotrop");
-			return cat ? cat->operator[]({ { "id", m_id } }) : row_handle{};
+			if (cat)
+				result = cat->operator[]({ { "id", m_id } });
+			return result;
 		}
 
 		const row_handle row_aniso() const
 		{
+			row_handle result{};
 			auto cat = m_db.get("atom_site_anisotrop");
-			return cat ? cat->operator[]({ { "id", m_id } }) : row_handle{};
+			if (cat)
+				result = cat->operator[]({ { "id", m_id } });
+			return result;
 		}
 
 		const datablock &m_db;
@@ -183,7 +190,7 @@ class atom
 	 * @param row The row containing the data for this atom
 	 */
 	atom(const datablock &db, const row_handle &row)
-		: atom(std::make_shared<atom_impl>(db, row["id"].as<std::string>()))
+		: atom(std::make_shared<atom_impl>(db, row["id"].get<std::string>()))
 	{
 	}
 
@@ -1059,12 +1066,30 @@ class structure
 	/// \return				The newly create asym ID
 	std::string create_non_poly(const std::string &entity_id, std::vector<row_initializer> atoms);
 
+	/// \brief Create a new NonPolymer struct_asym for a compound of type \a compound_id, returns asym_id.
+	/// This method creates new atom records filled with info from the CCD compound info.
+	///
+	/// \param compound_id	 The compound ID of the new nonpoly
+	/// \param skip_hydrogen Do not create hydrogen atoms when true
+	/// \return				 The newly create asym ID
+	std::string create_non_poly(const std::string &compound_id, bool skip_hydrogen);
+
 	/// \brief Create a new water with atom constructed from info in \a atom_info
 	/// This method creates a new atom record filled with info from the info.
 	///
 	/// \param atom			The set of item data containing the data for the atoms.
 	void create_water(row_initializer atom);
 
+	/// \brief Create a link, a struct_conn record for two atoms.
+	///
+	/// \param a1			Atom 1
+	/// \param a2			Atom 2
+	/// \param link_type	The struct_conn_type ID for the link
+	/// \param role			The pdbx_role field value
+	/// \return 			The ID of the struct_conn record created
+
+	std::string create_link(atom a1, atom a2, const std::string &link_type, const std::string &role);
+
 	/// \brief Create a new and empty (sugar) branch
 	branch &create_branch();
 

include/cif++/parser.hpp

@@ -164,7 +164,13 @@ class sac_parser
 		SAVE_NAME,
 		STOP,
 		ITEM_NAME,
-		VALUE
+
+		VALUE_INAPPLICABLE,
+		VALUE_UNKNOWN,
+		VALUE_NUMERIC_INTEGER,
+		VALUE_NUMERIC_FLOAT,
+		VALUE_CHARSTRING,
+		VALUE_TEXTFIELD
 	};
 
 	static constexpr const char *get_token_name(CIFToken token)
@@ -180,7 +186,15 @@ class sac_parser
 			case CIFToken::SAVE_NAME: return "SAVE+name";
 			case CIFToken::STOP: return "STOP";
 			case CIFToken::ITEM_NAME: return "Tag";
-			case CIFToken::VALUE: return "Value";
+			// case CIFToken::VALUE: return "Value";
+
+			case CIFToken::VALUE_INAPPLICABLE: return "Inapplicable value";
+			case CIFToken::VALUE_UNKNOWN: return "'Unknown' value (=null)";
+			case CIFToken::VALUE_NUMERIC_INTEGER: return "Integer value";
+			case CIFToken::VALUE_NUMERIC_FLOAT: return "Float value";
+			case CIFToken::VALUE_CHARSTRING: return "Charstring value";
+			case CIFToken::VALUE_TEXTFIELD: return "Textfield value";
+
 			default: return "Invalid token parameter";
 		}
 	}
@@ -262,7 +276,7 @@ class sac_parser
 	virtual void produce_datablock(std::string_view name) = 0;
 	virtual void produce_category(std::string_view name) = 0;
 	virtual void produce_row() = 0;
-	virtual void produce_item(std::string_view category, std::string_view item, std::string_view value) = 0;
+	virtual void produce_item(std::string_view category, std::string_view item, item_value value) = 0;
 
   protected:
 
@@ -281,7 +295,12 @@ class sac_parser
 		TextItem,
 		TextItemNL,
 		Reserved,
-		Value
+		Value,
+
+		Numeric_Integer,
+		Numeric_Float,
+		Numeric_Exponent1,
+		Numeric_Exponent2
 	};
 
 	std::streambuf &m_source;
@@ -294,6 +313,8 @@ class sac_parser
 	// token buffer
 	std::vector<char> m_token_buffer;
 	std::string_view m_token_value;
+	int64_t m_token_value_int;
+	double m_token_value_float;
 
 	/** @endcond */
 };
@@ -331,7 +352,7 @@ class parser : public sac_parser
 
 	void produce_row() override;
 
-	void produce_item(std::string_view category, std::string_view item, std::string_view value) override;
+	void produce_item(std::string_view category, std::string_view item, item_value value) override;
 
   protected:
 	file &m_file;

include/cif++/point.hpp

@@ -30,7 +30,9 @@
 #include <cmath>
 #include <complex>
 #include <cstdint>
+#include <format>
 #include <functional>
+#include <optional>
 #include <valarray>
 
 #if __has_include(<clipper/core/coords.h>)
@@ -365,11 +367,18 @@ class quaternion_type
 	}
 
 	/// \brief test for all zero values
-	constexpr operator bool() const
+	constexpr explicit operator bool() const
 	{
 		return a != 0 or b != 0 or c != 0 or d != 0;
 	}
 
+	/// \brief for debugging e.g.
+	friend std::ostream &operator<<(std::ostream &os, const quaternion_type &rhs)
+	{
+		os << std::format("{{ a: {}, b: {}, c: {}, d: {} }}", rhs.a, rhs.b, rhs.c, rhs.d);
+		return os;
+	}
+
   private:
 	value_type a, b, c, d;
 };
@@ -743,6 +752,55 @@ inline constexpr auto cross_product(const point_type<F1> &a, const point_type<F2
 		a.m_x * b.m_y - b.m_x * a.m_y);
 }
 
+/// \brief return the squared norm of point @a p
+template <typename F>
+constexpr F norm_squared(const point_type<F> &p)
+{
+	return p.m_x * p.m_x + p.m_y * p.m_y + p.m_z * p.m_z;
+}
+
+/// \brief return the norm of point @a p
+template <typename F>
+constexpr point_type<F> norm(const point_type<F> &p)
+{
+	return std::sqrt(norm_squared(p));
+}
+
+/// \brief return the point where two lines intersect, or an empty value if they don't intersect at all
+template <typename F>
+std::optional<cif::point> line_line_intersection(const point_type<F> &p1,
+	const point_type<F> &p2, const point_type<F> &p3, const point_type<F> &p4)
+{
+	auto p13 = p1 - p3;
+	auto p43 = p4 - p3;
+	if (std::abs(p43.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p43.m_z) < std::numeric_limits<F>::epsilon())
+		return {};
+
+	auto p21 = p2 - p1;
+	if (std::abs(p21.m_x) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_y) < std::numeric_limits<F>::epsilon() and std::abs(p21.m_z) < std::numeric_limits<F>::epsilon())
+		return {};
+
+	auto d1343 = cif::dot_product(p43, p13);
+	auto d4321 = cif::dot_product(p43, p21);
+	auto d1321 = cif::dot_product(p13, p21);
+	auto d4343 = cif::dot_product(p43, p43);
+	auto d2121 = cif::dot_product(p21, p21);
+
+	auto denom = d2121 * d4343 - d4321 * d4321;
+	if (std::abs(denom) < std::numeric_limits<F>::epsilon())
+		return {};
+
+	auto numer = d1343 * d4321 - d1321 * d4343;
+
+	auto mua = numer / denom;
+	auto mub = (d1343 + d4321 * mua) / d4343;
+
+	auto pa = p1 + mua * p21;
+	auto pb = p3 + mub * p43;
+
+	return { (pa + pb) / 2 };
+}
+
 /// \brief return the angle in degrees between the vectors from point @a p2 to @a p1 and @a p2 to @a p3
 template <typename F>
 constexpr auto angle(const point_type<F> &p1, const point_type<F> &p2, const point_type<F> &p3)
@@ -806,6 +864,9 @@ constexpr auto distance_point_to_line(const point_type<F> &l1, const point_type<
 	return cross.length() / line.length();
 }
 
+/// \brief return the smallest sphere around the points in @a pts
+std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts);
+
 // --------------------------------------------------------------------
 /**
  * @brief For e.g. simulated annealing, returns a new point that is moved in

include/cif++/row.hpp

@@ -41,7 +41,7 @@
  * that's not the case.
  * 
  * You can access the values of stored items by name or index.
- * The return value of operator[] is an cif::item_handle object.
+ * The return value of operator[] is a reference to a cif::item_value object.
  * 
  * @code {.cpp}
  * cif::category &atom_site = my_db["atom_site"];
@@ -93,7 +93,7 @@ namespace detail
 		{
 		}
 
-		const item_handle operator[](uint16_t ix) const
+		const item_value &operator[](uint16_t ix) const
 		{
 			return m_row[m_items[ix]];
 		}
@@ -107,7 +107,7 @@ namespace detail
 		template <typename... Ts, std::size_t... Is>
 		std::tuple<Ts...> get(std::index_sequence<Is...>) const
 		{
-			return std::tuple<Ts...>{ m_row[m_items[Is]].template as<Ts>()... };
+			return std::tuple<Ts...>{ m_row[m_items[Is]].template get<Ts>()... };
 		}
 
 		const row_handle &m_row;
@@ -173,14 +173,14 @@ class row : public std::vector<item_value>
 		return ix < size() ? &data()[ix] : nullptr;
 	}
 
-  private:
+//   private:
 	friend class category;
 	friend class category_index;
 
 	template <typename, typename...>
 	friend class iterator_impl;
 
-	void append(uint16_t ix, item_value &&iv)
+	void append(uint16_t ix, item_value iv)
 	{
 		if (ix >= size())
 			resize(ix + 1);
@@ -204,7 +204,6 @@ class row_handle
 {
   public:
 	/** @cond */
-	friend struct item_handle;
 	friend class category;
 	friend class category_index;
 	friend class row_initializer;
@@ -245,29 +244,20 @@ class row_handle
 		return not empty();
 	}
 
-	/// \brief return a cif::item_handle to the item in item @a item_ix
-	item_handle operator[](uint16_t item_ix)
-	{
-		return empty() ? item_handle::s_null_item : item_handle(item_ix, *this);
-	}
+	/// \brief return the count of the items
+	[[nodiscard]] size_t size() const { return m_row->size(); }
 
-	/// \brief return a const cif::item_handle to the item in item @a item_ix
-	const item_handle operator[](uint16_t item_ix) const
-	{
-		return empty() ? item_handle::s_null_item : item_handle(item_ix, const_cast<row_handle &>(*this));
-	}
+	/// \brief return a reference to a cif::item_value to the item in item @a item_ix
+	item_value &operator[](uint16_t item_ix);
 
-	/// \brief return a cif::item_handle to the item in the item named @a item_name
-	item_handle operator[](std::string_view item_name)
-	{
-		return empty() ? item_handle::s_null_item : item_handle(add_item(item_name), *this);
-	}
+	/// \brief return a const reference to a cif::item_value to the item in item @a item_ix
+	const item_value &operator[](uint16_t item_ix) const;
 
-	/// \brief return a const cif::item_handle to the item in the item named @a item_name
-	const item_handle operator[](std::string_view item_name) const
-	{
-		return empty() ? item_handle::s_null_item : item_handle(get_item_ix(item_name), const_cast<row_handle &>(*this));
-	}
+	/// \brief return a reference to a cif::item_value to the item in the item named @a item_name
+	item_value &operator[](std::string_view item_name);
+
+	/// \brief return a const reference to a cif::item_value to the item in the item named @a item_name
+	const item_value &operator[](std::string_view item_name) const;
 
 	/// \brief Return an object that can be used in combination with cif::tie
 	/// to assign the values for the items @a items
@@ -288,14 +278,14 @@ class row_handle
 	template <typename T>
 	T get(const char *item) const
 	{
-		return operator[](get_item_ix(item)).template as<T>();
+		return operator[](get_item_ix(item)).template get<T>();
 	}
 
 	/// \brief Get the value of item @a item cast to type @a T
 	template <typename T>
 	T get(std::string_view item) const
 	{
-		return operator[](get_item_ix(item)).template as<T>();
+		return operator[](get_item_ix(item)).template get<T>();
 	}
 
 	/// \brief assign each of the items named in @a values to their respective value
@@ -316,9 +306,9 @@ class row_handle
 	 * checked to see if it conforms to the rules defined in the dictionary
 	 */
 
-	void assign(std::string_view name, std::string_view value, bool updateLinked, bool validate = true)
+	void assign(std::string_view name, item_value value, bool updateLinked, bool validate = true)
 	{
-		assign(add_item(name), value, updateLinked, validate);
+		assign(add_item(name), std::move(value), updateLinked, validate);
 	}
 
 	/** \brief assign the value @a value to item at index @a item
@@ -332,7 +322,7 @@ class row_handle
 	 * checked to see if it conforms to the rules defined in the dictionary
 	 */
 
-	void assign(uint16_t item, std::string_view value, bool updateLinked, bool validate = true);
+	void assign(uint16_t item, item_value value, bool updateLinked, bool validate = true);
 
 	/// \brief compare two rows
 	bool operator==(const row_handle &rhs) const { return m_category == rhs.m_category and m_row == rhs.m_row; }
@@ -356,10 +346,10 @@ class row_handle
 		return m_row;
 	}
 
-	void assign(const item &i, bool updateLinked)
+	void assign(const item &i, bool updateLinked);/* 
 	{
 		assign(i.name(), i.value(), updateLinked);
-	}
+	} */
 
 	void swap(uint16_t item, row_handle &r);
 
@@ -408,7 +398,7 @@ class row_initializer : public std::vector<item>
 
 
 	/// \brief set the value for item name @a name to @a value
-	void set_value(std::string_view name, std::string_view value);
+	void set_value(std::string name, item_value value);
 
 	/// \brief set the value for item based on @a i
 	void set_value(const item &i)
@@ -417,7 +407,7 @@ class row_initializer : public std::vector<item>
 	}
 
 	/// \brief set the value for item name @a name to @a value, but only if the item did not have a value already
-	void set_value_if_empty(std::string_view name, std::string_view value);
+	void set_value_if_empty(std::string name, item_value value);
 
 	/// \brief set the value for item @a i, but only if the item did not have a value already
 	void set_value_if_empty(const item &i)

include/cif++/text.hpp

@@ -355,279 +355,35 @@ std::string cif_id_for_number(int number);
 std::vector<std::string> word_wrap(const std::string &text, std::size_t width);
 
 // --------------------------------------------------------------------
-/// \brief std::from_chars for floating point types.
-///
-/// These are optional, there's a selected_charconv class below that selects
-/// the best option to use based on support by the stl library.
-///
-/// I.e. that in case of GNU < 12 (or something) the cif implementation will
-/// be used, all other cases will use the stl version.
 
-template <typename FloatType, std::enable_if_t<std::is_floating_point_v<FloatType>, int> = 0>
-std::from_chars_result from_chars(const char *first, const char *last, FloatType &value)
-{
-	std::from_chars_result result{ first, {} };
-
-	enum State
-	{
-		IntegerSign,
-		Integer,
-		Fraction,
-		ExponentSign,
-		Exponent
-	} state = IntegerSign;
-	int sign = 1;
-	unsigned long long vi = 0;
-	int fl = 0, tz = 0;
-	int exponent_sign = 1;
-	int exponent = 0;
-	bool done = false;
-
-	while (not done and not (bool)result.ec)
-	{
-		char ch = result.ptr != last ? *result.ptr : 0;
-		++result.ptr;
-
-		switch (state)
-		{
-			case IntegerSign:
-				if (ch == '-')
-				{
-					sign = -1;
-					state = Integer;
-				}
-				else if (ch == '+')
-					state = Integer;
-				else if (ch >= '0' and ch <= '9')
-				{
-					vi = ch - '0';
-					state = Integer;
-				}
-				else if (ch == '.')
-					state = Fraction;
-				else
-					result.ec = std::errc::invalid_argument;
-				break;
-
-			case Integer:
-				if (ch >= '0' and ch <= '9')
-					vi = 10 * vi + (ch - '0');
-				else if (ch == 'e' or ch == 'E')
-					state = ExponentSign;
-				else if (ch == '.')
-					state = Fraction;
-				else
-				{
-					done = true;
-					--result.ptr;
-				}
-				break;
-
-			case Fraction:
-				if (ch >= '0' and ch <= '9')
-				{
-					vi = 10 * vi + (ch - '0');
-
-					if (ch == '0')
-						tz += 1;
-					else
-					{
-						fl += tz + 1;
-						tz = 0;
-					}
-				}
-				else if (ch == 'e' or ch == 'E')
-					state = ExponentSign;
-				else
-				{
-					done = true;
-					--result.ptr;
-				}
-				break;
-
-			case ExponentSign:
-				if (ch == '-')
-				{
-					exponent_sign = -1;
-					state = Exponent;
-				}
-				else if (ch == '+')
-					state = Exponent;
-				else if (ch >= '0' and ch <= '9')
-				{
-					exponent = ch - '0';
-					state = Exponent;
-				}
-				else
-					result.ec = std::errc::invalid_argument;
-				break;
-
-			case Exponent:
-				if (ch >= '0' and ch <= '9')
-					exponent = 10 * exponent + (ch - '0');
-				else
-				{
-					done = true;
-					--result.ptr;
-				}
-				break;
-		}
-	}
-
-	if (not (bool)result.ec)
-	{
-		while (tz-- > 0)
-			vi /= 10;
-
-		long double v = std::pow(10, -fl) * vi * sign;
-		if (exponent != 0)
-			v *= std::pow(10, exponent * exponent_sign);
-
-		if (std::isnan(v))
-			result.ec = std::errc::invalid_argument;
-		else if (std::abs(v) > std::numeric_limits<FloatType>::max())
-			result.ec = std::errc::result_out_of_range;
-
-		value = static_cast<FloatType>(v);
-	}
-
-	return result;
-}
-
-/// \brief duplication of std::chars_format for deficient STL implementations
-enum class chars_format
-{
-	scientific = 1,
-	fixed = 2,
-	// hex,
-	general = fixed | scientific
-};
-
-/// \brief a simplistic implementation of std::to_chars for old STL implementations
-template <typename FloatType, std::enable_if_t<std::is_floating_point_v<FloatType>, int> = 0>
-std::to_chars_result to_chars(char *first, char *last, FloatType &value, chars_format fmt)
-{
-	int size = static_cast<int>(last - first);
-	int r = 0;
-
-	switch (fmt)
-	{
-		case chars_format::scientific:
-			if constexpr (std::is_same_v<FloatType, long double>)
-				r = snprintf(first, last - first, "%le", value);
-			else
-				r = snprintf(first, last - first, "%e", value);
-			break;
-
-		case chars_format::fixed:
-			if constexpr (std::is_same_v<FloatType, long double>)
-				r = snprintf(first, last - first, "%lf", value);
-			else
-				r = snprintf(first, last - first, "%f", value);
-			break;
-
-		case chars_format::general:
-			if constexpr (std::is_same_v<FloatType, long double>)
-				r = snprintf(first, last - first, "%lg", value);
-			else
-				r = snprintf(first, last - first, "%g", value);
-			break;
-	}
-
-	std::to_chars_result result;
-	if (r < 0 or r >= size)
-		result = { first, std::errc::value_too_large };
-	else
-		result = { first + r, std::errc() };
-
-	return result;
-}
-
-/// \brief a simplistic implementation of std::to_chars for old STL implementations
-template <typename FloatType, std::enable_if_t<std::is_floating_point_v<FloatType>, int> = 0>
-std::to_chars_result to_chars(char *first, char *last, FloatType &value, chars_format fmt, int precision)
-{
-	int size = static_cast<int>(last - first);
-	int r = 0;
-
-	switch (fmt)
-	{
-		case chars_format::scientific:
-			if constexpr (std::is_same_v<FloatType, long double>)
-				r = snprintf(first, last - first, "%.*le", precision, value);
-			else
-				r = snprintf(first, last - first, "%.*e", precision, value);
-			break;
-
-		case chars_format::fixed:
-			if constexpr (std::is_same_v<FloatType, long double>)
-				r = snprintf(first, last - first, "%.*lf", precision, value);
-			else
-				r = snprintf(first, last - first, "%.*f", precision, value);
-			break;
-
-		case chars_format::general:
-			if constexpr (std::is_same_v<FloatType, long double>)
-				r = snprintf(first, last - first, "%.*lg", precision, value);
-			else
-				r = snprintf(first, last - first, "%.*g", precision, value);
-			break;
-	}
-
-	std::to_chars_result result;
-	if (r < 0 or r >= size)
-		result = { first, std::errc::value_too_large };
-	else
-		result = { first + r, std::errc() };
-
-	return result;
-}
-
-/// \brief class that uses our implementation of std::from_chars and std::to_chars
 template <typename T>
-struct my_charconv
-{
-	/// @brief Simply call our version of std::from_chars
-	static std::from_chars_result from_chars(const char *a, const char *b, T &d)
-	{
-		return cif::from_chars(a, b, d);
-	}
+using from_chars_function = decltype(std::from_chars(std::declval<const char *>(), std::declval<const char *>(), std::declval<T &>()));
 
-	/// @brief Simply call our version of std::to_chars
-	static std::to_chars_result to_chars(char *first, char *last, T &value, chars_format fmt)
-	{
-		return cif::to_chars(first, last, value, fmt);
-	}
-};
-
-/// \brief class that uses the STL implementation of std::from_chars and std::to_chars
 template <typename T>
 struct std_charconv
 {
-	/// @brief Simply call std::from_chars
 	static std::from_chars_result from_chars(const char *a, const char *b, T &d)
 	{
 		return std::from_chars(a, b, d);
 	}
-
-	/// @brief Simply call std::to_chars
-	static std::to_chars_result to_chars(char *first, char *last, T &value, chars_format fmt)
-	{
-		return std::to_chars(first, last, value, fmt);
-	}
 };
 
-/// \brief helper to find a from_chars function
-template <typename T>
-using from_chars_function = decltype(std::from_chars(std::declval<const char *>(), std::declval<const char *>(), std::declval<T &>()));
+template <typename T, typename = void>
+struct ff_charconv;
 
-/**
- * @brief Helper to select the best implementation of charconv based on availability of the
- * function in the std:: namespace
- *
- * @tparam T The type for which we want to find a from_chars/to_chars function
- */
 template <typename T>
-using selected_charconv = typename std::conditional_t<std_experimental::is_detected_v<from_chars_function, T>, std_charconv<T>, my_charconv<T>>;
+struct ff_charconv<T, typename std::enable_if_t<std::is_floating_point_v<T>>>
+{
+	static std::from_chars_result from_chars(const char *a, const char *b, T &v);
+};
+
+template <typename T>
+using charconv = typename std::conditional_t<std_experimental::is_detected_v<from_chars_function, T>, std_charconv<T>, ff_charconv<T>>;
+
+template <typename T>
+constexpr auto from_chars(const char *s, const char *e, T &v)
+{
+	return charconv<T>::from_chars(s, e, v);
+}
 
 } // namespace cif

include/cif++/utilities.hpp

@@ -296,6 +296,11 @@ class progress_bar
 	 */
 	void message(const std::string &inMessage);
 
+	/**
+	 * @brief Flush the progress bar to the output stream
+	 */
+	void flush();
+
   private:
 	progress_bar(const progress_bar &) = delete;
 	progress_bar &operator=(const progress_bar &) = delete;

include/cif++/validate.hpp

@@ -329,10 +329,10 @@ struct item_validator
 
 	/// @brief Validate the value in @a value for this item
 	/// Will throw a std::system_error exception if it fails
-	void operator()(std::string_view value) const;
+	void operator()(const item_value &value) const;
 
 	/// @brief A more gentle version of value validation
-	bool validate_value(std::string_view value, std::error_code &ec) const noexcept;
+	bool validate_value(const item_value &value, std::error_code &ec) const noexcept;
 };
 
 /**

src/compound.cpp

@@ -24,14 +24,38 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include "cif++.hpp"
+#include "cif++/compound.hpp" // for compound_atom, compound_bond, compoun...
 
-#include <filesystem>
-#include <fstream>
-#include <map>
-#include <mutex>
-#include <numeric>
-#include <shared_mutex>
+#include "cif++/atom_type.hpp" // for atom_type_traits
+#include "cif++/category.hpp"  // for category
+#include "cif++/datablock.hpp" // for datablock
+#include "cif++/file.hpp"      // for file
+#include "cif++/item.hpp"      // for item
+#include "cif++/iterator.hpp"  // for iterator_proxy
+#include "cif++/parser.hpp"    // for parser
+#include "cif++/point.hpp"     // for distance, point
+#include "cif++/row.hpp"       // for tie, row_initializer, tie_wrap
+#include "cif++/text.hpp"      // for iequals, replace_all, iset
+#include "cif++/utilities.hpp" // for load_resource, VERBOSE, colour_type
+
+#include <algorithm>    // for find_if
+#include <cstddef>      // for size_t
+#include <exception>    // for exception, throw_with_nested
+#include <filesystem>   // for path, exists
+#include <fstream>      // for char_traits, basic_ostream, operator<<
+#include <iomanip>      // for operator<<, quoted
+#include <iostream>     // for clog, cout, cerr
+#include <limits>       // for numeric_limits
+#include <list>         // for _List_iterator
+#include <map>          // for allocator, map, _Rb_tree_iterator
+#include <memory>       // for shared_ptr, unique_ptr, __shared_ptr_...
+#include <optional>     // for optional
+#include <shared_mutex> // for shared_lock, shared_timed_mutex
+#include <stdexcept>    // for runtime_error, invalid_argument, out_...
+#include <string>       // for basic_string, string, operator==, ope...
+#include <string_view>  // for string_view, basic_string_view
+#include <utility>      // for pair, exchange, move
+#include <vector>       // for vector
 
 namespace fs = std::filesystem;
 
@@ -140,7 +164,7 @@ compound::compound(cif::datablock &db)
 
 	cif::tie(m_id, m_name, m_type, m_formula, m_formula_weight, m_formal_charge, one_letter_code, m_parent_id) =
 		chemComp.front().get("id", "name", "type", "formula", "formula_weight", "pdbx_formal_charge", "one_letter_code", "mon_nstd_parent_comp_id");
-	
+
 	if (one_letter_code.length() == 1)
 		m_one_letter_code = one_letter_code.front();
 
@@ -159,7 +183,7 @@ compound::compound(cif::datablock &db)
 		if (stereo_config.empty())
 			atom.stereo_config = stereo_config_type::N;
 		else
-		atom.stereo_config = parse_stereo_config_from_string(stereo_config);
+			atom.stereo_config = parse_stereo_config_from_string(stereo_config);
 		m_atoms.push_back(std::move(atom));
 	}
 
@@ -172,7 +196,7 @@ compound::compound(cif::datablock &db)
 		if (valueOrder.empty())
 			bond.type = bond_type::sing;
 		else
-		bond.type = parse_bond_type_from_string(valueOrder);
+			bond.type = parse_bond_type_from_string(valueOrder);
 		m_bonds.push_back(std::move(bond));
 	}
 }
@@ -231,12 +255,12 @@ float compound::bond_length(const std::string &atomId_1, const std::string &atom
 
 bool compound::is_peptide() const
 {
-	return iequals(m_type, "l-peptide linking")	or iequals(m_type, "peptide linking");
+	return iequals(m_type, "l-peptide linking") or iequals(m_type, "peptide linking");
 }
 
 bool compound::is_base() const
 {
-	return iequals(m_type, "dna linking")	or iequals(m_type, "rna linking");
+	return iequals(m_type, "dna linking") or iequals(m_type, "rna linking");
 }
 
 // --------------------------------------------------------------------
@@ -294,12 +318,31 @@ class compound_factory_impl : public std::enable_shared_from_this<compound_facto
 			delete c;
 	}
 
+	virtual bool exists_self(const std::string &id) const
+	{
+		if (m_missing.contains(id))
+			return false;
+
+		if (std::find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c)
+				{ return c->id() == id; }) != m_compounds.end())
+			return true;
+
+		return m_next and m_next->exists_self(id);
+	}
+
+	bool exists(std::string_view id)
+	{
+		std::shared_lock lock(mMutex);
+
+		return exists_self(std::string{ id });
+	}
+
 	compound *get(std::string id)
 	{
 		std::shared_lock lock(mMutex);
 
 		compound *result = nullptr;
-		
+
 		for (auto impl = shared_from_this(); impl; impl = impl->m_next)
 		{
 			result = impl->create(id);
@@ -363,7 +406,9 @@ compound *compound_factory_impl::create(const std::string &id)
 	if (m_missing.contains(id))
 		return nullptr;
 
-	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) { return c->id() == id; }); i != m_compounds.end())
+	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c)
+			{ return c->id() == id; });
+		i != m_compounds.end())
 		return *i;
 
 	compound *result = nullptr;
@@ -454,7 +499,6 @@ class local_compound_factory_impl : public compound_factory_impl
 	compound *create(const std::string &id) override;
 
   private:
-
 	compound *construct_compound(const datablock &db, const std::string &id, const std::string &name, const std::string &three_letter_code, const std::string &group);
 
 	cif::file m_local_file;
@@ -465,7 +509,9 @@ compound *local_compound_factory_impl::create(const std::string &id)
 	if (m_missing.contains(id))
 		return nullptr;
 
-	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c) { return c->id() == id; }); i != m_compounds.end())
+	if (auto i = find_if(m_compounds.begin(), m_compounds.end(), [id](compound *c)
+			{ return c->id() == id; });
+		i != m_compounds.end())
 		return *i;
 
 	compound *result = nullptr;
@@ -510,12 +556,10 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 
 	float formula_weight = 0;
 	int formal_charge = 0;
-	std::map<std::string,std::size_t> formula_data;
+	std::map<std::string, std::size_t> formula_data;
 
 	for (std::size_t ord = 1; const auto &[atom_id, type_symbol, type, charge, x, y, z, xi, yi, zi] :
-		rdb["chem_comp_atom"].rows<std::string, std::string, std::string, int,
-			std::optional<float>, std::optional<float>, std::optional<float>,
-			std::optional<float>, std::optional<float>, std::optional<float>>(
+		rdb["chem_comp_atom"].rows<std::string, std::string, std::string, int, std::optional<float>, std::optional<float>, std::optional<float>, std::optional<float>, std::optional<float>, std::optional<float>>(
 			"atom_id", "type_symbol", "type", "charge",
 			"model_Cartn_x", "model_Cartn_y", "model_Cartn_z",
 			"pdbx_model_Cartn_x_ideal", "pdbx_model_Cartn_y_ideal", "pdbx_model_Cartn_z_ideal"))
@@ -525,16 +569,14 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 
 		formula_data[type_symbol] += 1;
 
-		db["chem_comp_atom"].emplace({
-			{ "comp_id", id },
+		db["chem_comp_atom"].emplace({ { "comp_id", id },
 			{ "atom_id", atom_id },
 			{ "type_symbol", type_symbol },
 			{ "charge", charge },
-			{ "model_Cartn_x",  x.has_value() ? x : xi, 3 },
-			{ "model_Cartn_y",  y.has_value() ? y : yi, 3 },
-			{ "model_Cartn_z",  z.has_value() ? z : zi, 3 },
-			{ "pdbx_ordinal", ord++ }
-		});
+			{ "model_Cartn_x", x.has_value() ? x : xi, 3 },
+			{ "model_Cartn_y", y.has_value() ? y : yi, 3 },
+			{ "model_Cartn_z", z.has_value() ? z : zi, 3 },
+			{ "pdbx_ordinal", ord++ } });
 
 		formal_charge += charge;
 	}
@@ -551,21 +593,19 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 		else if (cif::iequals(type, "triple") or cif::iequals(type, "trip"))
 			value_order = "TRIP";
 
-		db["chem_comp_bond"].emplace({
-			{ "comp_id", id },
+		db["chem_comp_bond"].emplace({ { "comp_id", id },
 			{ "atom_id_1", atom_id_1 },
 			{ "atom_id_2", atom_id_2 },
 			{ "value_order", value_order },
 			{ "pdbx_aromatic_flag", aromatic },
 			// TODO: fetch stereo_config info from chem_comp_chir
-			{ "pdbx_ordinal", ord++ }
-		});
+			{ "pdbx_ordinal", ord++ } });
 	}
 
 	db.emplace_back(rdb["pdbx_chem_comp_descriptor"]);
 
 	std::string formula;
-	for (bool first = true; const auto &[symbol, count]: formula_data)
+	for (bool first = true; const auto &[symbol, count] : formula_data)
 	{
 		if (std::exchange(first, false))
 			formula += ' ';
@@ -584,15 +624,13 @@ compound *local_compound_factory_impl::construct_compound(const datablock &rdb,
 	else
 		type = "NON-POLYMER";
 
-	db["chem_comp"].emplace({
-		{ "id", id },
+	db["chem_comp"].emplace({ { "id", id },
 		{ "name", name },
 		{ "type", type },
 		{ "formula", formula },
 		{ "pdbx_formal_charge", formal_charge },
 		{ "formula_weight", formula_weight },
-		{ "three_letter_code", three_letter_code }
-	});
+		{ "three_letter_code", three_letter_code } });
 
 	std::shared_lock lock(mMutex);
 
@@ -698,6 +736,11 @@ void compound_factory::pop_dictionary()
 		m_impl = m_impl->next();
 }
 
+bool compound_factory::exists(std::string_view id) const
+{
+	return m_impl and m_impl->exists(id);
+}
+
 const compound *compound_factory::create(std::string_view id)
 {
 	auto result = m_impl ? m_impl->get(std::string{ id }) : nullptr;
@@ -722,7 +765,7 @@ bool compound_factory::is_peptide(std::string_view res_name) const
 	bool result = is_std_peptide(res_name);
 	if (not result and m_impl)
 	{
-		auto compound = const_cast<compound_factory&>(*this).create(res_name);
+		auto compound = const_cast<compound_factory &>(*this).create(res_name);
 		result = compound != nullptr and compound->is_peptide();
 	}
 	return result;
@@ -734,7 +777,7 @@ bool compound_factory::is_base(std::string_view res_name) const
 	bool result = is_std_base(res_name);
 	if (not result and m_impl)
 	{
-		auto compound = const_cast<compound_factory&>(*this).create(res_name);
+		auto compound = const_cast<compound_factory &>(*this).create(res_name);
 		result = compound != nullptr and compound->is_base();
 	}
 	return result;
@@ -784,7 +827,7 @@ void compound_factory::report_missing_compound(std::string_view compound_id)
 				  << "in /var/cache/libcifpp using the following commands:\n\n"
 				  << "curl -o " << CACHE_DIR << "/components.cif https://files.wwpdb.org/pub/pdb/data/monomers/components.cif\n"
 				  << "curl -o " << CACHE_DIR << "/mmcif_pdbx.dic https://mmcif.wwpdb.org/dictionaries/ascii/mmcif_pdbx_v50.dic\n"
-				  << "curl -o " << CACHE_DIR << "/mmcif_ma.dic https://github.com/ihmwg/ModelCIF/raw/master/dist/mmcif_ma.dic\n\n";
+				  << "curl -o " << CACHE_DIR << "/mmcif_ma.dic https://mmcif.wwpdb.org/dictionaries/ascii/mmcif_ma.dic\n\n";
 #endif
 
 		if (m_impl)

src/condition.cpp

@@ -123,21 +123,6 @@ namespace detail
 		return this;
 	}
 
-	condition_impl *key_equals_number_condition_impl::prepare(const category &c)
-	{
-		m_item_ix = c.get_item_ix(m_item_name);
-
-		if (c.get_cat_validator() != nullptr and
-			c.key_item_indices().contains(m_item_ix) and
-			c.key_item_indices().size() == 1)
-		{
-			item v(m_item_name, m_value);
-			m_single_hit = c[{ { m_item_name, std::string{ v.value() }, false } }];
-		}
-
-		return this;
-	}
-
 	bool found_in_range(condition_impl *c, std::vector<and_condition_impl *>::iterator b, std::vector<and_condition_impl *>::iterator e)
 	{
 		bool result = true;
@@ -234,17 +219,6 @@ namespace detail
 					continue;
 				}
 
-				if (auto s = dynamic_cast<const key_equals_number_condition_impl *>(sub); s != nullptr)
-				{
-					if (keys.contains(s->m_item_name))
-					{
-						item v{ s->m_item_name, s->m_value };
-						lookup.emplace_back(s->m_item_name, std::string{ v.value() } );
-						subs.emplace_back(sub);
-					}
-					continue;
-				}
-
 				if (auto s = dynamic_cast<const key_equals_or_empty_condition_impl *>(sub); s != nullptr)
 				{
 					if (keys.contains(s->m_item_name))
@@ -255,17 +229,6 @@ namespace detail
 					}
 					continue;
 				}
-
-				if (auto s = dynamic_cast<const key_equals_number_or_empty_condition_impl *>(sub); s != nullptr)
-				{
-					if (keys.contains(s->m_item_name))
-					{
-						item v{ s->m_item_name, s->m_value };
-						lookup.emplace_back(s->m_item_name, std::string{ v.value() }, true );
-						subs.emplace_back(sub);
-					}
-					continue;
-				}
 			}
 
 			if (lookup.size() == keys.size())

src/datablock.cpp

@@ -25,8 +25,11 @@
  */
 
 #include "cif++/datablock.hpp"
+
 #include "cif++/validate.hpp"
 
+#include <exception>
+
 namespace cif
 {
 
@@ -42,7 +45,16 @@ datablock::datablock(const datablock &db)
 void datablock::load_dictionary()
 {
 	if (auto *audit_conform = get("audit_conform"); audit_conform and not audit_conform->empty())
-		set_validator(&validator_factory::instance().get(*audit_conform));
+	{
+		try
+		{
+			set_validator(&validator_factory::instance().get(*audit_conform));
+		}
+		catch (const std::exception &ex)
+		{
+			std::clog << ex.what() << '\n';
+		}
+	}
 }
 
 void datablock::set_validator(const validator *v)
@@ -96,7 +108,8 @@ bool datablock::strip()
 	bool result = true;
 
 	// remove all categories that have no validator
-	erase(std::remove_if(begin(), end(), [](category &c) {
+	erase(std::remove_if(begin(), end(), [](category &c)
+			  {
 		bool result = false;
 		if (c.get_cat_validator() == nullptr)
 		{
@@ -104,8 +117,8 @@ bool datablock::strip()
 				std::clog << "Dropping category " << c.name() << '\n';
 			result = true;
 		}
-		return result;
-	}), end());
+		return result; }),
+		end());
 
 	// then strip the remaining categories
 	for (auto &cat : *this)

src/dictionary_parser.cpp

@@ -28,6 +28,9 @@
 #include "cif++/dictionary_parser.hpp"
 #include "cif++/file.hpp"
 #include "cif++/parser.hpp"
+#include <exception>
+#include <iomanip>
+#include <stdexcept>
 
 namespace cif
 {
@@ -46,7 +49,7 @@ class dictionary_parser : public parser
 	void load_dictionary()
 	{
 		std::unique_ptr<datablock> dict;
-		auto savedDatablock = m_datablock;
+		auto savedDatablock = std::exchange(m_datablock, nullptr);
 
 		try
 		{
@@ -75,6 +78,9 @@ class dictionary_parser : public parser
 			error(ex.what());
 		}
 
+		if (m_datablock == nullptr)
+			throw std::runtime_error("Dictionary file is empty?");
+
 		// store all validators
 		for (auto &ic : mCategoryValidators)
 			m_validator.add_category_validator(std::move(ic));
@@ -147,7 +153,7 @@ class dictionary_parser : public parser
 					match(CIFToken::ITEM_NAME);
 				}
 
-				while (m_lookahead == CIFToken::VALUE)
+				while (m_lookahead >= CIFToken::VALUE_INAPPLICABLE)
 				{
 					cat->emplace({});
 					auto row = cat->back();
@@ -155,7 +161,7 @@ class dictionary_parser : public parser
 					for (auto item_name : item_names)
 					{
 						row[item_name] = m_token_value;
-						match(CIFToken::VALUE);
+						match(m_lookahead);
 					}
 				}
 
@@ -175,7 +181,7 @@ class dictionary_parser : public parser
 					cat->emplace({});
 				cat->back()[item_name] = m_token_value;
 
-				match(CIFToken::VALUE);
+				match(m_lookahead >= CIFToken::VALUE_INAPPLICABLE ? m_lookahead : CIFToken::VALUE_CHARSTRING);
 			}
 		}
 
@@ -187,11 +193,11 @@ class dictionary_parser : public parser
 
 			std::vector<std::string> keys;
 			for (auto k : dict["category_key"])
-				keys.push_back(std::get<1>(split_item_name(k["name"].as<std::string>())));
+				keys.push_back(std::get<1>(split_item_name(k["name"].get<std::string>())));
 
 			iset groups;
 			for (auto g : dict["category_group"])
-				groups.insert(g["id"].as<std::string>());
+				groups.insert(g["id"].get<std::string>());
 
 			mCategoryValidators.push_back(category_validator{ category, keys, groups });
 		}
@@ -206,7 +212,7 @@ class dictionary_parser : public parser
 
 			iset ess;
 			for (auto e : dict["item_enumeration"])
-				ess.insert(e["value"].as<std::string>());
+				ess.insert(e["value"].get<std::string>());
 
 			std::string defaultValue = dict["item_default"].front().get<std::string>("value");
 			// bool defaultIsNull = false;
@@ -399,7 +405,7 @@ class dictionary_parser : public parser
 			// look up the label
 			for (auto r : linkedGroup.find("category_id"_key == link.m_child_category and "link_group_id"_key == link.m_link_group_id))
 			{
-				link.m_link_group_label = r["label"].as<std::string>();
+				link.m_link_group_label = r["label"].get<std::string>();
 				break;
 			}
 

src/item.cpp

@@ -24,45 +24,109 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/item.hpp"
+
 #include "cif++/row.hpp"
 
 #include <cassert>
+#include <compare>
+#include <ios>
 
 namespace cif
 {
 
-const item_handle item_handle::s_null_item;
-row_handle s_null_row_handle;
-
-item_handle::item_handle()
-	: m_item_ix(std::numeric_limits<uint16_t>::max())
-	, m_row_handle(s_null_row_handle)
+int item_value::compare(const item_value &b, bool ignore_case) const noexcept
 {
-}
+	int d = static_cast<int>(m_data.m_type) - static_cast<int>(b.m_data.m_type);
 
-std::string_view item_handle::text() const
-{
-	if (not m_row_handle.empty())
+	if (d == 0)
 	{
-		auto iv = m_row_handle.m_row->get(m_item_ix);
-		if (iv != nullptr)
-			return iv->text();
+		switch (m_data.m_type)
+		{
+			case cif::item_value_type::BOOLEAN:
+				d = static_cast<int>(m_data.m_value.m_boolean) - static_cast<int>(b.m_data.m_value.m_boolean);
+				break;
+			case cif::item_value_type::INT:
+				d = m_data.m_value.m_integer - b.m_data.m_value.m_integer;
+				break;
+			case cif::item_value_type::FLOAT:
+			{
+				auto dp = (m_data.m_value.m_float <=> b.m_data.m_value.m_float);
+				if (dp == std::partial_ordering::less)
+					d = -1;
+				else if (dp == std::partial_ordering::greater)
+					d = 1;
+				break;
+			}
+			case cif::item_value_type::TEXT:
+				d = m_data.sv().compare(b.m_data.sv());
+				break;
+			default:;
+		}
 	}
 
-	return {};
+	return d;
 }
 
-void item_handle::assign_value(std::string_view value)
+// const item_handle item_handle::s_null_item;
+// row_handle s_null_row_handle;
+
+// item_handle::item_handle()
+// 	: m_item_ix(std::numeric_limits<uint16_t>::max())
+// 	, m_row_handle(s_null_row_handle)
+// {
+// }
+
+// std::string_view item_handle::text() const
+// {
+// 	if (not m_row_handle.empty())
+// 	{
+// 		auto iv = m_row_handle.m_row->get(m_item_ix);
+// 		if (iv != nullptr)
+// 			return iv->text();
+// 	}
+
+// 	return {};
+// }
+
+// void item_handle::assign_value(std::string_view value)
+// {
+// 	assert(not m_row_handle.empty());
+// 	m_row_handle.assign(m_item_ix, value, true);
+// }
+
+// void item_handle::swap(item_handle &b)
+// {
+// 	assert(m_item_ix == b.m_item_ix);
+// 	// assert(&m_row_handle.m_category == &b.m_row_handle.m_category);
+// 	m_row_handle.swap(m_item_ix, b.m_row_handle);
+// }
+
+std::ostream &operator<<(std::ostream &os, const item_value &v)
 {
-	assert(not m_row_handle.empty());
-	m_row_handle.assign(m_item_ix, value, true);
+	switch (v.type())
+	{
+		case cif::item_value_type::BOOLEAN:
+			os << std::boolalpha << v.m_data.m_value.m_boolean;
+			break;
+		case cif::item_value_type::INT:
+			os << v.m_data.m_value.m_integer;
+			break;
+		case cif::item_value_type::FLOAT:
+			os << v.m_data.m_value.m_float;
+			break;
+		case cif::item_value_type::TEXT:
+			os << v.m_data.sv();
+			break;
+		case cif::item_value_type::MISSING:
+			os << '?';
+			break;
+		case cif::item_value_type::EMPTY:
+			os << '.';
+			break;
+	}
+
+	return os;
 }
 
-void item_handle::swap(item_handle &b)
-{
-	assert(m_item_ix == b.m_item_ix);
-	// assert(&m_row_handle.m_category == &b.m_row_handle.m_category);
-	m_row_handle.swap(m_item_ix, b.m_row_handle);
-}
-
-}
+} // namespace cif

src/model.cpp

@@ -24,13 +24,17 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/model.hpp"
 #include "cif++.hpp"
+#include "cif++/point.hpp"
 
 #include <filesystem>
 #include <fstream>
+#include <initializer_list>
 #include <iomanip>
 #include <numeric>
 #include <stack>
+#include <stdexcept>
 
 namespace fs = std::filesystem;
 
@@ -58,7 +62,7 @@ void atom::atom_impl::moveTo(const point &p)
 
 std::string atom::atom_impl::get_property(std::string_view name) const
 {
-	return row()[name].as<std::string>();
+	return row()[name].get<std::string>();
 }
 
 int atom::atom_impl::get_property_int(std::string_view name) const
@@ -131,7 +135,7 @@ int atom::atom_impl::compare(const atom_impl &b) const
 
 int atom::atom_impl::get_charge() const
 {
-	auto formalCharge = row()["pdbx_formal_charge"].as<std::optional<int>>();
+	auto formalCharge = row()["pdbx_formal_charge"].get<std::optional<int>>();
 
 	if (not formalCharge.has_value())
 	{
@@ -348,17 +352,7 @@ std::tuple<point, float> residue::center_and_radius() const
 	for (auto &a : m_atoms)
 		pts.push_back(a.get_location());
 
-	auto center = centroid(pts);
-	float radius = 0;
-
-	for (auto &pt : pts)
-	{
-		float d = static_cast<float>(distance(pt, center));
-		if (radius < d)
-			radius = d;
-	}
-
-	return std::make_tuple(center, radius);
+	return smallest_sphere_around_points(pts);
 }
 
 bool residue::has_alternate_atoms() const
@@ -1275,28 +1269,28 @@ void structure::load_atoms_for_model(structure_open_options options)
 	else
 	{
 		std::vector<cif::mm::atom> atoms;
-		std::map<std::tuple<std::string,int>, std::map<std::string, float>> alts;
-	
+		std::map<std::tuple<std::string, int>, std::map<std::string, float>> alts;
+
 		for (auto id : atom_site.find<std::string>(std::move(c), "id"))
 		{
 			auto a = atoms.emplace_back(std::make_shared<atom::atom_impl>(m_db, id));
-	
+
 			if (a.is_alternate())
 			{
 				auto key = std::make_tuple(a.get_label_asym_id(), a.get_label_seq_id());
 				auto alt_id = a.get_label_alt_id();
-	
+
 				if (auto i = alts.find(key); i != alts.end())
 					i->second[alt_id] += a.get_occupancy();
 				else
 					alts[key][alt_id] = a.get_occupancy();
 			}
 		}
-	
+
 		for (auto &&[key, value] : alts)
 		{
 			// const auto &[asym_id, seq_id] = key;
-	
+
 			// select highest occupancy for this residue's alternates
 			std::string alt_id;
 			float occupancy = options.occupancy_mode == occupancy_policy::MAX ? 0.f : std::numeric_limits<float>::max();
@@ -1319,11 +1313,11 @@ void structure::load_atoms_for_model(structure_open_options options)
 					}
 				}
 			}
-	
+
 			value.clear();
 			value.emplace(alt_id, occupancy);
 		}
-	
+
 		for (auto a : atoms)
 		{
 			if (a.is_alternate())
@@ -1335,10 +1329,8 @@ void structure::load_atoms_for_model(structure_open_options options)
 			}
 			else
 				emplace_atom(a);
-	
 		}
 	}
-
 }
 
 void structure::load_data()
@@ -1898,13 +1890,8 @@ void structure::swap_atoms(atom a1, atom a2)
 		auto r1 = atomSites.find1(key("id") == a1.id());
 		auto r2 = atomSites.find1(key("id") == a2.id());
 
-		auto l1 = r1["label_atom_id"];
-		auto l2 = r2["label_atom_id"];
-		l1.swap(l2);
-
-		auto l3 = r1["auth_atom_id"];
-		auto l4 = r2["auth_atom_id"];
-		l3.swap(l4);
+		for (std::string fld : std::initializer_list<std::string>{ "label_atom_id", "auth_atom_id", "type_symbol" })
+			swap(r1[fld], r2[fld]);
 	}
 	catch (const std::exception &ex)
 	{
@@ -2348,6 +2335,36 @@ std::string structure::create_non_poly(const std::string &entity_id, std::vector
 	return asym_id;
 }
 
+std::string structure::create_non_poly(const std::string &compound_id, bool skip_hydrogen)
+{
+	auto compound = cif::compound_factory::instance().create(compound_id);
+	if (compound == nullptr)
+		throw std::runtime_error(std::format("{} is not a known compound", compound_id));
+	
+	std::vector<cif::row_initializer> atoms;
+	for (auto a : compound->atoms())
+	{
+		// We skip H-atoms, as fitting without H-atoms works better and we avoid conflicts in protonation states between CCD and MONLIB
+		if (skip_hydrogen and cif::atom_type_traits(a.type_symbol).symbol() == "H")
+			continue;
+
+		auto ax = a.get_location().get_x();
+		auto ay = a.get_location().get_y();
+		auto az = a.get_location().get_z();
+
+		atoms.emplace_back(cif::row_initializer{
+			{ "type_symbol", cif::atom_type_traits(a.type_symbol).symbol() },
+			{ "label_atom_id", a.id },
+			{ "auth_atom_id", a.id },
+			{ "Cartn_x", ax },
+			{ "Cartn_y", ay },
+			{ "Cartn_z", az },
+			{ "B_iso_or_equiv", 30.00 } });
+	}
+
+	return create_non_poly(create_non_poly_entity(compound_id), atoms);
+}
+
 void structure::create_water(row_initializer atom)
 {
 	using namespace literals;
@@ -2412,6 +2429,61 @@ void structure::create_water(row_initializer atom)
 	});
 }
 
+std::string structure::create_link(atom a1, atom a2, const std::string &link_type, const std::string &role)
+{
+	using namespace literals;
+
+	auto &struct_conn = m_db["struct_conn"];
+	auto &struct_conn_type = m_db["struct_conn_type"];
+
+	// This will validate link_type :-)
+	if (not struct_conn_type.contains("id"_key == link_type))
+		struct_conn_type.emplace({ { "id", link_type } });
+
+	std::string link_id = struct_conn.get_unique_id(link_type + '_');
+
+	item label_seq_id_1("ptnr1_label_seq_id");
+	if (int nr = a1.get_label_seq_id(); nr != 0)
+		label_seq_id_1.value(std::to_string(nr));
+
+	item label_seq_id_2("ptnr2_label_seq_id");
+	if (int nr = a2.get_label_seq_id(); nr != 0)
+		label_seq_id_2.value(std::to_string(nr));
+
+	struct_conn.emplace(
+		{ //
+			{ "id", link_id },
+			{ "conn_type_id", link_type },
+			{ "pdbx_leaving_atom_flag", "one" },
+
+			{ "ptnr1_label_asym_id", a1.get_label_asym_id() },
+			{ "ptnr1_label_comp_id", a1.get_label_comp_id() },
+			label_seq_id_1,
+			{ "ptnr1_label_atom_id", a1.get_label_atom_id() },
+			{ "pdbx_ptnr1_label_alt_id", a1.get_label_alt_id() },
+			{ "pdbx_ptnr1_PDB_ins_code", a1.get_pdb_ins_code() },
+			{ "ptnr1_auth_asym_id", a1.get_auth_asym_id() },
+			{ "ptnr1_auth_comp_id", a1.get_auth_comp_id() },
+			{ "ptnr1_auth_seq_id", a1.get_auth_seq_id() },
+			{ "ptnr1_symmetry", a1.symmetry() },
+
+			{ "ptnr2_label_asym_id", a2.get_label_asym_id() },
+			{ "ptnr2_label_comp_id", a2.get_label_comp_id() },
+			label_seq_id_2,
+			{ "ptnr2_label_atom_id", a2.get_label_atom_id() },
+			{ "pdbx_ptnr2_label_alt_id", a2.get_label_alt_id() },
+			{ "pdbx_ptnr2_PDB_ins_code", a2.get_pdb_ins_code() },
+			{ "ptnr2_auth_asym_id", a2.get_auth_asym_id() },
+			{ "ptnr2_auth_comp_id", a2.get_auth_comp_id() },
+			{ "ptnr2_auth_seq_id", a2.get_auth_seq_id() },
+			{ "ptnr2_symmetry", a2.symmetry() },
+
+			{ "pdbx_dist_value", distance(a1.get_location(), a2.get_location()), 3 },
+			{ "pdbx_role", role } });
+
+	return link_id;
+}
+
 branch &structure::create_branch()
 {
 	auto &entity = m_db["entity"];
@@ -2682,7 +2754,7 @@ void structure::cleanup_empty_categories()
 
 	for (auto chemComp : chem_comp)
 	{
-		std::string compID = chemComp["id"].as<std::string>();
+		std::string compID = chemComp["id"].get<std::string>();
 		if (atomSite.contains("label_comp_id"_key == compID or "auth_comp_id"_key == compID) or
 			pdbxPolySeqScheme.contains("mon_id"_key == compID or "auth_mon_id"_key == compID or "pdb_mon_id"_key == compID) or
 			entityPolySeq.contains("mon_id"_key == compID))
@@ -2703,7 +2775,7 @@ void structure::cleanup_empty_categories()
 
 	for (auto entity : entities)
 	{
-		std::string entityID = entity["id"].as<std::string>();
+		std::string entityID = entity["id"].get<std::string>();
 		if (atomSite.contains("label_entity_id"_key == entityID))
 			continue;
 
@@ -2845,8 +2917,8 @@ static int compare_numbers(std::string_view a, std::string_view b)
 
 	std::from_chars_result ra, rb;
 
-	ra = selected_charconv<double>::from_chars(a.data(), a.data() + a.length(), da);
-	rb = selected_charconv<double>::from_chars(b.data(), b.data() + b.length(), db);
+	ra = from_chars(a.data(), a.data() + a.length(), da);
+	rb = from_chars(b.data(), b.data() + b.length(), db);
 
 	if (not(bool) ra.ec and not(bool) rb.ec)
 	{

src/parser.cpp

@@ -24,15 +24,19 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include "cif++/utilities.hpp"
-#include "cif++/forward_decl.hpp"
 #include "cif++/parser.hpp"
+
 #include "cif++/file.hpp"
+#include "cif++/forward_decl.hpp"
+#include "cif++/item.hpp"
+#include "cif++/utilities.hpp"
 
 #include <cassert>
+#include <charconv>
 #include <iostream>
 #include <map>
 #include <stack>
+#include <system_error>
 
 namespace cif
 {
@@ -58,12 +62,12 @@ class reserved_words_automaton
 
 	constexpr bool finished() const
 	{
-		return m_state <= 0; 
+		return m_state <= 0;
 	}
 
 	constexpr bool matched() const
 	{
-		return m_state < 0; 
+		return m_state < 0;
 	}
 
 	constexpr move_result move(int ch)
@@ -75,7 +79,7 @@ class reserved_words_automaton
 			case 0:
 				break;
 
-			case -1:		// data_
+			case -1: // data_
 				if (sac_parser::is_non_blank(ch))
 					m_seen_trailing_chars = true;
 				else if (m_seen_trailing_chars)
@@ -84,15 +88,15 @@ class reserved_words_automaton
 					result = no_keyword;
 				break;
 
-			case -2:		// global_
+			case -2: // global_
 				result = sac_parser::is_non_blank(ch) ? no_keyword : global;
 				break;
 
-			case -3:		// loop_
+			case -3: // loop_
 				result = sac_parser::is_non_blank(ch) ? no_keyword : loop;
 				break;
 
-			case -4:		// save_
+			case -4: // save_
 				if (sac_parser::is_non_blank(ch))
 					m_seen_trailing_chars = true;
 				else if (m_seen_trailing_chars)
@@ -101,10 +105,10 @@ class reserved_words_automaton
 					result = save;
 				break;
 
-			case -5:		// stop_
+			case -5: // stop_
 				result = sac_parser::is_non_blank(ch) ? no_keyword : stop;
 				break;
-			
+
 			default:
 				assert(m_state > 0 and m_state < NODE_COUNT);
 
@@ -141,13 +145,13 @@ class reserved_words_automaton
 		int8_t next_nomatch;
 	} s_dag[] = {
 		{ 0 },
-		{ 'D',  5, 2 },
-		{ 'G',  9, 3 },
+		{ 'D', 5, 2 },
+		{ 'G', 9, 3 },
 		{ 'L', 15, 4 },
 		{ 'S', 19, 0 },
-		{ 'A',  6, 0 },
-		{ 'T',  7, 0 },
-		{ 'A',  8, 0 },
+		{ 'A', 6, 0 },
+		{ 'T', 7, 0 },
+		{ 'A', 8, 0 },
 		{ '_', -1, 0 },
 		{ 'L', 10, 0 },
 		{ 'O', 11, 0 },
@@ -155,7 +159,7 @@ class reserved_words_automaton
 		{ 'A', 13, 0 },
 		{ 'L', 14, 0 },
 		{ '_', -2, 0 },
-		{ 'O', 16, 0},
+		{ 'O', 16, 0 },
 		{ 'O', 17, 0 },
 		{ 'P', 18, 0 },
 		{ '_', -3, 0 },
@@ -238,7 +242,7 @@ int sac_parser::get_next_char()
 		}
 		else if (result == '\n')
 			++m_line_nr;
-		
+
 		m_token_buffer.push_back(std::char_traits<char>::to_char_type(result));
 	}
 
@@ -277,6 +281,8 @@ sac_parser::CIFToken sac_parser::get_next_token()
 	m_token_buffer.clear();
 	m_token_value = {};
 
+	bool negative = false;
+
 	reserved_words_automaton dag;
 
 	while (result == CIFToken::UNKNOWN)
@@ -310,6 +316,15 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				}
 				else if (dag.move(ch) == reserved_words_automaton::undefined)
 					state = State::Reserved;
+				else if (ch == '+' or ch == '-')
+				{
+					negative = true;
+					state = State::Numeric_Integer;
+				}
+				else if (ch >= '0' and ch <= '9')
+					state = State::Numeric_Integer;
+				else if (ch == '.')
+					state = State::Numeric_Float;
 				else
 					state = State::Value;
 				break;
@@ -326,7 +341,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else
 					m_bol = (ch == '\n');
 				break;
-			
+
 			case State::Comment:
 				if (ch == '\n')
 				{
@@ -339,12 +354,12 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else if (not is_any_print(ch))
 					error("invalid character in comment");
 				break;
-			
+
 			case State::QuestionMark:
 				if (not is_non_blank(ch))
 				{
 					retract();
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_UNKNOWN;
 				}
 				else
 					state = State::Value;
@@ -356,7 +371,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else if (ch == kEOF)
 					error("unterminated textfield");
 				else if (not is_any_print(ch) and cif::VERBOSE > 2)
-					warning("invalid character in text field '" + std::string({static_cast<char>(ch)}) + "' (" + std::to_string((int)ch) + ")");
+					warning("invalid character in text field '" + std::string({ static_cast<char>(ch) }) + "' (" + std::to_string((int)ch) + ")");
 				break;
 
 			case State::TextItemNL:
@@ -366,7 +381,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				{
 					assert(m_token_buffer.size() >= 2);
 					m_token_value = std::string_view(m_token_buffer.data() + 1, m_token_buffer.size() - 3);
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_TEXTFIELD;
 				}
 				else if (ch == kEOF)
 					error("unterminated textfield");
@@ -380,14 +395,14 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				else if (ch == quoteChar)
 					state = State::QuotedStringQuote;
 				else if (not is_any_print(ch) and cif::VERBOSE > 2)
-					warning("invalid character in quoted string: '" + std::string({static_cast<char>(ch)}) + "' (" + std::to_string((int)ch) + ")");
+					warning("invalid character in quoted string: '" + std::string({ static_cast<char>(ch) }) + "' (" + std::to_string((int)ch) + ")");
 				break;
 
 			case State::QuotedStringQuote:
 				if (is_white(ch))
 				{
 					retract();
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_CHARSTRING;
 					if (m_token_buffer.size() < 2)
 						error("Invalid quoted string token");
 
@@ -422,7 +437,7 @@ sac_parser::CIFToken sac_parser::get_next_token()
 						if (not is_non_blank(ch))
 						{
 							retract();
-							result = CIFToken::VALUE;
+							result = CIFToken::VALUE_CHARSTRING;
 							m_token_value = std::string_view(m_token_buffer.data(), m_token_buffer.size());
 						}
 						else
@@ -463,11 +478,65 @@ sac_parser::CIFToken sac_parser::get_next_token()
 				}
 				break;
 
+			case State::Numeric_Integer:
+				if (ch == '.')
+					state = State::Numeric_Float;
+				else if (ch == 'e' or ch == 'E')
+					state = State::Numeric_Exponent1;
+				else if (not is_non_blank(ch))
+				{
+					retract();
+					result = CIFToken::VALUE_NUMERIC_INTEGER;
+				}
+				else if (ch < '0' or ch > '9')
+					state = State::Value;
+				break;
+
+			case State::Numeric_Float:
+				if (not is_non_blank(ch))
+				{
+					retract();
+					if (m_token_buffer.size() == 1)
+						result = CIFToken::VALUE_INAPPLICABLE;
+					else
+						result = CIFToken::VALUE_NUMERIC_FLOAT;
+				}
+				else if (ch == 'e' or ch == 'E')
+					state = State::Numeric_Exponent1;
+				else if (ch < '0' or ch > '9')
+					state = State::Value;
+				break;
+
+			case State::Numeric_Exponent1:
+				if (ch == '+' or ch == '-' or (ch >= '0' and ch <= '9'))
+					state = State::Numeric_Exponent2;
+				else
+				{
+					if (VERBOSE > 0)
+						std::cerr << "parsing " << std::string_view{ m_token_buffer.data(), m_token_buffer.size() } << " Invalid floating point value, expected digit or sign character\n";
+					state = State::Value;
+				}
+				break;
+
+			case State::Numeric_Exponent2:
+				if (not is_non_blank(ch))
+				{
+					retract();
+					result = CIFToken::VALUE_NUMERIC_FLOAT;
+				}
+				else if (ch < '0' or ch > '9')
+				{
+					if (VERBOSE > 0)
+						std::cerr << "parsing " << std::string_view{ m_token_buffer.data(), m_token_buffer.size() } << " Invalid floating point value, expected exponent digit\n";
+					state = State::Value;
+				}
+				break;
+
 			case State::Value:
 				if (not is_non_blank(ch))
 				{
 					retract();
-					result = CIFToken::VALUE;
+					result = CIFToken::VALUE_CHARSTRING;
 					m_token_value = std::string_view(m_token_buffer.data(), m_token_buffer.size());
 					break;
 				}
@@ -480,12 +549,25 @@ sac_parser::CIFToken sac_parser::get_next_token()
 		}
 	}
 
-	if (VERBOSE >= 5)
+	// if (VERBOSE >= 5)
+	// {
+	// 	std::cerr << get_token_name(result);
+	// 	if (result != CIFToken::END_OF_FILE)
+	// 		std::cerr << " " << std::quoted(m_token_value);
+	// 	std::cerr << '\n';
+	// }
+
+	if (result == CIFToken::VALUE_NUMERIC_INTEGER)
 	{
-		std::cerr << get_token_name(result);
-		if (result != CIFToken::END_OF_FILE)
-			std::cerr << " " << std::quoted(m_token_value);
-		std::cerr << '\n';
+		auto [ptr, ec] = std::from_chars(m_token_buffer.data(), m_token_buffer.data() + m_token_buffer.size(), m_token_value_int);
+		if (ec != std::errc{})
+			error("Invalid integer value: " + std::make_error_code(ec).message());
+	}
+	else if (result == CIFToken::VALUE_NUMERIC_FLOAT)
+	{
+		auto [ptr, ec] = std::from_chars(m_token_buffer.data(), m_token_buffer.data() + m_token_buffer.size(), m_token_value_float);
+		if (ec != std::errc{})
+			error("Invalid integer value: " + std::make_error_code(ec).message());
 	}
 
 	return result;
@@ -661,7 +743,7 @@ sac_parser::datablock_index sac_parser::index_datablocks()
 			case data:
 				if (dblk[si] == 0 and is_non_blank(ch))
 				{
-					datablock = {static_cast<char>(ch)};
+					datablock = { static_cast<char>(ch) };
 					state = data_name;
 				}
 				else if (dblk[si++] != ch)
@@ -738,14 +820,17 @@ void sac_parser::parse_global()
 	while (m_lookahead == CIFToken::ITEM_NAME)
 	{
 		match(CIFToken::ITEM_NAME);
-		match(CIFToken::VALUE);
+		if (m_lookahead >= CIFToken::VALUE_INAPPLICABLE)
+			match(m_lookahead);
+		else
+			match(CIFToken::VALUE_CHARSTRING);
 	}
 }
 
 void sac_parser::parse_datablock()
 {
 	static const std::string kUnitializedCategory("<invalid>");
-	std::string cat = kUnitializedCategory;	// intial value acts as a guard for empty category names
+	std::string cat = kUnitializedCategory; // intial value acts as a guard for empty category names
 
 	while (m_lookahead == CIFToken::LOOP or m_lookahead == CIFToken::ITEM_NAME or m_lookahead == CIFToken::SAVE_NAME)
 	{
@@ -777,14 +862,38 @@ void sac_parser::parse_datablock()
 					match(CIFToken::ITEM_NAME);
 				}
 
-				while (m_lookahead == CIFToken::VALUE)
+				while (m_lookahead >= CIFToken::VALUE_INAPPLICABLE)
 				{
 					produce_row();
 
 					for (auto item_name : item_names)
 					{
-						produce_item(cat, item_name, m_token_value);
-						match(CIFToken::VALUE);
+						switch (m_lookahead)
+						{
+							case CIFToken::VALUE_INAPPLICABLE:
+								produce_item(cat, item_name, nullptr);
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_UNKNOWN:
+								produce_item(cat, item_name, std::optional<std::string>{});
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_NUMERIC_INTEGER:
+								produce_item(cat, item_name, m_token_value_int);
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_NUMERIC_FLOAT:
+								produce_item(cat, item_name, m_token_value_float);
+								match(m_lookahead);
+								break;
+							case CIFToken::VALUE_CHARSTRING:
+							case CIFToken::VALUE_TEXTFIELD:
+								produce_item(cat, item_name, m_token_value);
+								match(m_lookahead);
+								break;
+							default:;
+								match(CIFToken::VALUE_CHARSTRING);
+						}
 					}
 				}
 
@@ -806,9 +915,33 @@ void sac_parser::parse_datablock()
 
 				match(CIFToken::ITEM_NAME);
 
-				produce_item(cat, itemName, m_token_value);
+				switch (m_lookahead)
+				{
+					case CIFToken::VALUE_INAPPLICABLE:
+						produce_item(cat, itemName, nullptr);
+						match(CIFToken::VALUE_INAPPLICABLE);
+						break;
+					case CIFToken::VALUE_UNKNOWN:
+						produce_item(cat, itemName, item_value{ std::optional<std::string>{} });
+						match(CIFToken::VALUE_UNKNOWN);
+						break;
+					case CIFToken::VALUE_NUMERIC_INTEGER:
+						produce_item(cat, itemName, m_token_value_int);
+						match(CIFToken::VALUE_NUMERIC_INTEGER);
+						break;
+					case CIFToken::VALUE_NUMERIC_FLOAT:
+						produce_item(cat, itemName, m_token_value_float);
+						match(CIFToken::VALUE_NUMERIC_FLOAT);
+						break;
+					case CIFToken::VALUE_CHARSTRING:
+					case CIFToken::VALUE_TEXTFIELD:
+						produce_item(cat, itemName, m_token_value);
+						match(m_lookahead);
+						break;
+					default:
+						match(CIFToken::VALUE_CHARSTRING);
+				}
 
-				match(CIFToken::VALUE);
 				break;
 			}
 
@@ -864,7 +997,7 @@ void parser::produce_row()
 	// m_row.lineNr(m_line_nr);
 }
 
-void parser::produce_item(std::string_view category, std::string_view item, std::string_view value)
+void parser::produce_item(std::string_view category, std::string_view item, item_value value)
 {
 	if (VERBOSE >= 4)
 		std::cerr << "producing _" << category << '.' << item << " -> " << value << '\n';
@@ -872,7 +1005,7 @@ void parser::produce_item(std::string_view category, std::string_view item, std:
 	if (m_category == nullptr or not iequals(category, m_category->name()))
 		error("inconsistent categories in loop_");
 
-	m_row[item] = m_token_value;
+	m_row[item] = std::move(value);
 }
 
 } // namespace cif

src/pdb/pdb2cif.cpp

@@ -3142,7 +3142,6 @@ void PDBFileParser::ParseRemark350()
 	std::map<std::string, std::string> values;
 	std::vector<std::string> asymIdList;
 	std::smatch m;
-	cif::row_handle genR;
 
 	std::vector<double> mat, vec;
 

src/pdb/reconstruct.cpp

@@ -25,8 +25,12 @@
  */
 
 #include "cif++.hpp"
+#include "cif++/compound.hpp"
 #include "cif++/row.hpp"
 
+#include <stdexcept>
+#include <string>
+
 // --------------------------------------------------------------------
 
 namespace cif::pdb
@@ -186,6 +190,8 @@ void checkEntities(datablock &db)
 
 void createEntityIDs(datablock &db)
 {
+	using namespace literals;
+
 	// Suppose the file does not have entity ID's. We have to make up some
 
 	// walk the atoms. For each auth_asym_id we have a new struct_asym.
@@ -197,28 +203,44 @@ void createEntityIDs(datablock &db)
 	// that should cover it
 
 	auto &atom_site = db["atom_site"];
+	auto &entity = db["entity"];
 	auto &cf = compound_factory::instance();
 
-	std::vector<std::vector<residue_key_type>> entities;
+	std::vector<std::vector<row_handle>> entities;
 
 	std::string lastAsymID;
 	int lastSeqID = -1;
-	std::vector<residue_key_type> waters;
+	std::vector<row_handle> waters;
 
-	for (residue_key_type k : atom_site.rows<std::optional<std::string>,
-			 std::optional<int>,
-			 std::optional<std::string>,
-			 std::optional<std::string>,
-			 std::optional<int>,
-			 std::optional<std::string>>(
-			 "auth_asym_id", "auth_seq_id", "auth_comp_id",
-			 "label_asym_id", "label_seq_id", "label_comp_id"))
+	int nextEntityID;
+	try
 	{
+		if (entity.empty())
+			nextEntityID = 1;
+		else
+			nextEntityID = entity.find_max<int>("id") + 1;
+	}
+	catch (...)
+	{
+		nextEntityID = 1;
+	}
+
+	for (auto rh : atom_site.find("label_entity_id"_key == cif::null))
+	{
+		residue_key_type k = rh.get<std::optional<std::string>,
+			std::optional<int>,
+			std::optional<std::string>,
+			std::optional<std::string>,
+			std::optional<int>,
+			std::optional<std::string>>(
+			"auth_asym_id", "auth_seq_id", "auth_comp_id",
+			"label_asym_id", "label_seq_id", "label_comp_id");
+
 		std::string comp_id = get_comp_id(k);
 
 		if (cf.is_water(comp_id))
 		{
-			waters.emplace_back(k);
+			waters.emplace_back(rh);
 			continue;
 		}
 
@@ -227,19 +249,20 @@ void createEntityIDs(datablock &db)
 
 		bool is_monomer = cf.is_monomer(comp_id);
 
-		if (lastAsymID == asym_id and lastSeqID == seq_id and not is_monomer)
-			continue;
+		// if (lastAsymID == asym_id and lastSeqID == seq_id and not is_monomer)
+		// 	continue;
 
 		if (asym_id != lastAsymID or (not is_monomer and lastSeqID != seq_id))
 			entities.push_back({});
 
-		entities.back().emplace_back(k);
+		entities.back().emplace_back(rh);
 
 		lastAsymID = asym_id;
 		lastSeqID = seq_id;
 	}
 
 	std::map<std::size_t, std::string> entity_ids;
+	std::map<std::string, std::string> newEntitiesForCompound;
 
 	atom_site.add_item("label_entity_id");
 
@@ -248,7 +271,39 @@ void createEntityIDs(datablock &db)
 		if (entity_ids.contains(i))
 			continue;
 
-		auto entity_id = std::to_string(i + 1);
+		residue_key_type k = entities[i].front().get<std::optional<std::string>, std::optional<int>, std::optional<std::string>, std::optional<std::string>, std::optional<int>, std::optional<std::string>>(
+			"auth_asym_id", "auth_seq_id", "auth_comp_id",
+			"label_asym_id", "label_seq_id", "label_comp_id");
+
+		std::string comp_id = get_comp_id(k);
+
+		std::string entity_id;
+		if (auto v = db["pdbx_entity_nonpoly"].find_first("comp_id"_key == comp_id); v)
+			entity_id = v.get<std::string>("entity_id");
+		else if (auto i2 = newEntitiesForCompound.find(comp_id); i2 != newEntitiesForCompound.end())
+			entity_id = i2->second;
+		else
+		{
+			entity_id = std::to_string(nextEntityID++);
+
+			if (cf.is_monomer(comp_id))
+				entity.emplace({ //
+					{ "id", entity_id },
+					{ "type", "polymer" } });
+			else if (cf.is_water(comp_id))
+				entity.emplace({ //
+					{ "id", entity_id },
+					{ "type", "water" } });
+			else
+			{
+				entity.emplace({ //
+					{ "id", entity_id },
+					{ "type", "non-polymer" } });
+
+				newEntitiesForCompound[comp_id] = entity_id;
+			}
+		}
+
 		entity_ids[i] = entity_id;
 
 		for (std::size_t j = i + 1; j < entities.size(); ++j)
@@ -260,20 +315,17 @@ void createEntityIDs(datablock &db)
 
 	for (std::size_t ix = 0; auto &e : entities)
 	{
-		auto k = e.front();
 		const auto &entity_id = entity_ids[ix++];
 
-		std::string comp_id = get_comp_id(k);
-
-		for (auto &k2 : e)
-			atom_site.update_value(get_condition(k2), "label_entity_id", entity_id);
+		for (auto rh : e)
+			rh["label_entity_id"] = entity_id;
 	}
 
 	if (not waters.empty())
 	{
 		std::string waterEntityID = std::to_string(entities.size() + 1);
-		for (auto &k : waters)
-			atom_site.update_value(get_condition(k), "label_entity_id", waterEntityID);
+		for (auto rh : waters)
+			rh["label_entity_id"] = waterEntityID;
 	}
 }
 
@@ -440,9 +492,38 @@ void checkAtomRecords(datablock &db)
 	if (atom_site.contains(key("label_seq_id") < 0))
 		fixNegativeSeqID(atom_site);
 
-	std::set<int> polymer_entities;
-	for (int id : db["entity"].find<int>("type"_key == "polymer", "id"))
-		polymer_entities.insert(id);
+	std::set<std::string> polymer_entities;
+	if (db["entity"].empty())
+	{
+		// No entity, so we have to guess the types based on the content of atom_site
+
+		std::string last_entity_id;
+		std::optional<int> last_label_seq_id, last_auth_seq_id;
+
+		std::set<std::string> entityIDs;
+		for (auto &[entity_id, label_comp_id, label_seq_id, auth_comp_id, auth_seq_id] :
+			atom_site.rows<std::string, std::string, std::optional<int>, std::string, std::optional<int>>(
+				"label_entity_id", "label_comp_id", "label_seq_id", "auth_comp_id", "auth_seq_id"))
+		{
+			if (cf.is_water(label_comp_id) or cf.is_water(auth_comp_id))
+				continue;
+
+			if (polymer_entities.contains(entity_id))
+				continue;
+
+			if (last_entity_id == entity_id and (last_label_seq_id != label_seq_id or last_auth_seq_id != auth_seq_id))
+				polymer_entities.emplace(entity_id);
+
+			last_entity_id = entity_id;
+			last_label_seq_id = label_seq_id;
+			last_auth_seq_id = auth_seq_id;
+		}
+	}
+	else
+	{
+		for (std::string id : db["entity"].find<std::string>("type"_key == "polymer", "id"))
+			polymer_entities.insert(id);
+	}
 
 	std::set<std::string> missingCompounds;
 
@@ -479,7 +560,7 @@ void checkAtomRecords(datablock &db)
 		if (missingCompounds.contains(comp_id))
 			continue;
 
-		bool is_polymer = polymer_entities.contains(row["label_entity_id"].as<int>());
+		bool is_polymer = polymer_entities.contains(row["label_entity_id"].as<std::string>());
 		auto compound = cf.create(comp_id);
 
 		if (not compound)
@@ -533,14 +614,20 @@ void checkAtomRecords(datablock &db)
 		if (is_polymer and row["label_seq_id"].empty() and cf.is_monomer(comp_id))
 			row["label_seq_id"] = std::to_string(seq_id);
 
-		if (row["label_atom_id"].empty())
-			row["label_atom_id"] = row["auth_atom_id"].text();
 		if (row["label_asym_id"].empty())
 			row["label_asym_id"] = row["auth_asym_id"].text();
+		else if (row["auth_asym_id"].empty())
+			row["auth_asym_id"] = row["label_asym_id"].text();
+
 		if (row["label_comp_id"].empty())
 			row["label_comp_id"] = row["auth_comp_id"].text();
+		else if (row["auth_comp_id"].empty())
+			row["auth_comp_id"] = row["label_comp_id"].text();
+
 		if (row["label_atom_id"].empty())
 			row["label_atom_id"] = row["auth_atom_id"].text();
+		else if (row["auth_atom_id"].empty())
+			row["auth_atom_id"] = row["label_atom_id"].text();
 
 		// Rewrite the coordinates and other items that look better in a fixed format
 		// Be careful not to nuke invalidly formatted data here
@@ -563,7 +650,7 @@ void checkAtomRecords(datablock &db)
 			{
 				char b[12];
 
-				if (auto [ptr, ec] = cif::to_chars(b, b + sizeof(b), v, cif::chars_format::fixed, prec); ec == std::errc{})
+				if (auto [ptr, ec] = std::to_chars(b, b + sizeof(b), v, std::chars_format::fixed, prec); ec == std::errc{})
 					row.assign(item_name, { b, static_cast<std::string::size_type>(ptr - b) }, false, false);
 			}
 		}
@@ -605,19 +692,24 @@ void checkAtomAnisotropRecords(datablock &db)
 
 	std::vector<row_handle> to_be_deleted;
 
+	std::map<int, row_handle> atoms;
+	for (auto rh : atom_site)
+		atoms[rh.get<int>("id")] = rh;
+
 	bool warnReplaceTypeSymbol = true;
 	for (auto row : atom_site_anisotrop)
 	{
-		auto parents = atom_site_anisotrop.get_parents(row, atom_site);
-		if (parents.size() != 1)
+		auto ai = atoms.find(row.get<int>("id"));
+
+		if (ai == atoms.end())
 		{
 			to_be_deleted.emplace_back(row);
 			continue;
 		}
 
-		// this happens sometimes (Phenix):
+		auto parent = ai->second;
 
-		auto parent = parents.front();
+		// this happens sometimes (Phenix):
 
 		if (row["type_symbol"].empty())
 			row["type_symbol"] = parent["type_symbol"].text();
@@ -638,8 +730,6 @@ void checkAtomAnisotropRecords(datablock &db)
 			row["pdbx_label_atom_id"] = parent["label_atom_id"].text();
 		if (row["pdbx_label_comp_id"].empty() and not parent["label_comp_id"].empty())
 			row["pdbx_label_comp_id"] = parent["label_comp_id"].text();
-		// if (row["pdbx_PDB_model_num"].empty() and not parent["pdbx_PDB_model_num"].empty())
-		// 	row["pdbx_PDB_model_num"] = parent["pdbx_PDB_model_num"].text();
 	}
 
 	if (not to_be_deleted.empty())
@@ -652,23 +742,53 @@ void checkAtomAnisotropRecords(datablock &db)
 	}
 }
 
-void createStructAsym(datablock &db)
+void checkStructAsym(datablock &db)
 {
 	auto &atom_site = db["atom_site"];
 	auto &struct_asym = db["struct_asym"];
 
-	for (const auto &[label_asym_id, entity_id] : atom_site.rows<std::string, std::string>("label_asym_id", "label_entity_id"))
+	if (struct_asym.empty())
 	{
-		if (label_asym_id.empty())
-			throw std::runtime_error("File contains atom_site records without a label_asym_id");
-		if (struct_asym.count(key("id") == label_asym_id) == 0)
+		for (const auto &[label_asym_id, entity_id] : atom_site.rows<std::string, std::string>("label_asym_id", "label_entity_id"))
 		{
-			struct_asym.emplace({
-				// clang-format off
-				{ "id", label_asym_id },
-				{ "entity_id", entity_id }
-				//clang-format on
-			});
+			if (label_asym_id.empty())
+				throw std::runtime_error("File contains atom_site records without a label_asym_id");
+			if (struct_asym.count(key("id") == label_asym_id) == 0)
+			{
+				struct_asym.emplace({
+					// clang-format off
+					{ "id", label_asym_id },
+					{ "entity_id", entity_id }
+					//clang-format on
+				});
+			}
+		}
+	}
+	else
+	{
+		for (const auto &[label_asym_id, entity_id] :
+			atom_site.rows<std::string, std::string>("label_asym_id", "label_entity_id"))
+		{
+			if (label_asym_id.empty())
+				throw std::runtime_error("File contains atom_site records without a label_asym_id");
+			
+			auto sa = struct_asym.find_first(key("id") == label_asym_id);
+			if (sa)
+			{
+				if (sa["entity_id"].empty())
+					sa.assign("entity_id", entity_id, false, true);
+				else if (sa.get<std::string>("entity_id") != entity_id)
+					throw std::runtime_error("Inconsistent entity ID's in struct_asym");
+			}
+			else
+			{
+				struct_asym.emplace({
+					// clang-format off
+					{ "id", label_asym_id },
+					{ "entity_id", entity_id }
+					//clang-format on
+				});
+			}
 		}
 	}
 }
@@ -739,8 +859,11 @@ void createEntity(datablock &db)
 			auto c = cf.create(first_comp_id);
 
 			type = "non-polymer";
-			desc = c->name();
-			weight = c->formula_weight();
+			if (c)
+			{
+				desc = c->name();
+				weight = c->formula_weight();
+			}
 		}
 		else
 		{
@@ -811,28 +934,28 @@ void createEntityPoly(datablock &db)
 			if (type != "other")
 			{
 				std::string c_type;
-				if (cf.is_base(comp_id))
+				if (auto i = compound_factory::kBaseMap.find(comp_id); i != compound_factory::kBaseMap.end())
 				{
 					c_type = "polydeoxyribonucleotide";
-					letter_can = compound_factory::kBaseMap.at(comp_id);
+
+					letter_can = i->second;
+
 					if (comp_id.length() == 1)
 						letter = letter_can;
 					else
-						letter = '(' + letter_can + ')';
+						letter = '(' + comp_id + ')';
 				}
-				else if (cf.is_peptide(comp_id))
+				else if (auto i2 = compound_factory::kAAMap.find(comp_id); i2 != compound_factory::kAAMap.end())
 				{
 					c_type = "polypeptide(L)";
-					letter = letter_can = compound_factory::kAAMap.at(comp_id);
+
+					letter = letter_can = i2->second;
 				}
 				else if (iequals(c->type(), "D-PEPTIDE LINKING"))
 				{
 					c_type = "polypeptide(D)";
 
 					letter_can = c->one_letter_code();
-					if (letter_can == 0)
-						letter_can = 'X';
-
 					letter = '(' + comp_id + ')';
 
 					non_std_linkage = true;
@@ -843,9 +966,6 @@ void createEntityPoly(datablock &db)
 					c_type = "polypeptide(L)";
 
 					letter_can = c->one_letter_code();
-					if (letter_can == 0)
-						letter_can = 'X';
-
 					letter = '(' + comp_id + ')';
 
 					non_std_monomer = true;
@@ -855,9 +975,6 @@ void createEntityPoly(datablock &db)
 					// c_type = "other";
 
 					letter_can = c->one_letter_code();
-					if (letter_can == 0)
-						letter_can = 'X';
-
 					letter = '(' + comp_id + ')';
 
 					non_std_monomer = true;
@@ -870,7 +987,7 @@ void createEntityPoly(datablock &db)
 			}
 
 			seq[auth_asym_id] += letter;
-			seq_can[auth_asym_id] += letter_can;
+			seq_can[auth_asym_id] += letter_can ? letter_can : 'X';
 
 			if (find(pdb_strand_ids.begin(), pdb_strand_ids.end(), auth_asym_id) == pdb_strand_ids.end())
 				pdb_strand_ids.emplace_back(auth_asym_id);
@@ -1390,7 +1507,7 @@ bool reconstruct_pdbx(file &file, const validator &validator)
 	checkChemCompRecords(db);
 
 	// If the data is really horrible, it might not contain entities
-	if (not db["atom_site"].find_first(key("label_entity_id") != null))
+	if (db["atom_site"].find_first(key("label_entity_id") == null))
 		createEntityIDs(db);
 
 	// Now see if atom records make sense at all
@@ -1534,9 +1651,8 @@ bool reconstruct_pdbx(file &file, const validator &validator)
 		checkAtomAnisotropRecords(db);
 
 	// Now create any missing categories
-	// Next make sure we have struct_asym records
-	if (auto cat = db.get("struct_asym"); cat == nullptr or cat->empty())
-		createStructAsym(db);
+	// Next make sure we have good struct_asym records
+	checkStructAsym(db);
 
 	if (auto cat = db.get("entity"); cat == nullptr or cat->empty())
 		createEntity(db);
@@ -1645,9 +1761,8 @@ void fixup_pdbx(file &file, const validator &validator)
 	db.set_validator(&validator);
 
 	// Now create any missing categories
-	// Next make sure we have struct_asym records
-	if (auto cat = db.get("struct_asym"); cat == nullptr or cat->empty())
-		createStructAsym(db);
+	// Next make sure we have good struct_asym records
+	checkStructAsym(db);
 
 	if (auto cat = db.get("entity"); cat == nullptr or cat->empty())
 		createEntity(db);

src/pdb/validate-pdbx.cpp

@@ -69,7 +69,7 @@ bool is_valid_pdbx_file(const file &file, const validator &v)
 {
 	std::error_code ec;
 	bool result = is_valid_pdbx_file(file, v, ec);
-	return result and not (bool)ec;
+	return result and not(bool) ec;
 }
 
 bool is_valid_pdbx_file(const file &file, std::error_code &ec)
@@ -82,7 +82,7 @@ bool is_valid_pdbx_file(const file &file, std::error_code &ec)
 		result = is_valid_pdbx_file(file, validator_factory::instance().get(*ac), ec);
 	else
 		result = is_valid_pdbx_file(file, validator_factory::instance().get("mmcif_pdbx.dic"), ec);
-	
+
 	return result;
 }
 
@@ -167,7 +167,7 @@ bool is_valid_pdbx_file(const file &file, const validator &validator, std::error
 
 			const auto entity_poly_type = entity_poly.find1<std::string>("entity_id"_key == entity_id, "type");
 
-			std::map<int,std::set<std::string>> mon_per_seq_id;
+			std::map<int, std::set<std::string>> mon_per_seq_id;
 
 			for (const auto &[num, mon_id, hetero] : entity_poly_seq.find<int, std::string, bool>("entity_id"_key == entity_id, "num", "mon_id", "hetero"))
 			{
@@ -202,28 +202,37 @@ bool is_valid_pdbx_file(const file &file, const validator &validator, std::error
 					throw std::runtime_error("Mismatch between the hetero flag in the poly seq schemes and the number residues per seq_id");
 			}
 
-			for (const auto &[seq_id, mon_ids] : mon_per_seq_id)
-			{
-				for (auto asym_id : struct_asym.find<std::string>("entity_id"_key == entity_id, "id"))
-				{
-					condition cond;
-					
-					for (auto mon_id : mon_ids)
-						cond = std::move(cond) or "label_comp_id"_key == mon_id;
+			// This code proved to take too much time ...
 
-					cond = "label_entity_id"_key == entity_id and
-						"label_asym_id"_key == asym_id and
-						"label_seq_id"_key == seq_id and not std::move(cond);
-					
-					if (atom_site.contains(std::move(cond)))
-						throw std::runtime_error("An atom_site record exists that has no parent in the poly seq scheme categories");
-				}
+			// for (const auto &[seq_id, mon_ids] : mon_per_seq_id)
+			// {
+			// 	for (auto asym_id : struct_asym.find<std::string>("entity_id"_key == entity_id, "id"))
+			// 	{
+			// 		condition cond;
+
+			// 		for (auto mon_id : mon_ids)
+			// 			cond = std::move(cond) or "label_comp_id"_key == mon_id;
+
+			// 		cond = "label_entity_id"_key == entity_id and
+			// 			"label_asym_id"_key == asym_id and
+			// 			"label_seq_id"_key == seq_id and not std::move(cond);
+
+			// 		if (atom_site.contains(std::move(cond)))
+			// 			throw std::runtime_error("An atom_site record exists that has no parent in the poly seq scheme categories");
+			// 	}
+			// }
+
+			// ... so we're using this instead, should be almost the same...
+
+			for (const auto &[comp_id, seq_id] :
+				atom_site.find<std::string, int>("label_entity_id"_key == entity_id, "label_comp_id", "label_seq_id"))
+			{
+				if (not mon_per_seq_id[seq_id].contains(comp_id))
+					throw std::runtime_error("An atom_site record exists that has no parent in the poly seq scheme categories");
 			}
 
 			auto &&[seq, seq_can] = entity_poly.find1<std::optional<std::string>, std::optional<std::string>>("entity_id"_key == entity_id,
 				"pdbx_seq_one_letter_code", "pdbx_seq_one_letter_code_can");
-			
-			std::string::const_iterator si, sci, se, sce;
 
 			auto seq_match = [&](bool can, std::string::const_iterator si, std::string::const_iterator se)
 			{
@@ -260,8 +269,8 @@ bool is_valid_pdbx_file(const file &file, const validator &validator, std::error
 							else
 								letter = '(' + comp_id + ')';
 						}
-						
-						if (iequals(std::string{si, si + letter.length()}, letter))
+
+						if (iequals(std::string{ si, si + letter.length() }, letter))
 						{
 							match = true;
 							si += letter.length();
@@ -285,7 +294,9 @@ bool is_valid_pdbx_file(const file &file, const validator &validator, std::error
 			}
 			else
 			{
-				seq->erase(std::remove_if(seq->begin(), seq->end(), [](char ch) { return std::isspace(ch); }), seq->end());
+				seq->erase(std::remove_if(seq->begin(), seq->end(), [](char ch)
+							   { return std::isspace(ch); }),
+					seq->end());
 
 				if (not seq_match(false, seq->begin(), seq->end()))
 					throw std::runtime_error("Sequences do not match for entity " + entity_id);
@@ -298,7 +309,9 @@ bool is_valid_pdbx_file(const file &file, const validator &validator, std::error
 			}
 			else
 			{
-				seq_can->erase(std::remove_if(seq_can->begin(), seq_can->end(), [](char ch) { return std::isspace(ch); }), seq_can->end());
+				seq_can->erase(std::remove_if(seq_can->begin(), seq_can->end(), [](char ch)
+								   { return std::isspace(ch); }),
+					seq_can->end());
 
 				if (not seq_match(true, seq_can->begin(), seq_can->end()))
 					throw std::runtime_error("Canonical sequences do not match for entity " + entity_id);
@@ -315,11 +328,10 @@ bool is_valid_pdbx_file(const file &file, const validator &validator, std::error
 		ec = make_error_code(validation_error::not_valid_pdbx);
 	}
 
-	if (not result and (bool)ec)
+	if (not result and (bool) ec)
 		ec = make_error_code(validation_error::not_valid_pdbx);
 
 	return result;
 }
 
 } // namespace cif::pdb
-  

src/point.cpp

@@ -25,17 +25,19 @@
  */
 
 #include "cif++/point.hpp"
-#include "cif++/matrix.hpp"
 
-#include <cassert>
-#include <random>
+#include "cif++/matrix.hpp" // for matrix_subtraction, matrix_cofactors
+
+#include <initializer_list>
+#include <random> // for uniform_real_distribution, normal_distri...
+#include <stdexcept>
 
 namespace cif
 {
 
 // --------------------------------------------------------------------
 
-template<typename T>
+template <typename T>
 quaternion_type<T> normalize(quaternion_type<T> q)
 {
 	std::valarray<double> t(4);
@@ -126,10 +128,9 @@ quaternion construct_for_dihedral_angle(point p1, point p2, point p3, point p4,
 	p4 -= p3;
 	p3 -= p3;
 
-	quaternion q;
 	auto axis = -p2;
-
 	float dh = dihedral_angle(p1, p2, p3, p4);
+
 	return construct_from_angle_axis(angle - dh, axis);
 }
 
@@ -293,9 +294,9 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
 	}
 
 	quaternion q(
-		static_cast<float>(cf(maxR, 0)), 
-		static_cast<float>(cf(maxR, 1)), 
-		static_cast<float>(cf(maxR, 2)), 
+		static_cast<float>(cf(maxR, 0)),
+		static_cast<float>(cf(maxR, 1)),
+		static_cast<float>(cf(maxR, 2)),
 		static_cast<float>(cf(maxR, 3)));
 	q = normalize(q);
 
@@ -327,4 +328,251 @@ point nudge(point p, float offset)
 	return p + r;
 }
 
+// --------------------------------------------------------------------
+
+std::tuple<point, float> smallest_sphere_around_2_points(std::array<cif::point, 2> pts)
+{
+	return { (pts[0] + pts[1]) / 2, distance(pts[0], pts[1]) / 2 };
+}
+
+std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point, 3> pts)
+{
+	// Find two bisectors
+	auto vz = cross_product(pts[1] - pts[0], pts[2] - pts[0]);
+
+	auto bs1 = cross_product(vz, pts[1] - pts[0]);
+	bs1.normalize();
+
+	auto v1 = (pts[1] - pts[0]);
+	v1.normalize();
+
+	auto s1 = pts[0] + (distance(pts[1], pts[0]) / 2) * v1;
+
+	auto bs2 = cross_product(vz, pts[2] - pts[0]);
+	bs2.normalize();
+
+	auto v2 = (pts[2] - pts[0]);
+	v2.normalize();
+
+	auto s2 = pts[0] + (distance(pts[2], pts[0]) / 2) * v2;
+
+	auto c = line_line_intersection(s1, s1 + bs1, s2, s2 + bs2);
+	if (c)
+		return { *c, distance(*c, pts[0]) };
+
+	// Colinear points I guess, try something else
+	auto l1 = distance_squared(pts[0], pts[1]);
+	auto l2 = distance_squared(pts[0], pts[2]);
+	auto l3 = distance_squared(pts[1], pts[2]);
+
+	if (l1 > l2 and l1 > l3)
+		return smallest_sphere_around_2_points({ pts[0], pts[1] });
+	else if (l2 > l1 and l2 > l3)
+		return smallest_sphere_around_2_points({ pts[0], pts[2] });
+	else
+		return smallest_sphere_around_2_points({ pts[1], pts[2] });
+}
+
+std::tuple<point, float> smallest_sphere_around_4_points(std::array<cif::point, 4> pts)
+{
+	auto t0 = -norm_squared(pts[0]);
+	auto t1 = -norm_squared(pts[1]);
+	auto t2 = -norm_squared(pts[2]);
+	auto t3 = -norm_squared(pts[3]);
+
+	// clang-format off
+	matrix4x4<float> Tm({
+		pts[0].m_x, pts[0].m_y, pts[0].m_z, 1,
+		pts[1].m_x, pts[1].m_y, pts[1].m_z, 1,
+		pts[2].m_x, pts[2].m_y, pts[2].m_z, 1,
+		pts[3].m_x, pts[3].m_y, pts[3].m_z, 1
+	});
+	auto T = determinant(Tm);
+
+	if (T != 0)
+	{
+		matrix4x4<float> Dm({
+			t0, pts[0].m_y, pts[0].m_z, 1,
+			t1, pts[1].m_y, pts[1].m_z, 1,
+			t2, pts[2].m_y, pts[2].m_z, 1,
+			t3, pts[3].m_y, pts[3].m_z, 1
+		});
+		auto D = determinant(Dm) / T;
+		
+		matrix4x4<float> Em({
+			pts[0].m_x, t0, pts[0].m_z, 1,
+			pts[1].m_x, t1, pts[1].m_z, 1,
+			pts[2].m_x, t2, pts[2].m_z, 1,
+			pts[3].m_x, t3, pts[3].m_z, 1
+		});
+		auto E = determinant(Em) / T;
+		
+		matrix4x4<float> Fm({
+			pts[0].m_x, pts[0].m_y, t0, 1,
+			pts[1].m_x, pts[1].m_y, t1, 1,
+			pts[2].m_x, pts[2].m_y, t2, 1,
+			pts[3].m_x, pts[3].m_y, t3, 1
+		});
+		
+		auto F = determinant(Fm) / T;
+		
+		matrix4x4<float> Gm({
+			pts[0].m_x, pts[0].m_y, pts[0].m_z, t0,
+			pts[1].m_x, pts[1].m_y, pts[1].m_z, t1,
+			pts[2].m_x, pts[2].m_y, pts[2].m_z, t2,
+			pts[3].m_x, pts[3].m_y, pts[3].m_z, t3
+		});
+		auto G = determinant(Gm) / T;
+		
+		point center{ -D / 2, -E / 2, -F / 2 };
+		float radius = std::sqrt(D * D + E * E + F * F - 4 * G) / 2;
+
+		// clang-format on
+
+		return { center, radius };
+	}
+
+	// Perhaps some colinear points, try something else:
+
+	for (auto ix : std::initializer_list<std::array<size_t, 4>>{
+			 { 1, 2, 3, 0 },
+			 { 0, 2, 3, 1 },
+			 { 0, 1, 3, 2 },
+			 { 0, 1, 2, 3 },
+		 })
+	{
+		auto [center, radius] =
+			smallest_sphere_around_3_points({ pts[ix[0]], pts[ix[1]], pts[ix[2]] });
+
+		if (distance(pts[ix[3]], center) <= radius)
+			return { center, radius };
+	}
+
+	assert(false);
+	exit(1);
+}
+
+std::tuple<point, float> smallest_sphere_around_all_points(std::vector<point> P, std::vector<point> R)
+{
+	if (P.empty() or R.size() == 4)
+	{
+		switch (R.size())
+		{
+			case 1:
+				return { R[0], 0 };
+
+			case 2:
+				return smallest_sphere_around_2_points({ R[0], R[1] });
+
+			case 3:
+				return smallest_sphere_around_3_points({ R[0], R[1], R[2] });
+
+			case 4:
+				return smallest_sphere_around_4_points({ R[0], R[1], R[2], R[3] });
+
+			default:
+				assert(false);
+		}
+	}
+
+	auto p = P.back();
+	P.pop_back();
+
+	auto [c, r] = smallest_sphere_around_all_points(P, R);
+	assert(not std::isnan(r));
+	if (distance(c, p) <= r)
+		return { c, r };
+
+	R.emplace_back(p);
+	return smallest_sphere_around_all_points(P, R);
+}
+
+bool point_in_circle(point p, std::vector<point> c)
+{
+	switch (c.size())
+	{
+		case 0:
+			return false;
+
+		case 1:
+			return p == c.front();
+
+		case 2:
+		{
+			auto [center, radius] = smallest_sphere_around_2_points({ c[0], c[1] });
+			return cif::distance_squared(p, center) <= radius * radius;
+		}
+
+		case 3:
+		{
+			auto [center, radius] = smallest_sphere_around_3_points({ c[0], c[1], c[2] });
+			return cif::distance_squared(p, center) <= radius * radius;
+		}
+
+		case 4:
+		{
+			auto [center, radius] = smallest_sphere_around_4_points({ c[0], c[1], c[2], c[3] });
+			return cif::distance_squared(p, center) <= radius * radius;
+		}
+
+		default:
+			assert(false);
+			throw std::runtime_error("Error finding smallest sphere");
+	}
+}
+
+std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts)
+{
+	std::random_device rd;
+	std::mt19937 g(rd());
+
+	std::shuffle(pts.begin(), pts.end(), g);
+
+	std::vector<size_t> cix;
+
+	auto cirle_points = [&]()
+	{
+		std::vector<point> result;
+		for (auto ix : cix)
+			result.emplace_back(pts[ix]);
+		return result;
+	};
+
+	size_t i = 0;
+	while (i < pts.size())
+	{
+		if (std::find(cix.begin(), cix.end(), i) != cix.end() or
+			point_in_circle(pts[i], cirle_points()))
+		{
+			++i;
+		}
+		else
+		{
+			cix.erase(std::remove_if(cix.begin(), cix.end(), [i](size_t j)
+						  { return j < i; }),
+				cix.end());
+			cix.push_back(i);
+			if (cix.size() < 4)
+				i = 0;
+			else
+				++i;
+		}
+	}
+
+	switch (cix.size())
+	{
+		case 1:
+			return { pts[cix[0]], 0 };
+		case 2:
+			return smallest_sphere_around_2_points({ pts[cix[0]], pts[cix[1]] });
+		case 3:
+			return smallest_sphere_around_3_points({ pts[cix[0]], pts[cix[1]], pts[cix[2]] });
+		case 4:
+			return smallest_sphere_around_4_points({ pts[cix[0]], pts[cix[1]], pts[cix[2]], pts[cix[3]] });
+		default:
+			assert(false);
+			throw std::runtime_error("Error finding smallest sphere");
+	}
+}
+
 } // namespace cif

src/row.cpp

@@ -24,17 +24,50 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/row.hpp"
+
 #include "cif++/category.hpp"
+#include "cif++/item.hpp"
+#include <stdexcept>
 
 namespace cif
 {
+// --------------------------------------------------------------------
 
-void row_handle::assign(uint16_t item, std::string_view value, bool updateLinked, bool validate)
+item_value s_null_item;
+
+item_value &row_handle::operator[](uint16_t item_ix)
+{
+	return empty() or item_ix >= m_row->size() ? s_null_item : m_row->operator[](item_ix);
+}
+
+const item_value &row_handle::operator[](uint16_t item_ix) const
+{
+	return empty() or item_ix >= m_row->size() ? s_null_item : m_row->operator[](item_ix);
+}
+
+item_value &row_handle::operator[](std::string_view item_name)
+{
+	auto ix = add_item(item_name);
+	if (ix >= size())
+		m_row->resize(ix + 1);
+	
+	return m_row->operator[](ix);
+}
+
+const item_value &row_handle::operator[](std::string_view item_name) const
+{
+	return operator[](get_item_ix(item_name));
+}
+
+// --------------------------------------------------------------------
+
+void row_handle::assign(uint16_t item, item_value value, bool updateLinked, bool validate)
 {
 	if (not m_category)
 		throw std::runtime_error("uninitialized row");
 
-	m_category->update_value(m_row, item, value, updateLinked, validate);
+	m_category->update_value(m_row, item, std::move(value), updateLinked, validate);
 }
 
 uint16_t row_handle::get_item_ix(std::string_view name) const
@@ -86,28 +119,29 @@ row_initializer::row_initializer(row_handle rh)
 		auto &i = r->operator[](ix);
 		if (not i)
 			continue;
-		emplace_back(cat.get_item_name(ix), i.text());
+		emplace_back(cat.get_item_name(ix), i);
 	}
 }
 
-void row_initializer::set_value(std::string_view name, std::string_view value)
+void row_initializer::set_value(std::string name, item_value value)
 {
 	for (auto &i : *this)
 	{
 		if (i.name() == name)
 		{
-			i.value(value);
+			i.value(std::move(value));
 			return;
 		}
 	}
 
-	emplace_back(name, value);
+	emplace_back(std::move(name), std::move(value));
 }
 
-void row_initializer::set_value_if_empty(std::string_view name, std::string_view value)
+void row_initializer::set_value_if_empty(std::string name, item_value value)
 {
-	if (find_if(begin(), end(), [name](auto &i) { return i.name() == name; }) == end())
-		emplace_back(name, value);
+	if (std::ranges::find_if(*this, [name](auto &i)
+			{ return i.name() == name; }) == end())
+		emplace_back(std::move(name), std::move(value));
 }
 
 } // namespace cif

src/text.cpp

@@ -28,6 +28,11 @@
 
 #include <algorithm>
 #include <cassert>
+#include <charconv>
+
+#if __has_include("fast_float/fast_float.h")
+#include "fast_float/fast_float.h"
+#endif
 
 namespace cif
 {
@@ -512,4 +517,32 @@ std::vector<std::string> word_wrap(const std::string &text, std::size_t width)
 	return result;
 }
 
+#if __has_include("fast_float/fast_float.h")
+
+template <typename T>
+std::from_chars_result ff_charconv<T, typename std::enable_if_t<std::is_floating_point_v<T>>>::from_chars(const char *a, const char *b, T &v)
+{
+	auto r = fast_float::from_chars(a, b, v);
+	return { r.ptr, r.ec };
+}
+
+template struct ff_charconv<float>;
+template struct ff_charconv<double>;
+// template struct ff_charconv<long double>;
+
+#ifdef __STDCPP_FLOAT64_T__
+template struct ff_charconv<std::float64_t>;
+#endif
+#ifdef __STDCPP_FLOAT32_T__
+template struct ff_charconv<std::float32_t>;
+#endif
+#ifdef __STDCPP_FLOAT16_T__
+template struct ff_charconv<std::float16_t>;
+#endif
+#ifdef __STDCPP_BFLOAT16_T__
+template struct ff_charconv<std::bfloat16_t>;
+#endif
+
+#endif
+
 } // namespace cif

src/utilities.cpp

@@ -34,14 +34,20 @@
 #include <condition_variable>
 #include <cstring>
 #include <deque>
+#include <format>
 #include <fstream>
-#include <functional>
 #include <iomanip>
 #include <iostream>
 #include <map>
 #include <mutex>
 #include <sstream>
+#include <string>
 #include <thread>
+#include <utility>
+
+#if __cpp_lib_jthread >= 201911L
+#include <stop_token>
+#endif
 
 namespace fs = std::filesystem;
 
@@ -65,27 +71,50 @@ std::string get_version_nr()
 
 #if defined(_WIN32) or defined(__MINGW32__)
 }
-#include <windows.h>
-#include <libloaderapi.h>
-#include <wincon.h>
-
-#include <codecvt>
+// clang-format off
+# include <windows.h>
+# include <libloaderapi.h>
+# include <wincon.h>
+// clang-format on
 
 namespace cif
 {
 
 uint32_t get_terminal_width()
 {
-    CONSOLE_SCREEN_BUFFER_INFO csbi;
-    ::GetConsoleScreenBufferInfo(::GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
-    return csbi.srWindow.Right - csbi.srWindow.Left + 1;
+	CONSOLE_SCREEN_BUFFER_INFO csbi;
+	return ::GetConsoleScreenBufferInfo(::GetStdHandle(STD_OUTPUT_HANDLE), &csbi)
+	           ? csbi.srWindow.Right - csbi.srWindow.Left + 1
+	           : 80;
+}
+
+void write_to_console(const std::string &s)
+{
+	auto h = ::GetStdHandle(STD_OUTPUT_HANDLE);
+	CONSOLE_SCREEN_BUFFER_INFO csbi;
+
+	if (auto l = ::MultiByteToWideChar(CP_UTF8, 0, s.data(), s.length(), nullptr, 0);
+		l > 0 and ::GetConsoleScreenBufferInfo(::GetStdHandle(STD_OUTPUT_HANDLE), &csbi))
+	{
+		std::u16string ws(l, 0);
+
+		::MultiByteToWideChar(CP_UTF8, 0, s.data(), s.length(), (LPWSTR)ws.data(), l);
+
+		DWORD w;
+		::WriteConsoleW(h, ws.data(), ws.length(), &w, nullptr);
+	}
+	else
+	{
+		std::cout.write(s.data(), s.length());
+		std::cout.flush();
+	}
 }
 
 #else
 
-#include <sys/ioctl.h>
-#include <termios.h>
-#include <limits.h>
+# include <limits.h>
+# include <sys/ioctl.h>
+# include <termios.h>
 
 uint32_t get_terminal_width()
 {
@@ -100,59 +129,220 @@ uint32_t get_terminal_width()
 	return result;
 }
 
+inline void write_to_console(const std::string &s)
+{
+	std::cout << s << std::flush;
+}
+
 #endif
 
 // --------------------------------------------------------------------
 
 struct progress_bar_impl
 {
-	progress_bar_impl(int64_t inMax, const std::string &inAction)
-		: m_max_value(inMax)
+	progress_bar_impl(uint64_t max_value, const std::string &message)
+		: m_max_value(max_value)
 		, m_consumed(0)
-		, m_action(inAction)
-		, m_message(inAction)
-		, m_thread(std::bind(&progress_bar_impl::run, this))
+		, m_action(message)
+		, m_message(message)
 	{
 	}
 
-	progress_bar_impl(const progress_bar_impl&) = delete;
-	progress_bar_impl &operator=(const progress_bar_impl &) = delete;
+	virtual ~progress_bar_impl() {}
 
-	~progress_bar_impl();
-
-	void run();
-
-	void consumed(int64_t n);
-	void progress(int64_t p);
-	void message(const std::string &msg);
-
-	void print_progress();
-	void print_done();
+	virtual void consumed(uint64_t n);
+	virtual void progress(uint64_t p);
+	virtual void message(const std::string &msg);
+	virtual void print_done();
 
 	using time_point = std::chrono::time_point<std::chrono::system_clock>;
 
-	int64_t m_max_value;
-	std::atomic<int64_t> m_consumed;
-	int64_t m_last_consumed = 0;
-	int m_spinner_index = 0;
+	uint64_t m_max_value;
+	std::atomic<uint64_t> m_consumed;
 	std::string m_action, m_message;
-	std::mutex m_mutex;
-	std::thread m_thread;
 	time_point m_start = std::chrono::system_clock::now();
-	time_point m_last = std::chrono::system_clock::now();
-	bool m_stop = false;
 };
 
-progress_bar_impl::~progress_bar_impl()
+void progress_bar_impl::consumed(uint64_t n)
+{
+	m_consumed += n;
+}
+
+void progress_bar_impl::progress(uint64_t p)
+{
+	m_consumed = p;
+}
+
+void progress_bar_impl::message(const std::string &msg)
+{
+	m_message = msg;
+}
+
+void progress_bar_impl::print_done()
+{
+	std::chrono::duration<double> elapsed = std::chrono::system_clock::now() - m_start;
+	std::string days, hours, minutes, seconds;
+
+	uint64_t s = static_cast<uint64_t>(std::trunc(elapsed.count()));
+	if (s > 24 * 60 * 60)
+	{
+		days = std::format("{:d}d ", s / (24 * 60 * 60));
+		s %= 24 * 60 * 60;
+	}
+
+	if (s > 60 * 60)
+	{
+		hours = std::format("{:2d}h ", s / (60 * 60));
+		s %= 60 * 60;
+	}
+
+	if (s > 60)
+	{
+		minutes = std::format("{:2d}m ", s / 60);
+		s %= 60;
+	}
+
+	std::string msg = std::format("{} done in {}{}{}{:.1f}s", m_action, days, hours, minutes, s + 1e-6 * (elapsed.count() - s));
+
+	uint32_t width = get_terminal_width();
+
+	if (msg.length() < width)
+		msg += std::string(width - msg.length(), ' ');
+
+	write_to_console(msg += '\n');
+}
+
+// --------------------------------------------------------------------
+
+struct simple_progress_bar_impl : public progress_bar_impl
+{
+	simple_progress_bar_impl(uint64_t max_value, const std::string &message)
+		: progress_bar_impl(max_value, message)
+	{
+	}
+
+	~simple_progress_bar_impl()
+	{
+		if (m_printed_any)
+			print_done();
+	}
+
+	void consumed(uint64_t n) override
+	{
+		using namespace std::literals;
+
+		progress_bar_impl::consumed(n);
+
+		// print at most 10 steps, but only if it took long enough
+
+		int percentile = static_cast<int>(std::floor(10.f * m_consumed / m_max_value));
+		if (percentile > m_last_percentile and (m_printed_any or std::chrono::system_clock::now() - m_start >= 1s))
+		{
+			if (not std::exchange(m_printed_any, true))
+				write_to_console(m_action + ": ");
+
+			write_to_console(std::format("...{:d}0%", percentile));
+			m_last_percentile = percentile;
+		}
+	}
+
+	void progress(uint64_t p) override
+	{
+		consumed(p - m_consumed);
+	}
+
+	void print_done() override
+	{
+		if (m_printed_any)
+		{
+			write_to_console("\n");
+			progress_bar_impl::print_done();
+		}
+	}
+
+	bool m_printed_any = false;
+	int m_last_percentile = 0;
+};
+
+// --------------------------------------------------------------------
+
+struct fancy_progress_bar_impl : public progress_bar_impl
+{
+	fancy_progress_bar_impl(uint64_t max_value, const std::string &message)
+		: progress_bar_impl(max_value, message)
+		, m_thread(
+#if __cpp_lib_jthread >= 201911L
+			  [this](std::stop_token stoken)
+			  { this->run(stoken); }
+#else
+			  [this]()
+			  { this->run(); }
+#endif
+		  )
+	{
+	}
+
+	~fancy_progress_bar_impl();
+
+#if __cpp_lib_jthread >= 201911L
+	void run(std::stop_token stoken);
+#else
+	void run();
+#endif
+
+	void consumed(uint64_t n) override;
+	void progress(uint64_t p) override;
+	void message(const std::string &msg) override;
+
+	void print_progress();
+
+	std::mutex m_mutex;
+	std::condition_variable m_cv;
+
+	float m_progress;
+	uint32_t m_width, m_bar_width;
+	uint32_t m_steps, m_last_steps = 0;
+	uint64_t m_last_consumed = 0;
+#if __cpp_lib_jthread >= 201911L
+	std::jthread m_thread;
+#else
+	std::thread m_thread;
+	bool m_stop = false;
+#endif
+};
+
+const char *kBlocks[] = {
+	" ",
+	"▏",
+	"▎",
+	"▍",
+	"▌",
+	"▋",
+	"▊",
+	"▉",
+	"█",
+};
+
+const size_t kBlockCount = sizeof(kBlocks) / sizeof(void *) - 1;
+
+fancy_progress_bar_impl::~fancy_progress_bar_impl()
 {
 	using namespace std::literals;
 	assert(m_thread.joinable());
 
+#if __cpp_lib_jthread >= 201911L
+	m_thread.request_stop();
+#else
 	m_stop = true;
+#endif
 	m_thread.join();
 }
 
-void progress_bar_impl::run()
+#if __cpp_lib_jthread >= 201911L
+void fancy_progress_bar_impl::run(std::stop_token stoken)
+#else
+void fancy_progress_bar_impl::run()
+#endif
 {
 	using namespace std::literals;
 
@@ -160,25 +350,44 @@ void progress_bar_impl::run()
 
 	try
 	{
-		while (not m_stop)
+		for (;;)
 		{
+			std::unique_lock lock(m_mutex);
+
+			m_cv.wait_for(lock, 25ms);
+
+#if __cpp_lib_jthread >= 201911L
+			if (stoken.stop_requested())
+				break;
+#else
+			if (m_stop)
+				break;
+#endif
+
 			auto now = std::chrono::system_clock::now();
 
-			if (now - m_start < 2s or now - m_last < 100ms)
-			{
-				std::this_thread::sleep_for(10ms);
+			if (m_consumed == m_last_consumed or now - m_start < 1s)
 				continue;
-			}
 
-			std::lock_guard lock(m_mutex);
+			m_last_consumed = m_consumed;
 
-			if (not printedAny and isatty(STDOUT_FILENO))
-				std::cout << "\x1b[?25l";
+			// See if we need to do work
+			m_width = get_terminal_width();
+			m_progress = static_cast<float>(m_consumed) / m_max_value;
+			m_bar_width = 7 * m_width / 10; // 70% of the width of the terminal
+			m_steps = static_cast<uint32_t>(std::ceil(m_progress * m_bar_width * kBlockCount));
+
+			if (m_steps == m_last_steps)
+				continue;
+
+			m_last_steps = m_steps;
+
+			if (not printedAny)
+				write_to_console("\x1b[?25l");
 
 			print_progress();
 
 			printedAny = true;
-			m_last = std::chrono::system_clock::now();
 		}
 	}
 	catch (...)
@@ -187,161 +396,94 @@ void progress_bar_impl::run()
 
 	if (printedAny)
 	{
+		write_to_console("\r\x1b[?25h");
 		print_done();
-		if (isatty(STDOUT_FILENO))
-			std::cout << "\x1b[?25h";
 	}
 }
 
-void progress_bar_impl::consumed(int64_t n)
+void fancy_progress_bar_impl::consumed(uint64_t n)
 {
-	m_consumed += n;
+	progress_bar_impl::consumed(n);
+	// m_cv.notify_one();
 }
 
-void progress_bar_impl::progress(int64_t p)
+void fancy_progress_bar_impl::progress(uint64_t p)
 {
-	m_consumed = p;
+	progress_bar_impl::progress(p);
+	// m_cv.notify_one();
 }
 
-void progress_bar_impl::message(const std::string &msg)
+void fancy_progress_bar_impl::message(const std::string &msg)
 {
 	std::unique_lock lock(m_mutex);
-	m_message = msg;
+	progress_bar_impl::message(msg);
+	// m_cv.notify_one();
 }
 
-const char* kSpinner[] = {
-	// ".", "o", "O", "0", "O", "o", ".", " "
-	// "⢄", "⢂", "⢁", "⡁", "⡈", "⡐", "⡠"
-	 ".", "o", "O", "0", "@", "*", " "
-};
-
-const std::size_t kSpinnerCount = sizeof(kSpinner) / sizeof(char*);
-
-const int kSpinnerTimeInterval = 100;
-
 const uint32_t kMinBarWidth = 40, kMinMsgWidth = 12;
 
-void progress_bar_impl::print_progress()
+void fancy_progress_bar_impl::print_progress()
 {
-	const char *kBlocks[] = {
-		// "▯", // 0
-		// "▮", // 1
-		"=",
-		"-"
-	};
+	const uint32_t pct_width = 5;
+	uint32_t msg_width = m_width - m_bar_width - pct_width - 1;
 
-	uint32_t width = get_terminal_width();
-
-	float progress = static_cast<float>(m_consumed) / m_max_value;
-	
-	if (width < kMinBarWidth)
-		std::cout << (100 * progress) << "%\n";
-	else
+	if (msg_width < kMinMsgWidth)
 	{
-		uint32_t bar_width = 7 * width / 10;
-		uint32_t pct_width = 7;
-		uint32_t msg_width = width - bar_width - pct_width - 1;
+		m_bar_width += kMinMsgWidth - msg_width;
+		msg_width = kMinMsgWidth;
+	}
 
-		if (msg_width < kMinMsgWidth)
-		{
-			bar_width += kMinMsgWidth - msg_width;
-			msg_width = kMinMsgWidth;
-		}
+	std::string bar;
+	bar.reserve(m_bar_width * 4);
 
-		std::ostringstream msg;
-
-		if (m_message.length() <= msg_width)
-		{
-			msg << m_message;
-			if (m_message.length() < msg_width)
-				msg << std::string(msg_width - m_message.length(), ' ');
-		}
+	for (uint32_t i = 0; i < m_bar_width; ++i)
+	{
+		if (i * kBlockCount <= m_steps)
+			bar += kBlocks[kBlockCount];
+		else if (i * kBlockCount > m_steps + kBlockCount)
+			bar += kBlocks[0];
 		else
-			msg << m_message.substr(0, msg_width - 3) << "...";
-
-		msg << ' ';
-
-		uint32_t pi = static_cast<uint32_t>(std::ceil(progress * bar_width));
-
-		for (uint32_t i = 0; i < bar_width; ++i)
-			msg << kBlocks[i > pi ? 1 : 0];
-
-		msg << ' ';
-
-		msg << std::setw(3) << static_cast<int>(std::ceil(progress * 100)) << "% ";
-
-		auto now = std::chrono::system_clock::now();
-		m_spinner_index = (std::chrono::duration_cast<std::chrono::milliseconds>(now - m_start).count() / kSpinnerTimeInterval) % kSpinnerCount;
-
-		msg << kSpinner[m_spinner_index];
-
-		std::cout << '\r' << msg.str();
-		std::cout.flush();
+			bar += kBlocks[1 + m_steps % kBlockCount];
 	}
-}
 
-namespace
-{
-
-	std::ostream &operator<<(std::ostream &os, const std::chrono::duration<double> &t)
+	// make the bar more colorfull
+	struct color_type
 	{
-		uint64_t s = static_cast<uint64_t>(std::trunc(t.count()));
-		if (s > 24 * 60 * 60)
-		{
-			auto days = s / (24 * 60 * 60);
-			os << days << "d ";
-			s %= 24 * 60 * 60;
-		}
+		uint8_t r, g, b;
+	} fg{ 0, 3, 5 }, bg{ 0, 1, 2 };
 
-		if (s > 60 * 60)
-		{
-			auto hours = s / (60 * 60);
-			os << hours << "h ";
-			s %= 60 * 60;
-		}
+	auto esc_1 = std::format("\x1b[38;5;{}m\x1b[48;5;{}m",
+		16 + (fg.r * 36) + (fg.g * 6) + fg.b,
+		16 + (bg.r * 36) + (bg.g * 6) + bg.b);
+	std::string esc_2("\x1b[0m");
 
-		if (s > 60)
-		{
-			auto minutes = s / 60;
-			os << minutes << "m ";
-			s %= 60;
-		}
+	bar = esc_1 + bar + esc_2;
 
-		double ss = s + 1e-6 * (t.count() - s);
+	std::string msg = m_message.length() <= msg_width
+	                      ? m_message
+	                      : m_message.substr(0, msg_width - 3) + "...";
 
-		os << std::fixed << std::setprecision(1) << ss << 's';
-
-		return os;
-	}
-
-} // namespace
-
-void progress_bar_impl::print_done()
-{
-	std::chrono::duration<double> elapsed = std::chrono::system_clock::now() - m_start;
-
-	std::ostringstream msgstr;
-	msgstr << m_action << " done in " << elapsed << " seconds";
-	auto msg = msgstr.str();
-
-	uint32_t width = get_terminal_width();
-
-	if (msg.length() < width)
-		msg += std::string(width - msg.length(), ' ');
-
-	std::cout << '\r' << msg << '\n';
+	write_to_console(std::format("{:{}} {} {:3d}%\r", msg, msg_width, bar,
+		static_cast<int>(std::ceil(m_progress * 100))));
 }
 
-progress_bar::progress_bar(int64_t inMax, const std::string &inAction)
+// --------------------------------------------------------------------
+
+progress_bar::progress_bar(int64_t max_value, const std::string &message)
 	: m_impl(nullptr)
 {
-	if (isatty(STDOUT_FILENO) and VERBOSE >= 0)
-		m_impl = new progress_bar_impl(inMax, inAction);
+	if (VERBOSE >= 0)
+	{
+		if (isatty(STDOUT_FILENO) and get_terminal_width() > kMinBarWidth)
+			m_impl = new fancy_progress_bar_impl(max_value, message);
+		else
+			m_impl = new simple_progress_bar_impl(max_value, message);
+	}
 }
 
 progress_bar::~progress_bar()
 {
-	delete m_impl;
+	flush();
 }
 
 void progress_bar::consumed(int64_t inConsumed)
@@ -350,16 +492,25 @@ void progress_bar::consumed(int64_t inConsumed)
 		m_impl->consumed(inConsumed);
 }
 
-void progress_bar::progress(int64_t inProgress)
+void progress_bar::progress(int64_t value)
 {
 	if (m_impl != nullptr)
-		m_impl->progress(inProgress);
+		m_impl->progress(value);
 }
 
-void progress_bar::message(const std::string &inMessage)
+void progress_bar::message(const std::string &message)
 {
 	if (m_impl != nullptr)
-		m_impl->message(inMessage);
+		m_impl->message(message);
+}
+
+void progress_bar::flush()
+{
+	if (m_impl)
+	{
+		delete m_impl;
+		m_impl = nullptr;
+	}
 }
 
 } // namespace cif
@@ -387,13 +538,13 @@ struct rsrc_imp
 
 #if _WIN32
 
-#if __MINGW32__
+# if __MINGW32__
 
 extern "C" __attribute__((weak, alias("gResourceIndexDefault"))) const mrsrc::rsrc_imp gResourceIndex[];
 extern "C" __attribute__((weak, alias("gResourceDataDefault"))) const char gResourceData[];
 extern "C" __attribute__((weak, alias("gResourceNameDefault"))) const char gResourceName[];
 
-#else
+# else
 
 extern "C" const mrsrc::rsrc_imp *gResourceIndexDefault[1] = {};
 extern "C" const char *gResourceDataDefault[1] = {};
@@ -403,11 +554,11 @@ extern "C" const mrsrc::rsrc_imp gResourceIndex[];
 extern "C" const char gResourceData[];
 extern "C" const char gResourceName[];
 
-#pragma comment(linker, "/alternatename:gResourceIndex=gResourceIndexDefault")
-#pragma comment(linker, "/alternatename:gResourceData=gResourceDataDefault")
-#pragma comment(linker, "/alternatename:gResourceName=gResourceNameDefault")
+#  pragma comment(linker, "/alternatename:gResourceIndex=gResourceIndexDefault")
+#  pragma comment(linker, "/alternatename:gResourceData=gResourceDataDefault")
+#  pragma comment(linker, "/alternatename:gResourceName=gResourceNameDefault")
 
-#endif
+# endif
 
 #else
 extern const __attribute__((weak)) mrsrc::rsrc_imp gResourceIndex[];

src/validate.cpp

@@ -27,12 +27,11 @@
 #include "cif++/validate.hpp"
 #include "cif++/category.hpp"
 #include "cif++/dictionary_parser.hpp"
-#include "cif++/gzio.hpp"
 #include "cif++/utilities.hpp"
 
 #include <cassert>
-#include <fstream>
 #include <iostream>
+#include <mutex>
 
 // The validator depends on regular expressions. Unfortunately,
 // the implementation of std::regex in g++ is buggy and crashes
@@ -75,6 +74,7 @@ struct regex_impl
   private:
 	pcre2_code *m_rx = nullptr;
 	pcre2_match_data *m_data = nullptr;
+	mutable std::mutex m_mutex;
 };
 
 regex_impl::regex_impl(std::string_view rx)
@@ -95,6 +95,8 @@ regex_impl::regex_impl(std::string_view rx)
 
 regex_impl::~regex_impl()
 {
+	std::unique_lock lock(m_mutex);
+
 	if (m_data)
 		pcre2_match_data_free(m_data);
 
@@ -104,6 +106,8 @@ regex_impl::~regex_impl()
 
 bool regex_impl::match(std::string_view v) const
 {
+	std::unique_lock lock(m_mutex);
+
 	bool result = false;
 
 	if (int rc = pcre2_match(m_rx, (PCRE2_SPTR)v.data(), v.length(), 0, 0, m_data, nullptr); rc >= 0)
@@ -181,8 +185,8 @@ int type_validator::compare(std::string_view a, std::string_view b) const
 
 				std::from_chars_result ra, rb;
 
-				ra = selected_charconv<double>::from_chars(a.data(), a.data() + a.length(), da);
-				rb = selected_charconv<double>::from_chars(b.data(), b.data() + b.length(), db);
+				ra = from_chars(a.data(), a.data() + a.length(), da);
+				rb = from_chars(b.data(), b.data() + b.length(), db);
 
 				if (not(bool) ra.ec and not(bool) rb.ec)
 				{
@@ -259,24 +263,24 @@ int type_validator::compare(std::string_view a, std::string_view b) const
 
 // --------------------------------------------------------------------
 
-void item_validator::operator()(std::string_view value) const
+void item_validator::operator()(const item_value &value) const
 {
-	std::error_code ec;
-	if (not validate_value(value, ec))
-		throw std::system_error(ec, std::string{ value } + " does not match rx for " + m_item_name);
+	// std::error_code ec;
+	// if (not validate_value(value, ec))
+	// 	throw std::system_error(ec, std::string{ value } + " does not match rx for " + m_item_name);
 }
 
-bool item_validator::validate_value(std::string_view value, std::error_code &ec) const noexcept
+bool item_validator::validate_value(const item_value &value, std::error_code &ec) const noexcept
 {
 	ec.clear();
 
-	if (not value.empty() and value != "?" and value != ".")
-	{
-		if (m_type != nullptr and not m_type->m_rx->match(value))
-			ec = make_error_code(validation_error::value_does_not_match_rx);
-		else if (not m_enums.empty() and m_enums.count(std::string{ value }) == 0)
-			ec = make_error_code(validation_error::value_is_not_in_enumeration_list);
-	}
+	// if (not value.empty() and value != "?" and value != ".")
+	// {
+	// 	if (m_type != nullptr and not m_type->m_rx->match(value))
+	// 		ec = make_error_code(validation_error::value_does_not_match_rx);
+	// 	else if (not m_enums.empty() and m_enums.count(std::string{ value }) == 0)
+	// 		ec = make_error_code(validation_error::value_is_not_in_enumeration_list);
+	// }
 
 	return not(bool) ec;
 }
@@ -561,8 +565,10 @@ const validator &validator_factory::get(const category &audit_conform)
 	// If the audit conform contains only one record, this is easy
 	if (audit_conform.size() == 1)
 	{
-		const auto &[name, version] = audit_conform.front().get<std::string, std::optional<std::string>>("dict_name", "dict_version");
-		return m_validators.emplace_back(construct_validator(name, version));
+		const auto &[name, version] =
+			audit_conform.front().get<std::string, std::optional<std::string>>("dict_name", "dict_version");
+		if (not name.empty())
+			return m_validators.emplace_back(construct_validator(name, version));
 	}
 
 	// A new, merged dictionary
@@ -570,6 +576,9 @@ const validator &validator_factory::get(const category &audit_conform)
 	std::optional<validator> v;
 	for (const auto &[name, version] : audit_conform.rows<std::string, std::optional<std::string>>("dict_name", "dict_version"))
 	{
+		if (name.empty())
+			continue;
+
 		if (not v)
 			v = construct_validator(name, version);
 		else
@@ -604,7 +613,7 @@ validator validator_factory::construct_validator(std::string_view name, std::opt
 		not v.matches_audit_conform(category{ "audit_conform", //
 			{ { "dict_name", name }, { "dict_version", version } } }))
 	{
-		std::clog << "Invalid dictionary?\n";
+		std::clog << "Loaded dictionary does not match name=" << name << " and version=" << version.value_or("''") << "\n";
 	}
 
 	return v;

test/CMakeLists.txt

@@ -1,19 +1,39 @@
-# We're using the older version 2 of Catch2
+# SPDX-License-Identifier: BSD-2-Clause
+# 
+# Copyright (c) 2025  NKI/AVL, Netherlands Cancer Institute
+# 
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# 
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+# 
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-if(NOT(Catch2_FOUND OR TARGET Catch2))
-	find_package(Catch2 QUIET)
+if(NOT (Catch2_FOUND OR TARGET Catch2))
+	find_package(Catch2 3 QUIET)
 
 	if(NOT Catch2_FOUND)
-		include(FetchContent)
-
 		FetchContent_Declare(
 			Catch2
 			GIT_REPOSITORY https://github.com/catchorg/Catch2.git
-			GIT_TAG v2.13.9)
+			GIT_TAG v3.4.0)
 
 		FetchContent_MakeAvailable(Catch2)
 
-		set(Catch2_VERSION "2.13.9")
+		target_compile_features(Catch2 PRIVATE cxx_std_20)
 	endif()
 endif()
 
@@ -21,26 +41,21 @@ list(
 	APPEND
 	CIFPP_tests
 	unit-v2
-	unit-3d
-	model
-	query
-	rename-compound
-	sugar
-	spinner
+	# unit-3d
+	# model
+	# query
+	# rename-compound
+	# sugar
+	# spinner
 	# reconstruction
-	validate-pdbx
-	)
+	# validate-pdbx
+	# matrix
+)
 
 add_library(test-main OBJECT "${CMAKE_CURRENT_SOURCE_DIR}/test-main.cpp")
 
 target_link_libraries(test-main cifpp::cifpp Catch2::Catch2)
 
-if("${Catch2_VERSION}" VERSION_LESS 3.0.0)
-	target_compile_definitions(test-main PUBLIC CATCH22=1)
-else()
-	target_compile_definitions(test-main PUBLIC CATCH22=0)
-endif()
-
 foreach(CIFPP_TEST IN LISTS CIFPP_tests)
 	set(CIFPP_TEST "${CIFPP_TEST}-test")
 	set(CIFPP_TEST_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/${CIFPP_TEST}.cpp")
@@ -48,12 +63,6 @@ foreach(CIFPP_TEST IN LISTS CIFPP_tests)
 	add_executable(
 		${CIFPP_TEST} ${CIFPP_TEST_SOURCE} $<TARGET_OBJECTS:test-main>)
 
-	if(${Catch2_VERSION} VERSION_GREATER_EQUAL 3.0.0)
-		target_compile_definitions(${CIFPP_TEST} PUBLIC CATCH22=0)
-	else()
-		target_compile_definitions(${CIFPP_TEST} PUBLIC CATCH22=1)
-	endif()
-
 	target_link_libraries(${CIFPP_TEST} PRIVATE cifpp::cifpp Catch2::Catch2)
 	target_include_directories(${CIFPP_TEST} PRIVATE "${EIGEN_INCLUDE_DIR}")
 
@@ -62,15 +71,8 @@ foreach(CIFPP_TEST IN LISTS CIFPP_tests)
 		target_compile_options(${CIFPP_TEST} PRIVATE /EHsc)
 	endif()
 
-	add_custom_target(
-		"run-${CIFPP_TEST}"
-		DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/Run${CIFPP_TEST}.touch ${CIFPP_TEST})
-
-	add_custom_command(
-		OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/Run${CIFPP_TEST}.touch
-		COMMAND $<TARGET_FILE:${CIFPP_TEST}> --data-dir
-		${CMAKE_CURRENT_SOURCE_DIR})
-
-	add_test(NAME ${CIFPP_TEST} COMMAND $<TARGET_FILE:${CIFPP_TEST}> --data-dir
-		${CMAKE_CURRENT_SOURCE_DIR})
-endforeach()
+	if(NOT (CIFPP_TEST STREQUAL "spinner-test"))
+		add_test(NAME ${CIFPP_TEST}
+			COMMAND $<TARGET_FILE:${CIFPP_TEST}> --data-dir ${CMAKE_CURRENT_SOURCE_DIR})
+	endif()
+endforeach()

test/matrix-test.cpp

@@ -0,0 +1,103 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2025 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "cif++/matrix.hpp"
+#include "test-main.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+#include <cif++.hpp>
+
+TEST_CASE("m1")
+{
+	cif::matrix3x3<int> m = cif::identity_matrix<int>(3);
+
+	CHECK(cif::determinant(m) == 1);
+}
+
+TEST_CASE("m2")
+{
+	cif::matrix4x4<int> m = cif::identity_matrix<int>(4);
+
+	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
+	CHECK(ms == cif::identity_matrix<int>(3));
+}
+
+TEST_CASE("m3")
+{
+	cif::matrix4x4<int> m{
+		{ 1, 2, 3, 4,      //
+			5, 6, 7, 8,    //
+			9, 10, 11, 12, //
+			13, 14, 15, 16 }
+	};
+	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
+
+	cif::matrix3x3<int> t{
+		{ 1, 3, 4, 9, 11, 12, 13, 15, 16 }
+	};
+
+	CHECK(ms == t);
+}
+
+TEST_CASE("m4")
+{
+	cif::matrix4x4<int> m{
+		{
+			-2,
+			3,
+			1,
+			0,
+			4,
+			1,
+			-3,
+			2,
+			0,
+			-1,
+			2,
+			5,
+			3,
+			2,
+			0,
+			-4,
+		}
+	};
+
+	// std::cout << m << "\n\n";
+
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0)) << "\n\n";
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1)) << "\n\n";
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2)) << "\n\n";
+	// std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3)) << "\n\n";
+
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0))) << "\n\n";
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1))) << "\n\n";
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2))) << "\n\n";
+	// std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3))) << "\n\n";
+
+	CHECK(cif::determinant(m) == 332);
+}
+
+

test/query-test.cpp

@@ -49,6 +49,4 @@ TEST_CASE("q-1")
 
 		CHECK(pdbx_poly_seq_scheme.count("asym_id"_key == "A" and "entity_id"_key == 1 and "seq_id"_key == 1 and "mon_id"_key == "PRO" and "hetero"_key == false) == 1);
 	}
-	
-
 }

test/reconstruction-test.cpp

@@ -28,6 +28,7 @@
 
 #include <cif++.hpp>
 
+#include <filesystem>
 #include <iostream>
 #include <fstream>
 

test/test-main.cpp

@@ -2,7 +2,7 @@
 
 #include "test-main.hpp"
 
-#include <cif++.hpp>
+#include <cif++/utilities.hpp>
 
 std::filesystem::path gTestDir = std::filesystem::current_path();
 
@@ -11,12 +11,7 @@ int main(int argc, char *argv[])
 	Catch::Session session; // There must be exactly one instance
 
 	// Build a new parser on top of Catch2's
-#if CATCH22
-	using namespace Catch::clara;
-#else
-	// Build a new parser on top of Catch2's
 	using namespace Catch::Clara;
-#endif
 
 	auto cli = session.cli()                               // Get Catch2's command line parser
 	           | Opt(gTestDir, "data-dir")                 // bind variable to a new option, with a hint string
@@ -38,7 +33,7 @@ int main(int argc, char *argv[])
 	// initialize CCD location
 	cif::add_file_resource("components.cif", gTestDir / ".." / "rsrc" / "ccd-subset.cif");
 
-	cif::compound_factory::instance().push_dictionary(gTestDir / "HEM.cif");
+	// cif::compound_factory::instance().push_dictionary(gTestDir / "HEM.cif");
 
 	return session.run();
 }

test/test-main.hpp

@@ -26,11 +26,7 @@
 
 #pragma once
 
-#if CATCH22
-#include <catch2/catch.hpp>
-#else
 #include <catch2/catch_all.hpp>
-#endif
 
 #include <filesystem>
 

test/unit-3d-test.cpp

@@ -24,15 +24,17 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "cif++/point.hpp"
 #include "test-main.hpp"
 
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <cif++.hpp>
 #include <stdexcept>
 
-#include <cif++.hpp>
-
 #if defined(_MSC_VER)
-#pragma warning (disable : 5054)	// warning C5054: operator '&': deprecated between enumerations of different types
-#pragma warning (disable : 4127)	// conditional expression is constant
+# pragma warning(disable : 5054) // warning C5054: operator '&': deprecated between enumerations of different types
+# pragma warning(disable : 4127) // conditional expression is constant
 #endif
 
 #include <Eigen/Eigenvalues>
@@ -315,8 +317,7 @@ TEST_CASE("m2q_0a")
 		cif::point p2 = p1;
 		p2.rotate(q);
 
-		cif::matrix3x3<float> rot_c({
-			static_cast<float>(d[0]),
+		cif::matrix3x3<float> rot_c({ static_cast<float>(d[0]),
 			static_cast<float>(d[1]),
 			static_cast<float>(d[2]),
 			static_cast<float>(d[3]),
@@ -324,8 +325,7 @@ TEST_CASE("m2q_0a")
 			static_cast<float>(d[5]),
 			static_cast<float>(d[6]),
 			static_cast<float>(d[7]),
-			static_cast<float>(d[8])
-		});
+			static_cast<float>(d[8]) });
 
 		cif::point p3 = rot_c * p1;
 
@@ -610,3 +610,40 @@ TEST_CASE("volume_3bwh_1")
 
 	CHECK_THAT(c.get_cell().get_volume(), Catch::Matchers::WithinRel(741009.625f, 0.01f));
 }
+
+// --------------------------------------------------------------------
+
+TEST_CASE("smallest_sphere-1")
+{
+	std::vector<cif::point> pts{
+		{ 0.9295, 4.9006, 46.9706 },
+		{ -0.1215, 5.5936, 46.0726 },
+		{ -0.7975, 4.7046, 45.0796 },
+		{ -1.4875, 3.5486, 45.7196 },
+		{ -0.6535, 2.8816, 46.8186 },
+		{ 0.3825, 3.5156, 47.4496 },
+		{ 1.1995, 2.9206, 48.5286 },
+		{ 0.8255, 2.0466, 49.4716 },
+		{ 1.6625, 1.5036, 50.5176 },
+		{ 1.1165, 0.6056, 51.3626 },
+		{ 1.8325, -0.0064, 52.4656 },
+		{ 1.1945, -0.9044, 53.2216 },
+		{ 1.8135, -1.5534, 54.3566 },
+		{ 1.0925, -2.4574, 55.0656 },
+		{ 1.5205, -3.2204, 56.2476 },
+		{ 1.1955, 5.8066, 48.1796 },
+		{ 2.2495, 4.6896, 46.1796 },
+		{ -1.2515, 1.5186, 47.1786 },
+		{ 3.1385, 1.9106, 50.6166 },
+		{ 3.2605, -1.1834, 54.7206 },
+		{ 2.5975, -3.8554, 56.2096 },
+		{ 0.7975, -3.2184, 57.2686 }
+	};
+
+	for (int i = 0; i < 1000; ++i)
+	{
+		auto [c, r] = cif::smallest_sphere_around_points(pts);
+		CHECK_THAT(cif::distance(c, cif::point{ 0, 0.743099928, 51.1741028 }), Catch::Matchers::WithinAbs(0.f, 0.01f));
+		CHECK_THAT(r, Catch::Matchers::WithinAbs(7.31248331f, 0.01f));
+	}
+}

tools/update-libcifpp-data.in

@@ -63,7 +63,7 @@ update_dictionary() {
 
 update_dictionary "@CIFPP_CACHE_DIR@/components.cif" "https://files.wwpdb.org/pub/pdb/data/monomers/components.cif.gz"
 update_dictionary "@CIFPP_CACHE_DIR@/mmcif_pdbx.dic" "https://mmcif.wwpdb.org/dictionaries/ascii/mmcif_pdbx_v50.dic.gz"
-update_dictionary "@CIFPP_CACHE_DIR@/mmcif_ma.dic" "https://github.com/ihmwg/ModelCIF/raw/master/dist/mmcif_ma.dic"
+update_dictionary "@CIFPP_CACHE_DIR@/mmcif_ma.dic" "https://mmcif.wwpdb.org/dictionaries/ascii/mmcif_ma.dic"
 
 # notify subscribers, using find instead of run-parts to make it work on FreeBSD as well