text:
#pragma once

#include <atomic>
#include <array>
#include <cstdint>
#include <future>
#include <functional>
#include <optional>
#include <bitset>

#ifndef SPDLOG_FMT_EXTERNAL
#define SPDLOG_FMT_EXTERNAL 1
#endif
#include <spdlog/spdlog.h>

#include <gsl/gsl-lite.hpp>

// TODO: leaf pulls in Windows.h via common.hpp, so we bypass that here to avoid the mess
// if something did pull in windows.h, this is fine. custom_formatmessage is just a wrapper.
#ifndef _WINDOWS_

#ifndef LPVOID
#define LPVOID void*
#define LPVOID_redefined
#endif
#ifndef LPCSTR
#define LPCSTR char*
#define LPCSTR_redefined
#endif

namespace mu
{
	auto custom_formatmessage(
		unsigned long dwFlags,
		const void*	  lpSource,
		unsigned long dwMessageId,
		unsigned long dwLanguageId,
		char*		  lpBuffer,
		unsigned long nSize,
		va_list*	  Arguments) noexcept -> unsigned long;
}

#define FormatMessageA(DWORD_dwFlags, LPCVOID_lpSource, DWORD_dwMessageId, DWORD_dwLanguageId, LPSTR_lpBuffer, DWORD_nSize, va_list_Arguments)                                     \
	mu::custom_formatmessage(DWORD_dwFlags, LPCVOID_lpSource, DWORD_dwMessageId, DWORD_dwLanguageId, LPSTR_lpBuffer, DWORD_nSize, va_list_Arguments)
#define _WINDOWS_
#define _WINDOWS_redefined			   true
#define FORMAT_MESSAGE_ALLOCATE_BUFFER 0x00000100
#define FORMAT_MESSAGE_FROM_SYSTEM	   0x00001000
#define FORMAT_MESSAGE_IGNORE_INSERTS  0x00000200
#define MAKELANGID(p, s)			   ((((unsigned short)(s)) << 10) | (unsigned short)(p))
#define LANG_NEUTRAL				   0x00
#define SUBLANG_DEFAULT				   0x01
#define LPSTR						   char*
#else // #ifndef _WINDOWS_
#define _WINDOWS_redefined false
#endif // #else // #ifndef _WINDOWS_
#include <boost/leaf.hpp>
#if _WINDOWS_redefined
#ifdef LPVOID_redefined
#undef LPVOID
#undef LPVOID_redefined
#endif
#ifdef LPCSTR_redefined
#undef LPCSTR_redefined
#undef LPCSTR
#endif
#undef FormatMessageA
#undef _WINDOWS_
#undef _WINDOWS_redefined
#undef FORMAT_MESSAGE_ALLOCATE_BUFFER
#undef FORMAT_MESSAGE_FROM_SYSTEM
#undef FORMAT_MESSAGE_IGNORE_INSERTS
#undef MAKELANGID
#undef LANG_NEUTRAL
#undef SUBLANG_DEFAULT
#undef LPSTR
#endif // #if _WINDOWS_redefined

namespace mu
{
	template<typename E>
	constexpr auto underlying_cast(E e) noexcept -> typename std::underlying_type<E>::type
	{
		return static_cast<typename std::underlying_type<E>::type>(e);
	}
} // namespace mu

namespace mu
{
	namespace leaf = boost::leaf;

	struct common_error
	{
	};

	struct runtime_error
	{
		struct not_specified : common_error
		{
		};
	};
} // namespace mu

#define MU_LEAF_IF_REPORT(r)                                                                                \
	auto&& BOOST_LEAF_TMP = r;                                                                                \
	static_assert(::boost::leaf::is_result_type<typename std::decay<decltype(BOOST_LEAF_TMP)>::type>::value, "MU_LEAF_CHECK requires a result object (see is_result_type)");       \
	if (!BOOST_LEAF_TMP)                                                                                \
		/*TODO: REPORT*/;                                                                                \
	else

#define MU_LEAF_RETHROW(r)                                                                                \
	auto&& BOOST_LEAF_TMP = r;                                                                                \
	static_assert(::boost::leaf::is_result_type<typename std::decay<decltype(BOOST_LEAF_TMP)>::type>::value, "MU_LEAF_CHECK requires a result object (see is_result_type)");       \
	if (BOOST_LEAF_TMP)                                                                                \
		;                                                                                \
	else                                                                                \
	{                                                                                \
		throw BOOST_LEAF_TMP.error();                                                                                \
	}

#define MU_LEAF_ASSIGN(v, r)                                                                                \
	auto&& BOOST_LEAF_TMP = r;                                                                                \
	static_assert(                                                                                \
		::boost::leaf::is_result_type<std::decay<decltype(BOOST_LEAF_TMP)>::type>::value,                                                                                \
		"MU_LEAF_ASSIGN and MU_LEAF_AUTO require a result object as the second argument (see is_result_type)");                                                                    \
	if (!BOOST_LEAF_TMP)                                                                                \
	{                                                                                \
		return BOOST_LEAF_TMP.error();                                                                                \
	}                                                                                \
	v = std::forward<decltype(BOOST_LEAF_TMP)>(BOOST_LEAF_TMP).value()

#define MU_LEAF_PUSH_BACK(v, r)                                                                                \
	auto&& BOOST_LEAF_TMP = r;                                                                                \
	static_assert(                                                                                \
		::boost::leaf::is_result_type<std::decay<decltype(BOOST_LEAF_TMP)>::type>::value,                                                                                \
		"MU_LEAF_ASSIGN and MU_LEAF_AUTO require a result object as the second argument (see is_result_type)");                                                                    \
	if (!BOOST_LEAF_TMP)                                                                                \
	{                                                                                \
		return BOOST_LEAF_TMP.error();                                                                                \
	}                                                                                \
	v.push_back(std::forward<decltype(BOOST_LEAF_TMP)>(BOOST_LEAF_TMP).value())

#define MU_LEAF_EMPLACE_BACK(v, r)                                                                                \
	auto&& BOOST_LEAF_TMP = r;                                                                                \
	static_assert(                                                                                \
		::boost::leaf::is_result_type<std::decay<decltype(BOOST_LEAF_TMP)>::type>::value,                                                                                \
		"MU_LEAF_ASSIGN and MU_LEAF_AUTO require a result object as the second argument (see is_result_type)");                                                                    \
	if (!BOOST_LEAF_TMP)                                                                                \
	{                                                                                \
		return BOOST_LEAF_TMP.error();                                                                                \
	}                                                                                \
	v.emplace_back(std::forward<decltype(BOOST_LEAF_TMP)>(BOOST_LEAF_TMP).value())

#define MU_LEAF_AUTO(v, r) MU_LEAF_ASSIGN(auto v, r)

#define MU_LEAF_ASSIGN_THROW(v, r)                                                                                \
	auto&& BOOST_LEAF_TMP = r;                                                                                \
	static_assert(                                                                                \
		::boost::leaf::is_result_type<std::decay<decltype(BOOST_LEAF_TMP)>::type>::value,                                                                                \
		"MU_LEAF_ASSIGN and MU_LEAF_AUTO require a result object as the second argument (see is_result_type)");                                                                    \
	if (!BOOST_LEAF_TMP)                                                                                \
	{                                                                                \
		throw BOOST_LEAF_TMP.error();                                                                                \
	}                                                                                \
	v = std::forward<decltype(BOOST_LEAF_TMP)>(BOOST_LEAF_TMP).value()

#define MU_LEAF_AUTO_THROW(v, r) MU_LEAF_ASSIGN_THROW(auto v, r)

#define MU_LEAF_CHECK(r)                                                                                \
	auto&& BOOST_LEAF_TMP = r;                                                                                \
	static_assert(::boost::leaf::is_result_type<std::decay<decltype(BOOST_LEAF_TMP)>::type>::value, "MU_LEAF_CHECK requires a result object (see is_result_type)");                \
	if (BOOST_LEAF_TMP)                                                                                \
		;                                                                                \
	else                                                                                \
	{                                                                                \
		return BOOST_LEAF_TMP.error();                                                                                \
	}

#define MU_LEAF_NEW_ERROR		::boost::leaf::leaf_detail::inject_loc{__FILE__, __LINE__, __FUNCTION__} + ::boost::leaf::new_error
#define MU_LEAF_EXCEPTION		::boost::leaf::leaf_detail::inject_loc{__FILE__, __LINE__, __FUNCTION__} + ::boost::leaf::exception
#define MU_LEAF_THROW_EXCEPTION ::boost::leaf::leaf_detail::throw_with_loc{__FILE__, __LINE__, __FUNCTION__} + ::boost::leaf::exception
#define MU_LEAF_LOG_ERROR(...)	((void)0)

namespace mu
{
	template<typename T>
	using optional_future = std::future<std::optional<T>>;

	template<typename T>
	inline auto future_is_ready(T const& f) noexcept -> bool
	{
		return f.wait_for(std::chrono::seconds(0)) == std::future_status::ready;
	}
} // namespace mu

namespace mu
{
	namespace details
	{
		template<typename T>
		class static_root_singleton
		{
		public:
			auto operator->() noexcept -> T*
			{
				return s_instance;
			}

			auto operator->() const noexcept -> const T*
			{
				return s_instance;
			}

			auto operator*() noexcept -> T&
			{
				return *s_instance;
			}

			auto operator*() const noexcept -> const T&
			{
				return *s_instance;
			}

			static_root_singleton() noexcept
			{
				static bool static_init = []() -> bool
				{
					s_instance = new (&s_instance_memory[0]) T();
					std::atexit(destroy);
					return true;
				}();
			}

		protected:
			static void destroy() noexcept
			{
				if (s_instance)
				{
					s_instance->~T();
					s_instance = nullptr;
				}
			}

		private:
			static inline uint64_t s_instance_memory[1 + (sizeof(T) / sizeof(uint64_t))];
			static inline T*	   s_instance = nullptr;
		};

		template<size_t POOL_COUNT>
		class singleton_cleanup_list
		{
		public:
			using cleanup_func = std::function<void()>;

			void push(cleanup_func f) noexcept
			{
				s_cleanup_funcs[s_cleanup_index.fetch_add(1)] = f;
			}

			~singleton_cleanup_list() noexcept
			{
				for (auto i = s_cleanup_index.load(); i > 0; --i)
				{
					s_cleanup_funcs[i - 1]();
				}
				s_cleanup_index.store(0);
			}

		private:
			std::array<cleanup_func, POOL_COUNT> s_cleanup_funcs;
			std::atomic_size_t					 s_cleanup_index = 0;
		};

		using singleton_cleanup_root = details::static_root_singleton<details::singleton_cleanup_list<1024>>;

		template<typename... T_DEPS>
		struct singleton_dependencies
		{
		private:
			template<typename T_ARG, typename... T_ARGS>
			static inline void update_dependencies_impl() noexcept
			{
				const typename T_ARG::type* t = T_ARG().get();
				if constexpr (sizeof...(T_ARGS) > 0)
				{
					update_dependencies_impl<T_ARGS...>();
				}
			}

		public:
			static inline void update() noexcept
			{
				if constexpr (sizeof...(T_DEPS) > 0)
				{
					update_dependencies_impl<T_DEPS...>();
				}
			}
		};

		template<typename T, typename... T_DEPS>
		struct singleton_factory
		{
			static inline auto create() noexcept -> T*
			{
				if constexpr (sizeof...(T_DEPS) > 0)
				{
					singleton_dependencies<T_DEPS...>::update();
				}
				return new T();
			}
		};

		template<typename T>
		struct virtual_singleton_factory
		{
			static auto create() noexcept -> T*;
		};
	} // namespace details

#define MU_DEFINE_VIRTUAL_SINGLETON(T, T_DERIVED)                                                                                \
	namespace mu                                                                                \
	{                                                                                \
		namespace details                                                                                \
		{                                                                                \
			template<>                                                                                \
			auto virtual_singleton_factory<T>::create() noexcept -> T*                                                                                \
			{                                                                                \
				return new T_DERIVED;                                                                                \
			}                                                                                \
		}                                                                                \
	}                                                                                \
	/**/

#define MU_DEFINE_VIRTUAL_SINGLETON_DEPS(T, T_DERIVED, ...)                                                                                \
	namespace mu                                                                                \
	{                                                                                \
		namespace details                                                                                \
		{                                                                                \
			template<>                                                                                \
			auto virtual_singleton_factory<T>::create() noexcept -> T*                                                                                \
			{                                                                                \
				singleton_dependencies<__VA_ARGS__>::update();                                                                                \
				return new T_DERIVED;                                                                                \
			}                                                                                \
		}                                                                                \
	}                                                                                \
	/**/

	namespace details
	{
		template<typename T, typename T_FACTORY>
		class singleton_base
		{
		public:
			using factory = T_FACTORY;
			using type	  = T;

			auto operator->() noexcept -> T*
			{
				return s_instance;
			}

			auto operator->() const noexcept -> const T*
			{
				return s_instance;
			}

			auto operator*() noexcept -> T&
			{
				return *s_instance;
			}

			auto operator*() const noexcept -> const T&
			{
				return *s_instance;
			}

			auto get() noexcept -> T*
			{
				return s_instance;
			}

			auto get() const noexcept -> const T*
			{
				return s_instance;
			}

			singleton_base() noexcept
			{
				static bool static_init = []() -> bool
				{
					s_instance = factory::create();
					singleton_cleanup_root()->push(
						[]() -> void
						{
							delete s_instance;
							s_instance = nullptr;
						});
					return true;
				}();
			}

		private:
			static inline T* s_instance;
		};
	} // namespace details

	template<typename T, typename... T_DEPENDENCIES>
	using singleton = details::singleton_base<T, details::singleton_factory<T, T_DEPENDENCIES...>>;

	template<typename T, typename... T_DEPENDENCIES>
	using virtual_singleton = details::singleton_base<T, details::virtual_singleton_factory<T>>;

	template<typename T_SINGLETON>
	class exported_singleton
	{
	public:
		using singleton_type = T_SINGLETON;
		using type			 = typename singleton_type::type;

		auto operator->() noexcept -> type*
		{
			return s_instance;
		}

		auto operator->() const noexcept -> const type*
		{
			return s_instance;
		}

		auto operator*() noexcept -> type&
		{
			return *s_instance;
		}

		auto operator*() const noexcept -> const type&
		{
			return *s_instance;
		}

		auto get() noexcept -> type*
		{
			return s_instance;
		}

		auto get() const noexcept -> const type*
		{
			return s_instance;
		}

		exported_singleton() noexcept
		{
			static bool static_init = []() -> bool
			{
				s_instance = get_instance();
				return true;
			}();
		}

	private:
		static inline type* s_instance;
		static type*		get_instance() noexcept;
	};

#define MU_EXPORT_SINGLETON(T)                                                                                \
	template<>                                                                                \
	auto T::get_instance() noexcept->T::type*                                                                                \
	{                                                                                \
		return singleton_type().get();                                                                                \
	}                                                                                \
	/**/

#define MU_EXPORT_SINGLETON_DEPS(T, ...)                                                                                \
	template<>                                                                                \
	auto T::get_instance() noexcept->T::type*                                                                                \
	{                                                                                \
		::mu::details::singleton_dependencies<__VA_ARGS__>::update();                                                                                \
		return singleton_type().get();                                                                                \
	}                                                                                \
	/**/

	namespace details
	{
		template<typename T>
		class static_root_thread_local_singleton
		{
		public:
			auto operator->() noexcept -> T*
			{
				return s_instance;
			}

			auto operator->() const noexcept -> const T*
			{
				return s_instance;
			}

			auto operator*() noexcept -> T&
			{
				return *s_instance;
			}

			auto operator*() const noexcept -> const T&
			{
				return *s_instance;
			}

			static_root_thread_local_singleton() noexcept
			{
				static bool static_init = []() -> bool
				{
					s_instance = new (&s_instance_memory[0]) T();
					std::atexit(destroy);
					return true;
				}();
			}

		protected:
			static void destroy() noexcept
			{
				if (s_instance)
				{
					s_instance->~T();
					s_instance = nullptr;
				}
			}

		private:
			static inline thread_local uint64_t s_instance_memory[1 + (sizeof(T) / sizeof(uint64_t))];
			static inline thread_local T*		s_instance = nullptr;
		};

		using thread_local_singleton_cleanup_root = static_root_thread_local_singleton<singleton_cleanup_list<1024>>;

		template<typename T, typename T_FACTORY>
		class thread_local_singleton_base
		{
		public:
			using type = T;
			auto operator->() noexcept -> T*
			{
				return s_instance;
			}

			auto operator->() const noexcept -> const T*
			{
				return s_instance;
			}

			auto operator*() noexcept -> T&
			{
				return *s_instance;
			}

			auto operator*() const noexcept -> const T&
			{
				return *s_instance;
			}

			auto get() noexcept -> T*
			{
				return s_instance;
			}

			auto get() const noexcept -> const T*
			{
				return s_instance;
			}

			thread_local_singleton_base() noexcept
			{
				static thread_local bool static_init = []() -> bool
				{
					s_instance = new T();
					thread_local_singleton_cleanup_root()->push(
						[]() -> void
						{
							delete s_instance;
							s_instance = nullptr;
						});
					return true;
				}();
			}

			static inline thread_local T* s_instance;
		};
	} // namespace details

	template<typename T, typename... T_DEPENDENCIES>
	using thread_local_singleton = details::thread_local_singleton_base<T, details::singleton_factory<T, T_DEPENDENCIES...>>;

	template<typename T, typename... T_DEPENDENCIES>
	using thread_local_virtual_singleton = details::thread_local_singleton_base<T, details::virtual_singleton_factory<T>>;

	template<typename T_SINGLETON>
	class exported_thread_local_singleton
	{
	public:
		using singleton_type = T_SINGLETON;
		using type			 = typename singleton_type::type;

		auto operator->() noexcept -> type*
		{
			return s_instance;
		}

		auto operator->() const noexcept -> const type*
		{
			return s_instance;
		}

		auto operator*() noexcept -> type&
		{
			return *s_instance;
		}

		auto operator*() const noexcept -> const type&
		{
			return *s_instance;
		}

		auto get() noexcept -> type*
		{
			return s_instance;
		}

		auto get() const noexcept -> const type*
		{
			return s_instance;
		}

		exported_thread_local_singleton() noexcept
		{
			static thread_local bool static_init = []() -> bool
			{
				s_instance = get_instance();
				return true;
			}();
		}

	private:
		static inline thread_local type* s_instance;
		static type*					 get_instance() noexcept;
	};

#define MU_EXPORT_THREAD_LOCAL_SINGLETON(T)                                                                                \
	template<>                                                                                \
	auto ::mu::exported_thread_local_singleton<T>::get_instance() noexcept->::mu::exported_thread_local_singleton<T>::type*                                                        \
	{                                                                                \
		return singleton_type().get();                                                                                \
	}                                                                                \
	/**/
} // namespace mu

namespace mu
{
	namespace time
	{
		auto performance_frequency() noexcept -> int64_t;
		auto get_now() noexcept -> int64_t;
		void sleep(const int64_t milliseconds) noexcept;
		void micro_sleep(const int64_t tx) noexcept;
		void init() noexcept;
		void calibrate() noexcept;
		void set_high_resolution_timer() noexcept;
		auto release_high_resolution_timer() noexcept -> leaf::result<void>;

		class moment
		{
			friend class run_time;
			friend class long_clock;

		protected:
			int64_t value;

			moment(const int64_t v) noexcept : value(v) { }

			moment(const long double v) noexcept : value(static_cast<int64_t>(v)) { }

			friend auto now() noexcept -> moment;

			template<typename T>
			friend auto seconds(const T&) noexcept -> moment;

			template<typename T>
			friend auto milliseconds(const T&) noexcept -> moment;

			template<typename T>
			friend auto microseconds(const T&) noexcept -> moment;

			template<typename T>
			friend auto nanoseconds(const T&) noexcept -> moment;

		public:
			moment() noexcept : value(0) { }

			inline auto operator=(const moment& rhs) noexcept -> moment&
			{
				value = rhs.value;
				return *this;
			}

			inline auto operator+=(const moment& rhs) noexcept -> moment&
			{
				value += rhs.value;
				return *this;
			}

			inline auto operator-=(const moment& rhs) noexcept -> moment&
			{
				value += rhs.value;
				return *this;
			}

			template<typename T>
			inline auto operator*=(const T& rhs) noexcept -> moment&
			{
				value = static_cast<int64_t>(static_cast<long double>(value) * static_cast<long double>(rhs));
				return *this;
			}

			template<typename T>
			inline auto operator/=(const T& rhs) noexcept -> moment&
			{
				value = static_cast<int64_t>(static_cast<long double>(value) / static_cast<long double>(rhs));
				return *this;
			}

			inline auto add(const moment& rhs) const noexcept -> moment
			{
				return moment{value + rhs.value};
			}

			inline auto sub(const moment& rhs) const noexcept -> moment
			{
				return moment{value - rhs.value};
			}

			template<typename T>
			inline auto div(const T& rhs) const noexcept -> moment
			{
				return moment{static_cast<long double>(value) / static_cast<long double>(rhs)};
			}

			template<typename T>
			inline auto mul(const T& rhs) const noexcept -> moment
			{
				return moment{static_cast<long double>(value) * static_cast<long double>(rhs)};
			}

			template<typename T>
			inline auto as_seconds() const noexcept -> T
			{
				return static_cast<T>(static_cast<long double>(value) / static_cast<long double>(performance_frequency()));
			}

			template<typename T>
			inline auto as_milliseconds() const noexcept -> T
			{
				return static_cast<T>(static_cast<long double>(value) / ((static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull))));
			}

			template<typename T>
			inline auto as_microseconds() const noexcept -> T
			{
				return static_cast<T>(static_cast<long double>(value) / ((static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull * 1000ull))));
			}

			template<typename T>
			inline auto as_nanoseconds() const noexcept -> T
			{
				return static_cast<T>(
					static_cast<long double>(value) / ((static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull * 1000ull * 1000ull))));
			}

			template<typename T>
			inline auto as_ticks() const noexcept -> T
			{
				return static_cast<T>(value);
			}

			template<typename T>
			inline auto set_seconds(const T& seconds) noexcept -> moment&
			{
				value = static_cast<int64_t>((static_cast<long double>(performance_frequency())) * static_cast<long double>(seconds));
				return *this;
			}

			template<typename T>
			inline auto set_milliseconds(const T& seconds) noexcept -> moment&
			{
				value = static_cast<int64_t>((static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull)) * static_cast<long double>(seconds));
				return *this;
			}

			template<typename T>
			inline auto set_microseconds(const T& seconds) noexcept -> moment&
			{
				value = static_cast<int64_t>((static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull * 1000ull)) * static_cast<long double>(seconds));
				return *this;
			}

			template<typename T>
			inline auto set_nanoseconds(const T& seconds) noexcept -> moment&
			{
				value = static_cast<int64_t>(
					(static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull * 1000ull * 1000ull)) * static_cast<long double>(seconds));
				return *this;
			}

			template<typename T>
			inline auto set_ticks(const T& tx) noexcept -> moment&
			{
				value = static_cast<int64_t>(tx);
				return *this;
			}

			inline auto now() noexcept -> moment&
			{
				value = get_now();
				return *this;
			}
		};

		class long_clock
		{
			uint64_t value;
			moment	 last_moment;

		public:
			long_clock() noexcept : value(0) { }

			template<typename T>
			inline auto as_seconds() const noexcept -> T
			{
				return static_cast<T>(static_cast<long double>(value) / static_cast<long double>(1000ull));
			}

			template<typename T>
			inline auto as_milliseconds() const noexcept -> T
			{
				return static_cast<T>(static_cast<long double>(value));
			}

			template<typename T>
			inline auto as_microseconds() const noexcept -> T
			{
				return static_cast<T>(static_cast<long double>(value) * static_cast<long double>(1000ull));
			}

			void update() noexcept;
		};

		inline auto operator<(const moment& lhs, const moment& rhs) noexcept -> bool
		{
			return lhs.as_ticks<int64_t>() < rhs.as_ticks<int64_t>();
		}

		inline auto operator<=(const moment& lhs, const moment& rhs) noexcept -> bool
		{
			return lhs.as_ticks<int64_t>() <= rhs.as_ticks<int64_t>();
		}

		inline auto operator>(const moment& lhs, const moment& rhs) noexcept -> bool
		{
			return lhs.as_ticks<int64_t>() > rhs.as_ticks<int64_t>();
		}

		inline auto operator>=(const moment& lhs, const moment& rhs) noexcept -> bool
		{
			return lhs.as_ticks<int64_t>() >= rhs.as_ticks<int64_t>();
		}

		inline auto operator==(const moment& lhs, const moment& rhs) noexcept -> bool
		{
			return lhs.as_ticks<int64_t>() == rhs.as_ticks<int64_t>();
		}

		inline auto operator+(const moment& lhs, const moment& rhs) noexcept -> moment
		{
			return lhs.add(rhs);
		}

		inline auto operator-(const moment& lhs, const moment& rhs) noexcept -> moment
		{
			return lhs.sub(rhs);
		}

		template<typename T>
		inline auto operator/(const moment& lhs, const T& rhs) noexcept -> moment
		{
			return lhs.div(rhs);
		}

		template<typename T>
		inline auto operator*(const moment& lhs, const T& rhs) noexcept -> moment
		{
			return lhs.mul(rhs);
		}

		inline auto now() noexcept -> moment
		{
			moment t;
			return t.now();
		}

		template<typename T>
		inline auto seconds(const T& val) noexcept -> moment
		{
			moment t;
			return t.set_seconds(val);
		}

		template<typename T>
		inline auto milliseconds(const T& val) noexcept -> moment
		{
			moment t;
			return t.set_milliseconds(val);
		}

		template<typename T>
		inline auto microseconds(const T& val) noexcept -> moment
		{
			moment t;
			return t.set_microseconds(val);
		}

		template<typename T>
		inline auto nanoseconds(const T& val) noexcept -> moment
		{
			moment t;
			return t.set_nanoseconds(val);
		}

		template<typename T>
		inline auto ticks(const T& val) noexcept -> moment
		{
			moment t;
			return t.set_ticks(val);
		}

		template<typename T>
		inline void micro_sleep_seconds(const T& seconds) noexcept
		{
			micro_sleep(static_cast<int64_t>((static_cast<long double>(performance_frequency())) * static_cast<long double>(seconds)));
		}

		template<typename T>
		inline void micro_sleep_milliseconds(const T& seconds) noexcept
		{
			micro_sleep(static_cast<int64_t>((static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull)) * static_cast<long double>(seconds)));
		}

		template<typename T>
		inline void micro_sleep_microseconds(const T& seconds) noexcept
		{
			micro_sleep(
				static_cast<int64_t>((static_cast<long double>(performance_frequency()) / static_cast<long double>(1000ull * 1000ull)) * static_cast<long double>(seconds)));
		}

		template<typename T>
		inline void micro_sleep_ticks(const T& tx) noexcept
		{
			micro_sleep(static_cast<int64_t>(tx));
		}

		inline void micro_sleep(const moment& t) noexcept
		{
			micro_sleep(t.as_ticks<int64_t>());
		}

		template<typename T>
		inline void sleepseconds(const T& seconds) noexcept
		{
			sleep(moment().set_seconds(seconds).template as_milliseconds<int64_t>());
		}

		template<typename T>
		inline void sleepmilliseconds(const T& seconds) noexcept
		{
			sleep(moment().set_milliseconds(seconds).template as_milliseconds<int64_t>());
		}

		template<typename T>
		inline void sleepmicroseconds(const T& seconds) noexcept
		{
			sleep(moment().set_microseconds(seconds).template as_milliseconds<int64_t>());
		}

		template<typename T>
		inline void sleepticks(const T& tx) noexcept
		{
			sleep(moment().set_ticks(tx).template as_milliseconds<int64_t>());
		}

		inline void sleep(const moment& t) noexcept
		{
			sleep(t.template as_milliseconds<int64_t>());
		}

		inline void long_clock::update() noexcept
		{
			const moment n = now();
			const moment d = n - last_moment;
			last_moment	   = n;
			value += d.template as_milliseconds<int64_t>();
		}
	} // namespace time
} // namespace mu

namespace mu
{
	enum class messagebox_result : int
	{
		ok = 0,
		cancel,
		yes,
		no,
		quit,
		none,
		error
	};

	enum class messagebox_style : int
	{
		info = 0,
		warning,
		error,
		question
	};

	enum class messagebox_buttons : int
	{
		ok = 0,
		okcancel,
		yesno,
		quit
	};

	namespace details
	{
		template<typename T>
		struct future_helper
		{
			bool m_state = false;
			T	 m_future;

			using value_type = decltype(m_future.get());

			auto is_active() const noexcept -> bool
			{
				return m_state;
			}

			auto is_ready() const noexcept -> bool
			{
				if (m_state)
				{
					return future_is_ready(m_future);
				}
				return false;
			}

			auto acquire_value() noexcept -> value_type
			{
				m_state = false;
				return m_future.get();
			}

			auto get_value() noexcept -> std::optional<value_type>
			{
				if (is_ready())
				{
					return acquire_value();
				}
				return std::nullopt;
			}
		};
	} // namespace details

	namespace details
	{
		auto show_messagebox(const char* message, const char* title, messagebox_style style, messagebox_buttons buttons) noexcept -> std::future<messagebox_result>;
	}

	using messagebox_future = details::future_helper<decltype(details::show_messagebox("", "", messagebox_style(), messagebox_buttons()))>;
	inline auto show_messagebox(const char* message, const char* title, messagebox_style style, messagebox_buttons buttons) noexcept -> messagebox_future
	{
		return {true, details::show_messagebox(message, title, style, buttons)};
	}

	namespace details
	{
		auto show_file_open_dialog(std::string_view origin, std::string_view filter) noexcept -> optional_future<std::string>;
		auto show_file_open_multiple_dialog(std::string_view origin, std::string_view filter) noexcept -> optional_future<std::vector<std::string>>;
		auto show_file_save_dialog(std::string_view origin, std::string_view filter) noexcept -> optional_future<std::string>;
		auto show_path_dialog(std::string_view origin, std::string_view filter) noexcept -> optional_future<std::string>;
	} // namespace details

	using file_open_dialog_future = details::future_helper<decltype(details::show_file_open_dialog(std::string_view(), std::string_view()))>;
	inline auto show_file_open_dialog(std::string_view origin, std::string_view filter) noexcept -> file_open_dialog_future
	{
		return {true, details::show_file_open_dialog(origin, filter)};
	}

	using file_open_multiple_dialog_future = details::future_helper<decltype(details::show_file_open_multiple_dialog(std::string_view(), std::string_view()))>;
	inline auto show_file_open_multiple_dialog(std::string_view origin, std::string_view filter) noexcept -> file_open_multiple_dialog_future
	{
		return {true, details::show_file_open_multiple_dialog(origin, filter)};
	}

	using file_save_dialog_future = details::future_helper<decltype(details::show_file_save_dialog(std::string_view(), std::string_view()))>;
	inline auto show_file_save_dialog(std::string_view origin, std::string_view filter) noexcept -> file_save_dialog_future
	{
		return {true, details::show_file_save_dialog(origin, filter)};
	}

	using path_dialog_future = details::future_helper<decltype(details::show_path_dialog(std::string_view(), std::string_view()))>;
	inline auto show_path_dialog(std::string_view origin, std::string_view filter) noexcept -> path_dialog_future
	{
		return {true, details::show_path_dialog(origin, filter)};
	}
} // namespace mu

namespace backward
{
	class StackTrace;
}

namespace mu
{
	namespace debug
	{
		namespace details
		{
			struct logger_interface
			{
				logger_interface()			= default;
				virtual ~logger_interface() = default;

				virtual auto stdout_logger() noexcept -> std::shared_ptr<spdlog::logger> = 0;
				virtual auto stderr_logger() noexcept -> std::shared_ptr<spdlog::logger> = 0;
			};
		} // namespace details

		using logger = mu::exported_singleton<mu::virtual_singleton<details::logger_interface>>;

		void log_stack_trace(spdlog::logger& l, spdlog::level::level_enum lvl, unsigned int level_skip) noexcept;

	} // namespace debug

	static inline auto error_handlers = std::make_tuple(
		[](runtime_error::not_specified x, leaf::e_source_location sl)
		{
			debug::logger()->stderr_logger()->log(spdlog::level::err, "{0} :: {1} -> {2} : runtime_error :: not_specified", sl.line, sl.file, sl.function);
		},
		[](common_error x, leaf::e_source_location sl)
		{
			debug::logger()->stderr_logger()->log(spdlog::level::err, "{0} :: {1} -> {2} : common_error", sl.line, sl.file, sl.function);
		},
		[]
		{
			debug::logger()->stderr_logger()->log(spdlog::level::err, "???");
		});

	template<class TryBlock>
	constexpr inline auto try_handle(TryBlock&& try_block) -> typename std::decay<decltype(std::declval<TryBlock>()().value())>::type
	{
		return leaf::try_handle_all(try_block, error_handlers);
	}

} // namespace mu

namespace mu
{
	template<typename T_CONTAINER, typename T_FUNC, typename T_FUNC_RET>
	auto for_some_return(T_CONTAINER& container, T_FUNC func, T_FUNC_RET&& true_return_value, T_FUNC_RET&& false_return_value) noexcept -> T_FUNC_RET
	{
		for (auto& itor : container)
		{
			if (func(itor))
			{
				return true_return_value;
			}
		}

		return false_return_value;
	}

	template<typename T_CONTAINER, typename T_FUNC, typename T_FUNC_RET>
	auto for_some_optional_return(T_CONTAINER& container, T_FUNC func, T_FUNC_RET&& false_return_value) noexcept -> T_FUNC_RET
	{
		for (auto& itor : container)
		{
			if (auto res = func(itor))
			{
				return *res;
			}
		}

		return false_return_value;
	}
} // namespace mu

namespace mu
{
	void enable_dpi_awareness() noexcept;
	auto get_dpi_scale_for_monitor(void* monitor) noexcept -> float;
	auto get_dpi_scale_for_hwnd(void* hwnd) noexcept -> float;
} // namespace mu

namespace mu
{

} // namespace mu


text:
#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <cblas.h>
#include <string.h> // memcpy
#include <math.h>

#include "cgraph.h"
#include "cg_operation.h"
#include "cg_types.h"
#include "cg_variables.h"
#include "cg_errors.h"
#include "cg_constants.h"
#include "cg_enums.h"
#include "cg_factory.h"
#include "cg_math.h"

#if cg_float == float
#define cblas_dcopy cblas_scopy
#define cblas_dscal cblas_sscal
#define cblas_dgemm cblas_sgemm
#define cblas_ddot cblas_sdot
#define cblas_dgemv cblas_sgemv
#endif

int selectContext() {
    return 0;
}

void copyDataToHost(CGResultNode* res){}

/*
 * *********************
 * Multiplication
 * *********************
 */

/*
 * d.d
 */
CGResultNode* mulDD(CGDouble* M, CGDouble* V, CGraph* graph, CGNode* parentNode){
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = M->value * V->value;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * d.V == V.d
 */
CGResultNode* mulDV(CGDouble* a, CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = V->len;
    Y->data = y;

    if(a->value != 0){
        cblas_dcopy(V->len, V->data, 1, y, 1);
        cblas_dscal(V->len, a->value, y, 1);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * d.M == M.d
 */
CGResultNode* mulDM(CGDouble* a, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->cols*M->rows;
    cg_float value = a->value;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = y;

    if(a->value != 0){
        cblas_dcopy(size, M->data, 1, y, 1);
        cblas_dscal(size, a->value, y, 1);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M*v
 * broadcasting
 */
CGResultNode* mulMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->cols != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Add M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M->cols*M->rows;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M->data[i] * V->data[(i)%V->len];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * Mv
 */
/*
CGResultNode* mulMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
	if(M->cols != V->len){
		char msg[MAX_ERR_FMT_LEN];
		snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Add M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
		return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
	}

	uint64_t size = M->cols*M->rows;
	cg_float* res = calloc(size, sizeof(cg_float));

	CGMatrix* Y = calloc(1, sizeof(CGMatrix));
	Y->rows = M->rows;
	Y->cols = M->cols;
	Y->data = res;

	uint64_t i = 0;

	for(;i<size;i++){
		res[i] = M->data[i] * V->data[(i)%V->len];
	}

	CGResultNode* result = calloc(1, sizeof(CGResultNode));
	result->type = CGVT_MATRIX;
	result->value = Y;

	return result;
}
*/

/*
 * M*N
 * broadcasting
 */
CGResultNode* mulMM(CGMatrix* M1, CGMatrix* M2, CGraph* graph, CGNode* parentNode){
    if((M1->rows != M2->rows) && (M1->cols != M2->cols)){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot multiply M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M1->rows, M1->cols, M2->rows, M2->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M1->cols*M1->rows;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M1->rows;
    Y->cols = M1->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M1->data[i] * M2->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M.v
 */
CGResultNode* dotMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->cols != 1){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate M(%"PRIu64", %"PRIu64") DOT V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    cg_float* y = calloc(M->rows*V->len, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows= M->rows;
    Y->cols  = V->len;
    Y->data = y;

    /*
    cblas_dgemv(CblasRowMajor,
            CblasNoTrans, M->rows, M->cols, 1.0, M->data,
            M->cols, V->data, 1, 0.0, Y->data, 1);
    */

    cblas_dgemm(CblasRowMajor,
                CblasNoTrans, CblasNoTrans, M->rows, V->len, M->cols,
                1.0, M->data, M->cols, V->data, V->len, 0, y, Y->cols);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * v.M
 */
CGResultNode* dotVM(CGVector* V, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate V(%"PRIu64") DOT M(%"PRIu64", %"PRIu64")", V->len, M->rows, M->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    cg_float* y = calloc(M->cols, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = M->cols;
    Y->data = y;

    int m = 1;
    int n = M->cols;
    int k = V->len;

    cblas_dgemm(CblasRowMajor,
                CblasNoTrans, CblasNoTrans, m, n, k, 1, V->data, k, M->data, n, 0, Y->data, n);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M.N
 */
CGResultNode* dotMM(CGMatrix* M, CGMatrix* N, CGraph* graph, CGNode* parentNode){
    if(M->cols != N->rows){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot DOT M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M->rows, M->cols, N->rows, N->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    cg_float* z = calloc(M->rows*N->cols, sizeof(cg_float));

    CGMatrix* Z = calloc(1, sizeof(CGMatrix));
    Z->rows = M->rows;
    Z->cols = N->cols;

    cblas_dgemm(CblasRowMajor,
                CblasNoTrans, CblasNoTrans, M->rows, N->cols, M->cols,
                1.0, M->data, M->cols, N->data, N->cols, 0, z, Z->cols);

    Z->data = z;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Z;

    return result;
}

/*
 * V.V
 */
CGResultNode* crossVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot compute Cross Product of V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = V1->len;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = V1->data[i] * V2->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Dot Product
 * *********************
 */

/*
 * V.V
 */
CGResultNode* dotVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot compute Dot Product of V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;
    cg_float res = 0;

    CGDouble* Y = calloc(1, sizeof(CGDouble));

    Y->value = cblas_ddot(V1->len, V1->data, 1, V2->data, 1);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Division
 * *********************
 */

/*
 *  d/d
 */
CGResultNode* divDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    if(D2->value == 0.0){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Diving by zero D/(D=0)");
        return returnResultError(graph, CGET_DIVIDE_BY_ZERO, parentNode, msg);
    }

    cg_float res = D1->value / D2->value;

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 *  V/d
 */
CGResultNode* divVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    if(D->value == 0.0){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Diving by zero V/(D=0)");
        return returnResultError(graph, CGET_DIVIDE_BY_ZERO, parentNode, msg);
    }

    cg_float* res = calloc(V->len, sizeof(cg_float));
    cg_float value = D->value;

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = res;
    Y->len = V->len;

    cblas_dcopy(V->len, V->data, 1, res, 1);
    cblas_dscal(V->len, 1.0/D->value, res, 1);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 *  V/V
 */
CGResultNode* divVV(CGVector* V, CGVector* D, CGraph* graph, CGNode* parentNode){
    if(V->len != D->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate  V(%"PRIu64") DIV V(%"PRIu64")", V->len, D->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    cg_float* res = calloc(V->len, sizeof(cg_float));

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = res;
    Y->len = V->len;

    uint64_t i = 0;
    for(;i<V->len;i++){
        res[i] = V->data[i]/D->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 *  d/V element-wise
 */
CGResultNode* divDV(CGDouble* D, CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* res = calloc(V->len, sizeof(cg_float));
    cg_float value = D->value;

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = res;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        res[i] = value / V->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 *  M/d
 */
CGResultNode* divMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    if(D->value == 0.0){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Diving by zero M/(D=0)");
        return returnResultError(graph, CGET_DIVIDE_BY_ZERO, parentNode, msg);
    }

    uint64_t size = M->rows*M->cols;
    cg_float* res = calloc(size, sizeof(cg_float));
    cg_float value = D->value;

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    cblas_dcopy(size, M->data, 1, res, 1);
    cblas_dscal(size, 1.0/D->value, res, 1);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M/v
 */
CGResultNode* divMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot calculate  M(%"PRIu64", %"PRIu64") DIV V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M->rows*M->cols;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(; i < size;i++)
    {
        res[i] = M->data[i] / V->data[i/M->cols];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 *  d/M element-wise
 */
CGResultNode* divDM(CGDouble* D, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;
    cg_float* res = calloc(size, sizeof(cg_float));
    cg_float value = D->value;

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = value / M->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Addition
 * *********************
 */

/*
 * d+d
 */
CGResultNode* addDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    cg_float res = D1->value + D2->value;

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * V+d
 */
CGResultNode* addVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float* res = calloc(V->len, sizeof(cg_float));
    cg_float value = D->value;

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = res;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        res[i] = V->data[i] + value;
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M+d
 */
CGResultNode* addMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;
    cg_float* res = calloc(size, sizeof(cg_float));
    cg_float value = D->value;

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M->data[i] + value;
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * V+V
 */
CGResultNode* addVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Add V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = V1->len;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = V1->data[i] + V2->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M+V
 */
CGResultNode* addMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->cols != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Add M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M->cols*M->rows;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M->data[i] + V->data[(i)%V->len];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M+M
 */
CGResultNode* addMM(CGMatrix* M1, CGMatrix* M2, CGraph* graph, CGNode* parentNode){
    if((M1->rows != M2->rows) && (M1->cols != M2->cols)){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot add M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M1->rows, M1->cols, M2->rows, M2->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M1->cols*M1->rows;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M1->rows;
    Y->cols = M1->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M1->data[i] + M2->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Substraction
 * *********************
 */

/*
 * d-d
 */
CGResultNode* subDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    cg_float res = D1->value - D2->value;

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * V-d
 */
CGResultNode* subVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float* res = calloc(V->len, sizeof(cg_float));
    cg_float value = D->value;

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = res;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        res[i] = V->data[i] - value;
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * d-V
 */
CGResultNode* subDV(CGDouble* D, CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* res = calloc(V->len, sizeof(cg_float));
    cg_float value = D->value;

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = res;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        res[i] = value - V->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M-d
 */
CGResultNode* subMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;
    cg_float* res = calloc(size, sizeof(cg_float));
    cg_float value = D->value;

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M->data[i] - value;
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * D-M
 */
CGResultNode* subDM(CGDouble* D, CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;
    cg_float* res = calloc(size, sizeof(cg_float));
    cg_float value = D->value;

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = value - M->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * V-V
 */
CGResultNode* subVV(CGVector* V1, CGVector* V2, CGraph* graph, CGNode* parentNode){
    if(V1->len != V2->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract V(%"PRIu64") by V(%"PRIu64")", V1->len, V2->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = V1->len;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = V1->len;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = V1->data[i] - V2->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M-M
 */
CGResultNode* subMM(CGMatrix* M1, CGMatrix* M2, CGraph* graph, CGNode* parentNode){
    if((M1->rows != M2->rows) && (M1->cols != M2->cols)){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract M(%"PRIu64", %"PRIu64") by N(%"PRIu64", %"PRIu64")", M1->rows, M1->cols, M2->rows, M2->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M1->cols*M1->rows;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M1->rows;
    Y->cols = M1->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M1->data[i] - M2->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M-V
 */
CGResultNode* subMV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract M(%"PRIu64", %"PRIu64") by V(%"PRIu64")", M->rows, M->cols, V->len);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M->cols*M->rows;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = M->data[i] - V->data[(i)/M->cols];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * V-M
 */
CGResultNode* subVM(CGVector* V, CGMatrix* M, CGraph* graph, CGNode* parentNode){

    if(M->rows != V->len){
        char msg[MAX_ERR_FMT_LEN];
        snprintf(msg, MAX_ERR_FMT_LEN, "Cannot Substract V(%"PRIu64") by M(%"PRIu64", %"PRIu64")", V->len, M->rows, M->cols);
        return returnResultError(graph, CGET_INCOMPATIBLE_DIMENTIONS_EXCEPTION, parentNode, msg);
    }

    uint64_t size = M->cols*M->rows;
    cg_float* res = calloc(size, sizeof(cg_float));

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = V->data[(i)/M->cols] - M->data[i];
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Power
 * *********************
 */

/*
 * d^d
 */
CGResultNode* powDD(CGDouble* D1, CGDouble* D2, CGraph* graph, CGNode* parentNode){
    cg_float res = pow(D1->value, D2->value);

    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = res;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * V^d
 */
CGResultNode* powVD(CGVector* V, CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float* res = calloc(V->len, sizeof(cg_float));
    cg_float value = D->value;

    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = res;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        res[i] = pow(V->data[i], value);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * M^d
 */
CGResultNode* powMD(CGMatrix* M, CGDouble* D, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;
    cg_float* res = calloc(size, sizeof(cg_float));
    cg_float value = D->value;

    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->rows = M->rows;
    Y->cols = M->cols;
    Y->data = res;

    uint64_t i = 0;

    for(;i<size;i++){
        res[i] = pow(M->data[i], value);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Transpose Multiplication
 * *********************
 */
/*
 * M^T.v
 */
CGResultNode* mulMtV(CGMatrix* M, CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->len = V->len;
    Y->data = y;

    cblas_dgemv(CblasRowMajor,
                CblasTrans, M->rows, M->cols, 1.0, M->data,
                V->len, V->data, 1, 0.0, Y->data, 1);

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * exp(d)
 */
CGResultNode* expD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = exp(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * exp(v)
 */
CGResultNode* expV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = y;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        y[i] = exp(V->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * exp(M)
 */
CGResultNode* expM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = M->rows;
    Y->cols = M->cols;

    uint64_t i = 0;

    for(;i<size;i++){
        y[i] = exp(M->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * log(d)
 */
CGResultNode* logD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = log(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * log(V)
 */
CGResultNode* logV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = y;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        y[i] = log(V->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * exp(M)
 */
CGResultNode* logM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = M->rows;
    Y->cols = M->cols;

    uint64_t i = 0;

    for(;i<size;i++){
        y[i] = log(M->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Sin
 * *********************
 */

/*
 * sin(d)
 */
CGResultNode* sinD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = sin(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * sin(V)
 */
CGResultNode* sinV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = y;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        y[i] = sin(V->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * sin(M)
 */
CGResultNode* sinM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = M->rows;
    Y->cols = M->cols;

    uint64_t i = 0;

    for(;i<size;i++){
        y[i] = sin(M->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * cos
 * *********************
 */

/*
 * cos(d)
 */
CGResultNode* cosD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = cos(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * cos(V)
 */
CGResultNode* cosV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = y;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        y[i] = cos(V->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * cos(M)
 */
CGResultNode* cosM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = M->rows;
    Y->cols = M->cols;

    uint64_t i = 0;

    for(;i<size;i++){
        y[i] = cos(M->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * tan
 * *********************
 */

/*
 * tan(d)
 */
CGResultNode* tanD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = tan(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * tan(V)
 */
CGResultNode* tanV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = y;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        y[i] = tan(V->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * tan(M)
 */
CGResultNode* tanM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = M->rows;
    Y->cols = M->cols;

    uint64_t i = 0;

    for(;i<size;i++){
        y[i] = tan(M->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * tanh
 * *********************
 */

/*
 * tanh(d)
 */
CGResultNode* tanhD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = tanh(D->value);
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * tanh(V)
 */
CGResultNode* tanhV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    CGVector* Y = calloc(1, sizeof(CGVector));
    Y->data = y;
    Y->len = V->len;

    uint64_t i = 0;

    for(;i<V->len;i++){
        y[i] = tanh(V->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = Y;

    return result;
}

/*
 * tanh(M)
 */
CGResultNode* tanhM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = M->rows;
    Y->cols = M->cols;

    uint64_t i = 0;

    for(;i<size;i++){
        y[i] = tanh(M->data[i]);
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * Invert
 * *********************
 */

// TODO

/*
 * *********************
 * Transpose
 * *********************
 */

/*
 * D^t
 */
CGResultNode* transposeD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    CGDouble* R = calloc(1, sizeof(CGDouble));
    R->value = D->value;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = R;

    return result;
}

/*
 * V^t
 */
CGResultNode* transposeV(CGVector* V, CGraph* graph, CGNode* parentNode){
    uint64_t size = V->len;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = V->len;
    Y->cols = 1;
    Y->data = y;

    uint64_t i = 0;
    uint64_t j = 0;

    memcpy(Y->data, V->data, V->len*sizeof(cg_float));

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * M^t
 */
CGResultNode* transposeM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t size = M->rows*M->cols;

    cg_float* y = calloc(size, sizeof(cg_float));
    CGMatrix* Y = calloc(1, sizeof(CGMatrix));
    Y->data = y;
    Y->rows = M->cols;
    Y->cols = M->rows;
    Y->data = y;

    uint64_t i = 0;
    uint64_t j = 0;

    for(;i<M->rows;i++){
        for(j = 0;j<M->cols;j++){
            y[j*Y->cols+i] = M->data[i*M->cols +j];
        }
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_MATRIX;
    result->value = Y;

    return result;
}

/*
 * *********************
 * SUM
 * *********************
 */

/*
 * sum(D)
 */
CGResultNode* sumD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float y = D->value;
    CGDouble* Y = calloc(1, sizeof(CGDouble));
    Y->value = y;

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * sum(V)
 */
CGResultNode* sumV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float y = 0;
    CGDouble* Y = calloc(1, sizeof(CGDouble));

    uint64_t i = 0;

    for(;i<V->len;i++){
        y += V->data[i];
    }

    Y->value = y;
    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_DOUBLE;
    result->value = Y;

    return result;
}

/*
 * TODO:
 * sum(M, axis=axis)
 */

CGResultNode* sumM(CGMatrix* M, CGraph* graph, CGNode* parentNode, uint8_t axis){
    CGVector* V = calloc(1, sizeof(CGVector));

    if(axis == 0){
        uint64_t len = M->cols;
        cg_float* y = calloc(len, sizeof(cg_float));

        V->data = y;
        V->len = len;

        uint64_t i = 0;

        for(;i<M->rows*M->cols;i++){
            y[i%len] += M->data[i];
        }
    }
    else {
        uint64_t len = M->rows;
        cg_float* y = calloc(len, sizeof(cg_float));

        V->data = y;
        V->len = len;

        uint64_t i = 0;

        for(;i<M->rows*M->cols;i++){
            y[i/M->cols] += M->data[i];
        }
    }

    CGResultNode* result = calloc(1, sizeof(CGResultNode));
    result->type = CGVT_VECTOR;
    result->value = V;

    return result;
}

CGResultNode* max(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            cg_float m = v->data[0];

            for (;i < v->len; i++){
                if(v->data[i] > m)
                    m = v->data[i];
            }

            Y->value = m;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = M->data[i];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i%len] < M->data[i])
                        y[i%len] = M->data[i];
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = M->data[i*M->cols];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i/M->cols] < M->data[i])
                        y[i/M->cols] = M->data[i];
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* min(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            cg_float m = v->data[0];

            for (;i < v->len; i++){
                if(v->data[i] > m)
                    m = v->data[i];
            }

            Y->value = m;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = M->data[i];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i%len] > M->data[i])
                        y[i%len] = M->data[i];
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = M->data[i*M->cols];
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(y[i/M->cols] > M->data[i])
                        y[i/M->cols] = M->data[i];
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* argmax(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            uint64_t maxIdx = 0;

            for (;i < v->len; i++){
                if(v->data[i] > v->data[maxIdx])
                    maxIdx = i;
            }

            Y->value = maxIdx;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = i;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i%len]] < M->data[i])
                        y[i%len] = i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->cols;i++){
                    y[i] = ((int)y[i])/M->cols;
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = i*M->cols;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i/M->cols]] < M->data[i])
                        y[i/M->cols] = i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->rows;i++){
                    y[i] = ((int)y[i])%M->cols;
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* argmin(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = 0;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 1;
            uint64_t minIdx = 0;

            for (;i < v->len; i++){
                if(v->data[i] > v->data[minIdx])
                    minIdx = (cg_float)i;
            }

            Y->value = minIdx;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->cols;i++){
                    y[i] = i;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i%len]] > M->data[i])
                        y[i%len] = (cg_float)i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->cols;i++){
                    y[i] = ((int)y[i])/M->cols;
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows;i++){
                    y[i] = i*M->cols;
                }

                for(i=0;i<M->rows*M->cols;i++){
                    if(M->data[(uint64_t)y[i/M->cols]] > M->data[i])
                        y[i/M->cols] = (cg_float)i;
                }

                // transform vector representation into matrix
                for(i=0;i<M->rows;i++){
                    y[i] = ((int)y[i])%M->cols;
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

CGResultNode* mean(CGNode* X, CGraph* graph){
    CGResultNode* res = computeCGNode(graph, X->axop->uhs);

    switch(res->type){
        case CGVT_DOUBLE:
        {
            cg_float y = 0;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            Y->value = ((CGDouble*)res->value)->value;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_VECTOR:
        {
            CGVector* v = (CGVector*)res->value;
            CGDouble* Y = calloc(1, sizeof(CGDouble));

            uint64_t i = 0;
            cg_float m = v->data[0];

            for (;i < v->len; i++){
                m += v->data[i];
            }

            m /= v->len;

            Y->value = m;
            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_DOUBLE;
            result->value = Y;

            return result;
        }

        case CGVT_MATRIX: {
            CGMatrix* M = (CGMatrix*)res->value;
            CGVector* V = calloc(1, sizeof(CGVector));

            if(X->axop->axis == 0){
                uint64_t len = M->cols;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows*M->cols;i++){
                    y[i%len] += M->data[i];
                }

                for(i=0;i<M->cols;i++){
                    y[i] /= M->rows;
                }
            }
            else {
                uint64_t len = M->rows;
                cg_float* y = calloc(len, sizeof(cg_float));

                V->data = y;
                V->len = len;

                uint64_t i = 0;

                for(;i<M->rows*M->cols;i++){
                    y[i/M->cols] += M->data[i];
                }

                for(i=0;i<M->rows;i++){
                    y[i] /= M->cols;
                }
            }

            CGResultNode* result = calloc(1, sizeof(CGResultNode));
            result->type = CGVT_VECTOR;
            result->value = V;

            return result;

        }
    }
}

/*
 * relu(d)
 */
CGResultNode* reluD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float val = D->value;
    return makeDoubleResultNode(__cg_relu(val));
}

/*
 * relu(V)
 */
CGResultNode* reluV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    __cg_map_array(V->data, y, V->len, __cg_relu);
    return makeVectorResultNode(V->len, y);
}

/*
 * relu(M)
 */
CGResultNode* reluM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t len = M->cols*M->rows;
    cg_float* y = calloc(len, sizeof(cg_float));
    __cg_map_array(M->data, y, len, __cg_relu);
    return makeMatrixResultNode(M->rows, M->cols, y);
}

/*
 * softplus(d)
 */
CGResultNode* softplusD(CGDouble* D, CGraph* graph, CGNode* parentNode){
    cg_float val = D->value;
    return makeDoubleResultNode(__cg_softplus(val));
}

/*
 * softplus(V)
 */
CGResultNode* softplusV(CGVector* V, CGraph* graph, CGNode* parentNode){
    cg_float* y = calloc(V->len, sizeof(cg_float));
    __cg_map_array(V->data, y, V->len, __cg_softplus);
    return makeVectorResultNode(V->len, y);
}

/*
 * softplus(M)
 */
CGResultNode* softplusM(CGMatrix* M, CGraph* graph, CGNode* parentNode){
    uint64_t len = M->cols*M->rows;
    cg_float* y = calloc(len, sizeof(cg_float));
    __cg_map_array(M->data, y, len, __cg_softplus);
    return makeMatrixResultNode(M->rows, M->cols, y);
}

#if cg_float == float
#undef cblas_dcopy
#undef cblas_dscal
#undef cblas_dgemm
#undef cblas_ddot
#undef cblas_dgemv
#endif

ext: c

max_stars_repo_name: Enehcruon/cgraph

max_issues_repo_name: Enehcruon/cgraph

max_forks_repo_name: praisethemoon/ccgraph


text:
/* multilargenlin/convergence.c
 *
 * Copyright (C) 2015 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_multilarge_nlinear.h>

static double scaled_infnorm(const gsl_vector *x, const gsl_vector *g);

/*
gsl_multilarge_nlinear_test()
  Convergence tests for nonlinear least squares minimization

(1) |dx_i| <= xtol * (1 + |x_i|) for all i
(2) || g .* x ||_inf <= gtol ||f||^2
(3) ||f(x+dx) - f(x)|| <= ftol * max(||f(x)||, 1)

Inputs: xtol - tolerance for step size
        gtol - tolerance for gradient vector
        ftol - tolerance for residual vector
        info - (output)
               1 - stopped by small x step
               2 - stopped by small gradient
               3 - stopped by small residual vector change
        w    - workspace
*/

int
gsl_multilarge_nlinear_test (const double xtol, const double gtol,
                             const double ftol, int *info,
                             const gsl_multilarge_nlinear_workspace * w)
{
  int status;
  double gnorm, fnorm, phi;

  *info = 0;

  status = gsl_multifit_test_delta(w->dx, w->x, xtol*xtol, xtol);
  if (status == GSL_SUCCESS)
    {
      *info = 1;
      return GSL_SUCCESS;
    }

  /* compute gnorm = max_i( g_i * max(x_i, 1) ) */
  gnorm = scaled_infnorm(w->x, w->g);

  /* compute fnorm = ||f|| */
  fnorm = gsl_blas_dnrm2(w->f);
  phi = 0.5 * fnorm * fnorm;

#if 0
  fprintf(stderr, "gnorm = %.12e fnorm = %.12e gnorm/phi = %.12e\n", gnorm, fnorm, gnorm / phi);
#endif

  if (gnorm <= gtol * GSL_MAX(phi, 1.0))
    {
      *info = 2;
      return GSL_SUCCESS;
    }

#if 0
  if (dfnorm <= ftol * GSL_MAX(fnorm, 1.0))
    {
      *info = 3;
      return GSL_SUCCESS;
    }
#endif

  return GSL_CONTINUE;
}

static double
scaled_infnorm(const gsl_vector *x, const gsl_vector *g)
{
  const size_t n = x->size;
  size_t i;
  double norm = 0.0;

  for (i = 0; i < n; ++i)
    {
      double xi = GSL_MAX(gsl_vector_get(x, i), 1.0);
      double gi = gsl_vector_get(g, i);
      double tmp = fabs(xi * gi);

      if (tmp > norm)
        norm = tmp;
    }

  return norm;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multilarge.h>
#include <gsl/gsl_blas.h>

/* function to be fitted */
double
func(const double t)
{
  double x = sin(10.0 * t);
  return exp(x*x*x);
}

/* construct a row of the least squares matrix */
int
build_row(const double t, gsl_vector *row)
{
  const size_t p = row->size;
  double Xj = 1.0;
  size_t j;

  for (j = 0; j < p; ++j)
    {
      gsl_vector_set(row, j, Xj);
      Xj *= t;
    }

  return 0;
}

int
solve_system(const int print_data, const gsl_multilarge_linear_type * T,
             const double lambda, const size_t n, const size_t p,
             gsl_vector * c)
{
  const size_t nblock = 5;         /* number of blocks to accumulate */
  const size_t nrows = n / nblock; /* number of rows per block */
  gsl_multilarge_linear_workspace * w =
    gsl_multilarge_linear_alloc(T, p);
  gsl_matrix *X = gsl_matrix_alloc(nrows, p);
  gsl_vector *y = gsl_vector_alloc(nrows);
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  const size_t nlcurve = 200;
  gsl_vector *reg_param = gsl_vector_alloc(nlcurve);
  gsl_vector *rho = gsl_vector_alloc(nlcurve);
  gsl_vector *eta = gsl_vector_alloc(nlcurve);
  size_t rowidx = 0;
  double rnorm, snorm, rcond;
  double t = 0.0;
  double dt = 1.0 / (n - 1.0);

  while (rowidx < n)
    {
      size_t nleft = n - rowidx;         /* number of rows left to accumulate */
      size_t nr = GSL_MIN(nrows, nleft); /* number of rows in this block */
      gsl_matrix_view Xv = gsl_matrix_submatrix(X, 0, 0, nr, p);
      gsl_vector_view yv = gsl_vector_subvector(y, 0, nr);
      size_t i;

      /* build (X,y) block with 'nr' rows */
      for (i = 0; i < nr; ++i)
        {
          gsl_vector_view row = gsl_matrix_row(&Xv.matrix, i);
          double fi = func(t);
          double ei = gsl_ran_gaussian (r, 0.1 * fi); /* noise */
          double yi = fi + ei;

          /* construct this row of LS matrix */
          build_row(t, &row.vector);

          /* set right hand side value with added noise */
          gsl_vector_set(&yv.vector, i, yi);

          if (print_data && (i % 100 == 0))
            printf("%f %f\n", t, yi);

          t += dt;
        }

      /* accumulate (X,y) block into LS system */
      gsl_multilarge_linear_accumulate(&Xv.matrix, &yv.vector, w);

      rowidx += nr;
    }

  if (print_data)
    printf("\n\n");

  /* compute L-curve */
  gsl_multilarge_linear_lcurve(reg_param, rho, eta, w);

  /* solve large LS system and store solution in c */
  gsl_multilarge_linear_solve(lambda, c, &rnorm, &snorm, w);

  /* compute reciprocal condition number */
  gsl_multilarge_linear_rcond(&rcond, w);

  fprintf(stderr, "=== Method %s ===\n", gsl_multilarge_linear_name(w));
  fprintf(stderr, "condition number = %e\n", 1.0 / rcond);
  fprintf(stderr, "residual norm    = %e\n", rnorm);
  fprintf(stderr, "solution norm    = %e\n", snorm);

  /* output L-curve */
  {
    size_t i;
    for (i = 0; i < nlcurve; ++i)
      {
        printf("%.12e %.12e %.12e\n",
               gsl_vector_get(reg_param, i),
               gsl_vector_get(rho, i),
               gsl_vector_get(eta, i));
      }
    printf("\n\n");
  }

  gsl_matrix_free(X);
  gsl_vector_free(y);
  gsl_multilarge_linear_free(w);
  gsl_rng_free(r);
  gsl_vector_free(reg_param);
  gsl_vector_free(rho);
  gsl_vector_free(eta);

  return 0;
}

int
main(int argc, char *argv[])
{
  const size_t n = 50000;   /* number of observations */
  const size_t p = 16;      /* polynomial order + 1 */
  double lambda = 0.0;      /* regularization parameter */
  gsl_vector *c_tsqr = gsl_vector_alloc(p);
  gsl_vector *c_normal = gsl_vector_alloc(p);

  if (argc > 1)
    lambda = atof(argv[1]);

  /* solve system with TSQR method */
  solve_system(1, gsl_multilarge_linear_tsqr, lambda, n, p, c_tsqr);

  /* solve system with Normal equations method */
  solve_system(0, gsl_multilarge_linear_normal, lambda, n, p, c_normal);

  /* output solutions */
  {
    gsl_vector *v = gsl_vector_alloc(p);
    double t;

    for (t = 0.0; t <= 1.0; t += 0.01)
      {
        double f_exact = func(t);
        double f_tsqr, f_normal;

        build_row(t, v);
        gsl_blas_ddot(v, c_tsqr, &f_tsqr);
        gsl_blas_ddot(v, c_normal, &f_normal);

        printf("%f %e %e %e\n", t, f_exact, f_tsqr, f_normal);
      }

    gsl_vector_free(v);
  }

  gsl_vector_free(c_tsqr);
  gsl_vector_free(c_normal);

  return 0;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim