gctl/lib/core/sparray.h

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_SPARRAY_H
#define _GCTL_SPARRAY_H

// library's head files
#include "array.h"
#include "enum.h"

namespace gctl
{
	template <typename NodeValueType>
	struct arr_node
	{
		int id;
		NodeValueType val;
		arr_node(){};
		arr_node(int in_id, NodeValueType in_val)
		{
			id = in_id;
			val = in_val;
		}
	};

	template <typename ArrayValueType>
	class sparray
	{
	protected:
		std::vector<arr_node<ArrayValueType> > arr1d; ///< 稀疏数组存储向量
		int full_length; ///< 稀疏数组的全长
		ArrayValueType zero_val; ///< 稀疏数组的零值
	public:
		sparray(); ///< 默认构造函数
		sparray(int in_len, ArrayValueType null_val); ///< 构造函数：预留数组空间并制定零值
		sparray(const array<ArrayValueType> &b, ArrayValueType null_val, 
			ArrayValueType off_set); ///< 构造函数：从数组初始化
		sparray(const std::vector<ArrayValueType> &b, ArrayValueType null_val, 
			ArrayValueType off_set); ///< 构造函数：从向量初始化
		sparray(const sparray<ArrayValueType> &b); ///< 拷贝构造函数
		sparray<ArrayValueType> &operator= (const sparray<ArrayValueType> &b); ///< 重载赋值操作符
		~sparray(); ///< 析构函数

		void malloc(int in_len, ArrayValueType null_val); ///< 预留数组空间并制定零值
		void malloc(const sparray<ArrayValueType> &b); ///< 从稀疏数组初始化
		void malloc(const array<ArrayValueType> &b, ArrayValueType null_val, 
			ArrayValueType off_set); ///< 从数组初始化
		void malloc(const std::vector<ArrayValueType> &b, ArrayValueType null_val, 
			ArrayValueType off_set); ///< 从向量初始化
		void clear(); ///< 清空数组内容

		ArrayValueType value(unsigned int in_id) const; ///< 返回索引值位置的元素值
		arr_node<ArrayValueType> at(unsigned int in_id) const; ///< 返回非零元素在索引位置的值

		void set(int in_id, ArrayValueType in_val); ///< 在队中添加元素
		void remove(int in_id); ///< 删除索引值位置的元素

		bool empty() const; ///< 数组是否为空
		int size() const; ///< 返回非零元素的个数
		int capcity() const; ///< 返回数组的预定长度
		ArrayValueType zero_value() const; ///< 返回0值

		void show_list(int st_id = 0, int ed_id = 0) const; ///< 显示存储的元素
		void show_full() const; ///< 显示完整的数组
		void export_dense(array<ArrayValueType> &b, double multipler = 1.0, value_operator_e type = ReplaceVal) const; ///< 输出为全长度的数组

		void copy_to(sparray<ArrayValueType> &b, double multipler = 1.0) const; ///< 拷贝到新的稀疏数组
	private:
		typename std::vector<arr_node<ArrayValueType> >::iterator find_iterator(int in_id);
		typename std::vector<arr_node<ArrayValueType> >::const_iterator find_const_iterator(int in_id) const; ///< 寻找索引值位置的元素并返回其迭代器值
		typename std::vector<arr_node<ArrayValueType> >::iterator find_gap(int in_id); ///< 寻找包含索引值的左右两个元素并返回右侧元素的迭代器
	};

	template <typename ArrayValueType>
	sparray<ArrayValueType>::sparray()
	{
		full_length = -1;
	}

	template <typename ArrayValueType>
	sparray<ArrayValueType>::sparray(int in_len, ArrayValueType null_val)
	{
		full_length = -1;
		malloc(in_len, null_val);
	}

	template <typename ArrayValueType>
	sparray<ArrayValueType>::sparray(const array<ArrayValueType> &b, ArrayValueType null_val, 
		ArrayValueType off_set)
	{
		full_length = -1;
		malloc(b, null_val, off_set);
	}

	template <typename ArrayValueType>
	sparray<ArrayValueType>::sparray(const std::vector<ArrayValueType> &b, ArrayValueType null_val, 
		ArrayValueType off_set)
	{
		full_length = -1;
		malloc(b, null_val, off_set);
	}

	template <typename ArrayValueType>
	sparray<ArrayValueType>::sparray(const sparray<ArrayValueType> &b)
	{
		full_length = -1;
		malloc(b);
	}

	template <typename ArrayValueType>
	sparray<ArrayValueType> &sparray<ArrayValueType>::operator= (const sparray<ArrayValueType> &b)
	{
		if (!empty()) clear();

		zero_val = b.zero_value();
		full_length = b.capcity();
		arr1d.resize(b.size());
		for (int i = 0; i < b.size(); i++)
		{
			arr1d[i] = b.at(i);
		}
		return this;
	}

	template <typename ArrayValueType>
	sparray<ArrayValueType>::~sparray()
	{
		clear();
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::malloc(int in_len, ArrayValueType null_val)
	{
		if (in_len <= 0)
		{
			std::string err_str = "Invalid array size. From void sparray<T>::malloc(...)";
			throw runtime_error(err_str);
		}

		full_length = in_len;
		zero_val = null_val;
		arr1d.reserve(in_len);
		return;
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::malloc(const sparray<ArrayValueType> &b)
	{
		if (b.empty())
		{
			std::string err_str = "The input array is empty. From void sparray<T>::malloc(...)";
			throw runtime_error(err_str);
		}
		else
		{
			zero_val = b.zero_value();
			full_length = b.capcity();
			arr1d.resize(b.size());
			for (int i = 0; i < b.size(); i++)
			{
				arr1d[i] = b.at(i);
			}
		}
		return;
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::malloc(const array<ArrayValueType> &b, ArrayValueType null_val, 
		ArrayValueType off_set)
	{
		if (b.empty())
		{
			std::string err_str = "The input array is empty. From void sparray<T>::malloc(...)";
			throw runtime_error(err_str);
		}
		else
		{
			full_length = b.size();
			zero_val = null_val;
			arr1d.reserve(full_length);
			for (int i = 0; i < b.size(); i++)
			{
				if (b[i] >= zero_val + off_set || b[i] <= zero_val - off_set)
				{
					set(i, b[i]);
				}
			}
		}
		return;
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::malloc(const std::vector<ArrayValueType> &b, ArrayValueType null_val, 
		ArrayValueType off_set)
	{
		if (b.empty())
		{
			std::string err_str = "The input array is empty. From void sparray<T>::malloc(...)";
			throw runtime_error(err_str);
		}

		full_length = b.size();
		zero_val = null_val;
		arr1d.reserve(full_length);
		for (int i = 0; i < b.size(); i++)
		{
			if (b.at(i) >= zero_val + off_set || 
				b.at(i) <= zero_val - off_set)
				insert(i, b.at(i));
		}
		return;
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::clear()
	{
		//full_length = 0; // 注意清理稀疏数组时不清理此参数，否则需要重新声明此数组
		arr1d.clear();
		std::vector<arr_node<ArrayValueType> >().swap(arr1d);
		return;
	}

	template <typename ArrayValueType>
	ArrayValueType sparray<ArrayValueType>::value(unsigned int in_id) const
	{
		typename std::vector<arr_node<ArrayValueType> >::const_iterator iter = find_const_iterator(in_id);
		if (iter != arr1d.end())
			return iter->val;
		return zero_val;
	}

	template <typename ArrayValueType>
	arr_node<ArrayValueType> sparray<ArrayValueType>::at(unsigned int in_id) const
	{
		return arr1d[in_id];
	}

	/**
	 * @brief      在队中添加元素 替换已有元素的值
	 *
	 * @param[in]  in_id   输入索引值
	 * @param[in]  in_val  输入值
	 */
	template <typename ArrayValueType>
	void sparray<ArrayValueType>::set(int in_id, ArrayValueType in_val)
	{
		if (in_id < 0 || in_id >= full_length)
		{
			std::string err_str = "The index is out of range. From void sparray<T>::set(...)";
			throw runtime_error(err_str);
		}

		if (in_val == zero_val)
			return;

		if (arr1d.empty())
		{
			arr1d.push_back(arr_node<ArrayValueType>(in_id, in_val));
		}
		else if (in_id > arr1d.back().id) // 输入索引大于向量内的最大索引值
		{
			arr1d.push_back(arr_node<ArrayValueType>(in_id, in_val));
		}
		else if (in_id == arr1d.back().id)
		{
			arr1d.back().val = in_val;
		}
		else if (in_id < arr1d.front().id)
		{
			arr1d.insert(arr1d.begin(), arr_node<ArrayValueType>(in_id, in_val));
		}
		else if (in_id == arr1d.front().id)
		{
			arr1d.front().val = in_val;
		}
		// 在中间插入元素比较耗时
		else
		{
			typename std::vector<arr_node<ArrayValueType> >::iterator iter = find_iterator(in_id);
			if (iter != arr1d.end()) iter->val = in_val;
			else
			{
				iter = find_gap(in_id);
				arr1d.insert(iter, arr_node<ArrayValueType>(in_id, in_val));
			}
		}
		return;
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::remove(int in_id)
	{
		typename std::vector<arr_node<ArrayValueType> >::iterator iter = find_iterator(in_id);
		if (iter != arr1d.end())
			arr1d.erase(iter);
		return;
	}

	template <typename ArrayValueType>
	bool sparray<ArrayValueType>::empty() const
	{
		if (arr1d.empty()) return true;
		else return false;
	}

	template <typename ArrayValueType>
	int sparray<ArrayValueType>::size() const
	{
		if (arr1d.empty()) return 0;
		else return arr1d.size();
	}

	/**
	 * @brief     返回数组的预定长度
	 *
	 * @return     数组全长
	 */
	template <typename ArrayValueType>
	int sparray<ArrayValueType>::capcity() const
	{
		return full_length;
	}

	template <typename ArrayValueType>
	ArrayValueType sparray<ArrayValueType>::zero_value() const
	{
		return zero_val;
	}

	/**
	 * @brief      显示存储的元素
	 *
	 */
	template <typename ArrayValueType>
	void sparray<ArrayValueType>::show_list(int st_id, int ed_id) const
	{
		if (ed_id == 0)
			ed_id = arr1d.size();

		for (int i = st_id; i < ed_id; i++)
		{
			std::cout << "( " << arr1d[i].id << ", " << arr1d[i].val << ")" << std::endl;
		}
		return;
	}

	/**
	 * @brief      显示完整的数组
	 */
	template <typename ArrayValueType>
	void sparray<ArrayValueType>::show_full() const
	{
		int st = 0;
		for (int i = 0; i < arr1d.size(); i++)
		{
			for (int s = st; s < arr1d[i].id; s++)
				std::cout << zero_val << " ";
			std::cout << arr1d[i].val << " ";
			st = arr1d[i].id+1;
		}

		if (st <= full_length-1)
		{
			for (int s = st; s < full_length; s++)
				std::cout << zero_val << " ";
		}
		std::cout << std::endl;
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::export_dense(array<ArrayValueType> &b, double multipler, 
		value_operator_e type) const
	{
		std::string err_str;
		if (type == AppendVal)
		{
			if (b.size() != full_length)
			{
				err_str = "Invalid output array size. From void sparray<T>::export_dense(...)";
				throw runtime_error(err_str);
			}

			for (int i = 0; i < arr1d.size(); i++)
			{
				b.at(arr1d[i].id) += arr1d[i].val * multipler;
			}
		}
		else if (type == ReplaceVal)
		{
			b.resize(full_length, zero_val);

			for (int i = 0; i < arr1d.size(); i++)
			{
				b.at(arr1d[i].id) = arr1d[i].val * multipler;
			}
		}
		else
		{
			err_str = "Invalid value type for the output array. From void sparray<T>::export_dense(...)";
			throw runtime_error(err_str);
		}
		return;
	}

	template <typename ArrayValueType>
	void sparray<ArrayValueType>::copy_to(sparray<ArrayValueType> &b, double multipler) const
	{
		if (!b.empty()) b.clear();
		b.malloc(full_length, zero_val);

		for (int i = 0; i < arr1d.size(); i++)
		{
			b.insert(arr1d[i].id, multipler*arr1d[i].val);
		}
		return;
	}

	/******************************/
	/*      以下是类的私有函数       */
	/******************************/

	template <typename ArrayValueType>
	typename std::vector<arr_node<ArrayValueType> >::iterator sparray<ArrayValueType>::find_iterator(int in_id)
	{
		if (in_id == arr1d.front().id)
			return arr1d.begin();
		else if (in_id < arr1d.front().id)
			return arr1d.end();
		else if (in_id == arr1d.back().id)
			return --arr1d.end();
		else if (in_id > arr1d.back().id)
			return arr1d.end();
		else
		{
			typename std::vector<arr_node<ArrayValueType> >::iterator l_iter = arr1d.begin(), c_iter;

			int l_id = 0, u_id = arr1d.size()-1;

			int jump_size;
			while (u_id-l_id > 1)
			{
				jump_size = (u_id-l_id)/2;
				c_iter = l_iter + jump_size;
				if (c_iter->id == in_id)
				{
					return c_iter;
				}
				else if (c_iter->id > in_id)
				{
					u_id -= jump_size;
				}
				else if (c_iter->id < in_id)
				{
					l_id += jump_size;
					l_iter = c_iter;
				}
			}
			return arr1d.end();
		}
	}

	template <typename ArrayValueType>
	typename std::vector<arr_node<ArrayValueType> >::const_iterator 
		sparray<ArrayValueType>::find_const_iterator(int in_id) const
	{
		if (in_id == arr1d.front().id)
			return arr1d.begin();
		else if (in_id < arr1d.front().id)
			return arr1d.end();
		else if (in_id == arr1d.back().id)
			return --arr1d.end();
		else if (in_id > arr1d.back().id)
			return arr1d.end();
		else
		{
			typename std::vector<arr_node<ArrayValueType> >::const_iterator l_iter = arr1d.begin(), c_iter;

			int l_id = 0, u_id = arr1d.size()-1;

			int jump_size;
			while (u_id-l_id > 1)
			{
				jump_size = (u_id-l_id)/2;
				c_iter = l_iter + jump_size;
				if (c_iter->id == in_id)
				{
					return c_iter;
				}
				else if (c_iter->id > in_id)
				{
					u_id -= jump_size;
				}
				else if (c_iter->id < in_id)
				{
					l_id += jump_size;
					l_iter = c_iter;
				}
			}
			return arr1d.end();
		}
	}

	template <typename ArrayValueType>
	typename std::vector<arr_node<ArrayValueType> >::iterator sparray<ArrayValueType>::find_gap(int in_id)
	{
		if (in_id < arr1d.front().id)
			return arr1d.end();
		else if (in_id > arr1d.back().id)
			return arr1d.end();

		typename std::vector<arr_node<ArrayValueType> >::iterator l_iter = arr1d.begin(), c_iter;
		int l_id = 0, u_id = arr1d.size()-1;

		int jump_size;
		while (u_id-l_id > 0)
		{
			jump_size = (u_id-l_id)/2;
			c_iter = l_iter + jump_size;
			if (c_iter->id < in_id && in_id < (c_iter+1)->id)
			{
				return c_iter+1;
			}
			else if (c_iter->id > in_id)
			{
				u_id -= jump_size;
			}
			else if ((c_iter+1)->id < in_id)
			{
				l_id += jump_size;
				l_iter = c_iter;
			}
		}
		return arr1d.end();
	}

	/******************以下是基于sparray的sparray2d类型********************/
	/**
	 * @brief      基于向量的二维稀疏矩阵。
	 * 
	 * 与十字链表实现的二维稀疏矩阵不同，sparray2d主要基于向量实现。每行单独为一个稀疏
	 * 的向量数组，行与行之间无连接，因此无法做稀疏矩阵之间的乘法与向量前乘矩阵等操作。
	 * 但是从每行尾部插入元素的效率与每行从左至右访问元素的效率都很高。适合对列操作无需求，且（或）
	 * 矩阵元素经常变化或者被访问的情况下使用。
	 * 
	 * 在每行中间替换元素值的复杂度为O(logN)，插入新值的复杂度为O(2*logN)，
	 * 但由于向量需要顺序偏移插入点后的元素，实际效率会略慢。使用中应尽量从行尾插入元素。
	 * 
	 *
	 * @tparam     ArrayValueType  矩阵元素值类型
	 */
	template <typename ArrayValueType>
	class sparray2d
	{
	private:
		std::vector<sparray<ArrayValueType> > arr2d;
		ArrayValueType zero_val;
		int full_rlength, full_clength;
	public:
		sparray2d();
		sparray2d(int rnum, int cnum, ArrayValueType null_val);
		~sparray2d();

		void malloc(int rnum, int cnum, ArrayValueType null_val);

		void clear();

		sparray<ArrayValueType> *at(unsigned int index);

		int size();
	};

	template <typename ArrayValueType>
	sparray2d<ArrayValueType>::sparray2d()
	{
		full_rlength = full_clength = 0;
	}

	template <typename ArrayValueType>
	sparray2d<ArrayValueType>::sparray2d(int rnum, int cnum, ArrayValueType null_val)
	{
		malloc(rnum, cnum, null_val);
	}

	template <typename ArrayValueType>
	sparray2d<ArrayValueType>::~sparray2d()
	{
		clear();
	}

	template <typename ArrayValueType>
	void sparray2d<ArrayValueType>::malloc(int rnum, int cnum, ArrayValueType null_val)
	{
		if (rnum <= 0 || cnum <= 0)
		{
			std::string err_str = "Invalid matrix size. From void spmatrix<T>::malloc(...)";
			throw runtime_error(err_str);
		}

		full_rlength = rnum; full_clength = cnum;
		zero_val = null_val;
		arr2d.resize(rnum);
		for (int i = 0; i < rnum; i++)
		{
			arr2d[i].malloc(cnum, null_val);
		}
		return;
	}

	template <typename ArrayValueType>
	void sparray2d<ArrayValueType>::clear()
	{
		full_rlength = full_clength = 0;
		for (int i = 0; i < arr2d.size(); i++)
		{
			arr2d[i].clear();
		}
		arr2d.clear();
		std::vector<sparray<ArrayValueType> >().swap(arr2d);
		return;
	}

	template <typename ArrayValueType>
	sparray<ArrayValueType> *sparray2d<ArrayValueType>::at(unsigned int index)
	{
		return &arr2d[index];
	}

	template <typename ArrayValueType>
	int sparray2d<ArrayValueType>::size()
	{
		return arr2d.size();
	}
}

#endif // _GCTL_SPARRAY_H