gctl/lib/core/array.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_ARRAY_H
#define _GCTL_ARRAY_H

// library's head files
#include "../gctl_config.h"
#include "enum.h"
#include "macro.h"
#include "vector_t.h"
#include "exceptions.h"
#include "random"

#ifdef GCTL_EIGEN
/*The followings bypass a bug that VsCode's IntelliSense reports incorrect errors when using eigen3 library*/
#if __INTELLISENSE__
#undef __ARM_NEON
#undef __ARM_NEON__
#endif
#include "Eigen/Dense"
#include "Eigen/Sparse"
#endif // GCTL_EIGEN

#ifdef GCTL_OPENMP
#include "omp.h"
#endif // GCTL_OPENMP

namespace gctl
{
	template <typename ArrValType> class array;

	/*
	 * Here we define some commonly used array types.
	 */
	typedef array<int> _1i_array; ///< 1-D integer array
	typedef array<float> _1f_array; ///< 1-D single-precision floating-point array
	typedef array<double> _1d_array; ///< 1-D double-precision floating-point array
	typedef array<std::string> _1s_array; ///< 1-D string array
	typedef array<std::complex<float>> _1cf_array; ///< 1-D complex float array
	typedef array<std::complex<double>> _1cd_array; ///< 1-D complex double array

	/**
	 * @brief Iterator template. This could be used to enable the use of iteration algorithms for array-like objects.
	 * 
	 */
	template <typename T>
	class iterator
	{
	private:
		T *iter;

	public:
		iterator(T *para, size_t n) {iter = para + n;}
		virtual ~iterator(){}

		T &operator*() {return *iter;}
		bool operator!=(const iterator &that) {return this->iter != that.iter;}
		iterator &operator++() {++iter; return *this;}		
	};

	/**
	 * @brief      1-D array class template
	 *
	 * @tparam     ArrValType     template array
	 */
	template <typename ArrValType>
	class array
	{
	protected:
		ArrValType *val_; ///< Pointer of the array
		size_t length_; ///< Length of the array

	public:
		/**
		 * @brief      Default constructor.
		 */
		array();

		/**
		 * @brief      Construct an array with the given length.
		 *
		 * @param[in]  len   Length of the array. Must be equal to or bigger than zero.
		 */
		array(size_t len);

		/**
		 * @brief      Construct an array with the given length and initial values.
		 *
		 * @param[in]  len       Length of the array. Must be equal to or bigger than zero.
		 * @param[in]  init_val  Initial value of the elements.
		 */
		array(size_t len, ArrValType init_val);

		/**
		 * @brief      Construct an array with the given length and fill elements using a native c++ array.
		 * 
		 * @note       Note that if you use an native c++ array to initialize the array template, the pointer 
		 * must be first variable and the length is last variable. This design is mainly to avoid unclear function 
		 * calls. Consider an array of pointers array<double*> A(10, nullptr) means that the initial value of all 
		 * elements is nullptr and array<double*> A((double**)B, 10) means that the element value of A is initialized 
		 * with array B If the above expression is written as array<double*> A(10, B), the compiler cannot call the 
		 * function correctly. Because it cannot be judged whether we want to use the initial value or 
		 * the array to initialize the array.
		 *
		 * @param      init_array  Pointer of the native c++ array.
		 * @param[in]  len         Length of the input c++ array. Must be equal to or bigger than zero.
		 */
		array(ArrValType *init_array, size_t len);

		/**
		 * @brief      Copy constructor.
		 *
		 * @param[in]  b     Original array.
		 */
		array(const array<ArrValType> &b);

		/**
		 * @brief Construct a new array object from std::initializer
		 * 
		 * @param init_val Initial values
		 */
		array(std::initializer_list<ArrValType> init_val);

		/**
		 * @brief Construct a new array object as a sequence
		 * 
		 * @param s Size of the sequence.
		 * @param st Start value
		 * @param inc Invrease interval
		 */
		array(size_t s, ArrValType st, ArrValType inc);

		/**
		 * @brief      Overloaded assignment operator for the array template.
		 *
		 * @param[in]  b     Original array.
		 *
		 * @return     Target array.
		 */
		array<ArrValType>& operator= (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded assignment operator for the array template.
		 * 
		 * @param v          Initial value.
		 * @return     Target array.
		 */
		array<ArrValType>& operator= (ArrValType v);

		/**
		 * @brief      Overloaded summation operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType> operator+ (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded minus operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType> operator- (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded multiple operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType> operator* (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded division operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType> operator/ (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded summation operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType>& operator+= (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded summation operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType>& operator-= (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded multiple operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType>& operator*= (const array<ArrValType> &b);

		/**
		 * @brief      Overloaded division operator for the array template.
		 * 
		 * @param a    Input array
		 * @param b    Inpur array
		 * @return Target array
		 */
		array<ArrValType>& operator/= (const array<ArrValType> &b);

		/**
		 * @brief      Destructor
		 */
		virtual ~array();

		/**
		 * @brief      Set the length of the array to the given value.
		 * 
		 * @note       If the input length is not equal to the existing length of the array, 
		 * the previous array is destroyed and reconstructed.
		 *
		 * @param[in]  len   Length of the array. Must be equal to or bigger than zero.
		 */
		void resize(size_t len);

		/**
		 * @brief      Set the length of the array to the given value, and fill elements with initial values.
		 * 
		 * @note       If the input length is not equal to the existing length of the array, 
		 * the previous array is destroyed and reconstructed.
		 *
		 * @param[in]  len       Length of the array. Must be equal to or bigger than zero.
		 * @param[in]  init_val  Initial value of the elements.
		 */
		void resize(size_t len, ArrValType init_val);

		/**
		 * @brief      Set the length of the array to the given value, and fill elements using a native c++ array.
		 * 
		 * @note       If the input length is not equal to the existing length of the array, 
		 * the previous array is destroyed and reconstructed.
		 *
		 * @param      init_array  Pointer of the native c++ array.
		 * @param[in]  len         Length of the input c++ array. Must be equal to or bigger than zero.
		 */
		void resize(ArrValType *init_array, size_t len);

		/**
		 * @brief      Copy a array to the array.
		 * 
		 * @note       If the length of the input array is not equal to the existing length of the array, 
		 * the previous array is destroyed and reconstructed.
		 *
		 * @param[in]  b     Original array
		 */
		void resize(const array<ArrValType> &b);

		/**
		 * @brief Resize a new array object from std::initializer
		 * 
		 * @param init_val Initial values
		 */
		void resize(std::initializer_list<ArrValType> init_val);

		/**
		 * @brief Resize a new array object as a sequence
		 * 
		 * @param s Size of the sequence.
		 * @param st Start value
		 * @param inc Invrease interval
		 */
		void resize(size_t s, ArrValType st, ArrValType inc);

		/**
		 * @brief      Clear memory space and reset variables.
		 */
		void clear();

		/**
		 * @brief      Set all elements to the given value.
		 *
		 * @param[in]  in_val    Input value.
		 */
		void assign_all(ArrValType in_val);

		/**
		 * @brief Set value to segments of an array.
		 * 
		 * @param in_val Input value
		 * @param st Starting index
		 * @param range segment range
		 */
		void assign(ArrValType in_val, size_t st, size_t range);
		
		/**
		 * @brief Append an array to the end of the existing one
		 * 
		 * @param b Input array
		 */
		void append_array(const array<ArrValType> &b);

		/**
		 * @brief Extract an array
		 * 
		 * @param b   Output array. Must be initialized before use.
		 * @param st  Start index
		 */
		void extract_array(array<ArrValType> &b, size_t st = 0) const;

		/**
		 * @brief Extract an array
		 * 
		 * @param b   Output array
		 * @param ids Index of the extracting elements. Must be initialized before use.
		 */
		void extract_array(array<ArrValType> &b, const array<size_t> &ids) const;

		/**
		 * @brief Insert data from an input array
		 * 
		 * @param b The input array
		 * @param st The inserting point. The default is zero
		 */
		void insert_array(array<ArrValType> &b, size_t st = 0);

		/**
		 * @brief      Get element value at the given index.
		 * 
		 * @note       This function could be called by a pointer.
		 *
		 * @param[in]  index  index value.
		 *
		 * @return     element's value.
		 */
		ArrValType &at(size_t index);

		/**
		 * @brief      Get element value at the given index (Constant).
		 * 
		 * @note       This function could be called by a pointer.
		 *
		 * @param[in]  index  index value.
		 *
		 * @return     element's value.
		 */
		ArrValType &at(size_t index) const;

		/**
		 * @brief      Get element value at the given index.
		 * 
		 * @note       This function could not be called by a pointer.
		 *
		 * @param[in]  index  index value.
		 *
		 * @return     element's value.
		 */
		ArrValType &operator[](size_t index);

		/**
		 * @brief      Get element value at the given index  .
		 * 
		 * @note       This function could not be called by a pointer.
		 *
		 * @param[in]  index  index value.
		 *
		 * @return     element's value.
		 */
		ArrValType &operator[](size_t index) const;
		
		/**
		 * @brief      Get the first element
		 * 
		 * @return     element value
		 */
		ArrValType front() const;

		/**
		 * @brief      Get the last element
		 *
		 * @return     element value
		 */
		ArrValType back() const;

		/**
		 * @brief Return the begining of the iterator
		 * 
		 * @return iterator
		 */
		iterator<ArrValType> begin();

		/**
		 * @brief Return the ending of the iterator
		 * 
		 * @return iterator
		 */
		iterator<ArrValType> end();

		/**
		 * @brief      Get the pointer to a element at the given index.
		 *
		 * @param[in]  index  index value.
		 *
		 * @return     Pointer to the element
		 */
		ArrValType *get(size_t index = 0) const;

		/**
		 * @brief      Discriminate if the array is empty.
		 *
		 * @return     True for empty. False other wise.
		 */
		bool empty() const;

		/**
		 * @brief      Return the length of the class member array.
		 *
		 * @return     Length.
		 */
		size_t size() const;

		/**
		 * @brief      Copy the array to a vector.
		 *
		 * @param      b     Target vector.
		 */
		void export_vector(std::vector<ArrValType> &b) const;

		/**
		 * @brief      Copy the array from a vector.
		 *
		 * @param      b     Target vector.
		 */
		void import_vector(const std::vector<ArrValType> &b);

#ifdef GCTL_EIGEN
		/**
		 * @brief      Copy the array to a Eigen3 vector.
		 *
		 * @param      b     Target vector.
		 */
		void export_eigen_vector(Eigen::VectorXd &b) const;

		/**
		 * @brief      Copy the array from a Eigen3 vector.
		 *
		 * @param      b     Target vector.
		 */
		void import_eigen_vector(const Eigen::VectorXd &b);
#endif // GCTL_EIGEN

		/**
		 * @brief element operate function pointer
		 * 
		 */
		typedef void (*foreach_a_ptr)(ArrValType *ele_ptr, size_t id);
		
		/**
		 * @brief    Operate on each and every element
		 * 
		 * @param func  operation function
		 */
		void for_each(foreach_a_ptr func);

		/**
		 * @brief   Display the elements.
		 * 
		 */
		void show(std::ostream &os = std::cout, char sep = ' ');

		/**
		 * @brief      Initialize the array with selected random types.
		 *
		 * @param[in]  np1          Mean (Gauss) or low bound value (Even)
		 * @param[in]  np2          Standard deviation (Gauss) or hig bound value (Even).
		 * @param[in]  mode         Random types. 'RdNormal' for Gaussian distributed numbers and 
		 * 'RdUniform' for even distributed numbers.
		 * @param[in] seed          Random seed.
		 */
		void random_float(ArrValType np1, ArrValType np2, random_type_e mode = RdNormal, unsigned int seed = 0);

		/**
		 * @brief      Initialize the array with selected random types.
		 *
		 * @param[in]  np1          Mean (Gauss) or low bound value (Even)
		 * @param[in]  np2          Standard deviation (Gauss) or hig bound value (Even).
		 * @param[in] seed          Random seed.
		 */
		void random_int(ArrValType np1, ArrValType np2, unsigned int seed = 0);

		/**
		 * @brief Set elements' value as a sequent.
		 * 
		 * @param st_val  Start value.
		 * @param inc Increasement.
		 */
		void sequence(ArrValType st_val, ArrValType inc);

		/**
		 * @brief Assign values as a 2D sequence.
		 * 
		 * @param rs Start value of the row sequence
		 * @param rinc Increase value of the row sequence
		 * @param cs Start value of the column sequence
		 * @param cinc Increase value of the column sequence
		 * @param rn Row nnumber
		 * @param cn Column number
		 */
		void sequence2d(ArrValType rs, ArrValType rinc, ArrValType cs, ArrValType cinc, 
			size_t rn, size_t cn);

		/**
		 * @brief Assign values as a 3D sequence.
		 * 
		 * @param ls Start value of the layer sequence
		 * @param linc Increase value of the layer sequence
		 * @param rs Start value of the row sequence
		 * @param rinc Increase value of the row sequence
		 * @param cs Start value of the column sequence
		 * @param cinc Increase value of the column sequence
		 * @param ln Layer number
		 * @param rn Row nnumber
		 * @param cn Column number
		 */
		void sequence3d(ArrValType ls, ArrValType linc, 
			ArrValType rs, ArrValType rinc, 
			ArrValType cs, ArrValType cinc, 
			size_t ln, size_t rn, size_t cn);

		/**
		 * @brief Scale elements
		 * 
		 * @param in_val input factor
		 */
		void scale(ArrValType in_val);

		/**
		 * @brief  Dot product of two arrays.
		 * 
		 * @param a Input array. Must be the same length_.
		 * @return dot product
		 */
		ArrValType dot(const array<ArrValType> &a);

		/**
		 * @brief Return the mean value.
		 * 
		 * @return mean value.
		 */
		ArrValType mean() const;

		/**
		 * @brief Return the variance value.
		 * 
		 * @return Variance value.
		 */
		ArrValType variance() const;

		/**
		 * @brief Return the standard deviation value.
		 * 
		 * @return std value.
		 */
		ArrValType std() const;

		/**
		 * @brief Return the unbiased standard deviation value.
		 * 
		 * @return std value.
		 */
		ArrValType std_unbiased() const;

		/**
		 * @brief Return the root mean square value.
		 * 
		 * @return RMS value.
		 */
		ArrValType rms() const;

		/**
		 * @brief Return the maximal value.
		 * 
		 * @return maximal value.
		 */
		ArrValType max() const;

		/**
		 * @brief Return the minimal value.
		 * 
		 * @return minimal value.
		 */
		ArrValType min() const;

		/**
		 * @brief      计算数组的模长
		 *
		 * @param[in]  n_type      模长的计算类型
		 * 
		 * @return     返回的模长
		 */
		ArrValType module(norm_type_e n_type = L2);

		/**
		 * @brief Normalize the array to the given module value
		 * 
		 * @param norm targeting norm-value
		 * @param n_type Norm-type
		 * 
		 * @return The module
		 */
		ArrValType normalize(ArrValType norm = 1.0, norm_type_e n_type = L2);

		/**
		 * @brief      计算一个向量相对于另一个向量的正交向量
		 *
		 * @param      b     输入的数组，计算后的正交向量以引用的形式返回
		 */
		void orth(array<ArrValType> &b);
		
		/**
		 * @brief Assign values as a power sequence.
		 * 
		 * @param base Power base
		 */
		void log2linear(const ArrValType &base);

		/**
		 * @brief Convert values to a log root.
		 * 
		 * @param base Log base
		 */
		void linear2log(ArrValType base);

		/**
		 * @brief Set values to a range.
		 * 
		 * @param min Lower bound
		 * @param max Higher bound
		 * @param rt Bounder type
		 * 
		 */
		void set2range(ArrValType min, ArrValType max, range_type_e rt = HardScale);
	};

	template <typename ArrValType>
	array<ArrValType>::array()
	{
		length_ = 0;
		val_ = nullptr;
	}

	template <typename ArrValType>
	array<ArrValType>::array(size_t len)
	{
		length_ = 0;
		val_ = nullptr;
		resize(len);
	}

	template <typename ArrValType>
	array<ArrValType>::array(size_t len, ArrValType init_val)
	{
		length_ = 0;
		val_ = nullptr;
		resize(len, init_val);
	}

	template <typename ArrValType>
	array<ArrValType>::array(ArrValType *init_array, size_t len)
	{
		length_ = 0;
		val_ = nullptr;
		resize(init_array, len);
	}

	template <typename ArrValType>
	array<ArrValType>::array(const array<ArrValType> &b)
	{
		length_ = 0;
		val_ = nullptr;
		resize(b);
	}

	template <typename ArrValType>
	array<ArrValType>::array(std::initializer_list<ArrValType> init_val)
	{
		length_ = 0;
		val_ = nullptr;
		resize(init_val);
	}

	template <typename ArrValType>
	array<ArrValType>::array(size_t s, ArrValType st, ArrValType inc)
	{
		length_ = 0;
		val_ = nullptr;
		resize(s, st, inc);
	}

	template <typename ArrValType>
	array<ArrValType> &array<ArrValType>::operator= (const array<ArrValType> &b)
	{
		resize(b.size());
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = b[i];
		}
		return *this;
	}

	template <typename ArrValType>
	array<ArrValType> &array<ArrValType>::operator= (ArrValType v)
	{
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = v;
		}
		return *this;
	}

	template <typename ArrValType>
	array<ArrValType> array<ArrValType>::operator+ (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator+] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE

		array<ArrValType> out(length_);
		for (size_t i = 0; i < length_; i++)
		{
			out[i] = val_[i] + b[i];
		}
		return out;
	}

	template <typename ArrValType>
	array<ArrValType> array<ArrValType>::operator- (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator-] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE

		array<ArrValType> out(length_);
		for (size_t i = 0; i < length_; i++)
		{
			out[i] = val_[i] - b[i];
		}
		return out;
	}

	template <typename ArrValType>
	array<ArrValType> array<ArrValType>::operator* (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator*] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE

		array<ArrValType> out(length_);
		for (size_t i = 0; i < length_; i++)
		{
			out[i] = val_[i] * b[i];
		}
		return out;
	}

	template <typename ArrValType>
	array<ArrValType> array<ArrValType>::operator/ (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator/] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE

		array<ArrValType> out(length_);
		for (size_t i = 0; i < length_; i++)
		{
			out[i] = val_[i] / b[i];
		}
		return out;
	}

	template <typename ArrValType>
	array<ArrValType>& array<ArrValType>::operator+= (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator+=] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE

		for (size_t i = 0; i < length_; i++)
		{
			val_[i] += b[i];
		}
		return *this;
	}

	template <typename ArrValType>
	array<ArrValType>& array<ArrValType>::operator-= (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator-=] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE

		for (size_t i = 0; i < length_; i++)
		{
			val_[i] -= b[i];
		}
		return *this;
	}

	template <typename ArrValType>
	array<ArrValType>& array<ArrValType>::operator*= (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator*=] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE

		for (size_t i = 0; i < length_; i++)
		{
			val_[i] *= b[i];
		}
		return *this;
	}

	template <typename ArrValType>
	array<ArrValType>& array<ArrValType>::operator/= (const array<ArrValType> &b)
	{
#ifdef GCTL_CHECK_SIZE
		if (b.size() != length_)
		{
			throw std::runtime_error("[gctl::array<T>::operator/=] Incompatible array sizes.");
		}
#endif // GCTL_CHECK_SIZE
		
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] /= b[i];
		}
		return *this;
	}

	template <typename ArrValType>
	array<ArrValType>::~array()
	{
		clear();
	}

	template <typename ArrValType>
	void array<ArrValType>::resize(size_t len)
	{
		if (len == 0)
		{
			throw std::invalid_argument("[gctl::array<T>::resize] Invalid array size.");
		}

		if (length_ != len)
		{
			clear();

			length_ = len;
			val_ = new ArrValType[len];
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::resize(size_t len, ArrValType init_val)
	{
		resize(len);
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = init_val;
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::resize(ArrValType *init_array, size_t len)
	{
		if (init_array == nullptr)
		{
			throw std::domain_error("[gctl::array<T>::resize] Invalid pointer.");
		}

		resize(len);
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = init_array[i];
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::resize(const array<ArrValType> &b)
	{
		if (!b.empty())
		{
			resize(b.size());
			for (size_t i = 0; i < length_; i++)
			{
				val_[i] = b[i];
			}
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::resize(std::initializer_list<ArrValType> init_val)
	{
		resize(init_val.size());

		typename std::initializer_list<ArrValType>::iterator iter = init_val.begin();
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = *iter;
			iter++;
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::resize(size_t s, ArrValType st, ArrValType inc)
	{
		resize(s);
		sequence(st, inc);
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::clear()
	{
		if (val_ != nullptr)
		{
			delete[] val_;
			val_ = nullptr;

			length_ = 0;
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::assign_all(ArrValType in_val)
	{
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = in_val;
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::assign(ArrValType in_val, size_t st, size_t range)
	{
		if (st + range <= length_)
		{
			for (size_t i = st; i < st + range; i++)
			{
				val_[i] = in_val;
			}
			return;
		}
		
		throw std::runtime_error("[gctl::array<T>::assign] Invalid segment range.");
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::append_array(const array<ArrValType> &b)
	{
		ArrValType *t = new ArrValType [length_ + b.size()];

		for (size_t i = 0; i < length_; i++)
		{
			t[i] = val_[i];
		}

		for (size_t i = 0; i < b.size(); i++)
		{
			t[i + length_] = b[i];
		}

		delete[] val_;
		val_ = t;

		length_ += b.size();
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::extract_array(array<ArrValType> &b, size_t st) const
	{
		if (b.size() + st <= length_)
		{
			for (size_t i = 0; i < b.size(); i++)
			{
				b[i] = val_[i + st];
			}
		}
		else
		{
			size_t c = 0;
			for (size_t i = st; i < length_; i++)
			{
				b[c] = val_[i];
				c++;
			}
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::extract_array(array<ArrValType> &b, const array<size_t> &ids) const
	{
		b.resize(ids.size());
		for (size_t i = 0; i < ids.size(); i++)
		{
			b[i] = val_[ids[i]];
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::insert_array(array<ArrValType> &b, size_t st)
	{
		if (b.size() + st <= length_)
		{
			for (size_t i = 0; i < b.size(); i++)
			{
				val_[i + st] = b[i];
			}
		}
		else
		{
			size_t c = 0;
			for (size_t i = st; i < length_; i++)
			{
				val_[i] = b[c];
				c++;
			}
		}
		return;
	}

	template <typename ArrValType>
	ArrValType *array<ArrValType>::get(size_t index) const
	{
#ifdef GCTL_CHECK_BOUNDER
		if (index >= length_)
			throw std::out_of_range("[gctl::array<T>::at] Invalid index.");
#endif // GCTL_CHECK_BOUNDER

		return &val_[index];
	}

	template <typename ArrValType>
	ArrValType &array<ArrValType>::at(size_t index)
	{
#ifdef GCTL_CHECK_BOUNDER
		if (index >= length_)
			throw std::out_of_range("[gctl::array<T>::at] Invalid index.");
#endif // GCTL_CHECK_BOUNDER

		return val_[index];
	}

	template <typename ArrValType>
	ArrValType &array<ArrValType>::at(size_t index) const
	{
#ifdef GCTL_CHECK_BOUNDER
		if (index >= length_)
			throw std::out_of_range("[gctl::array<T>::at] Invalid index.");
#endif // GCTL_CHECK_BOUNDER

		return val_[index];
	}

	template <typename ArrValType>
	ArrValType &array<ArrValType>::operator[](size_t index)
	{
		return val_[index];
	}

	template <typename ArrValType>
	ArrValType &array<ArrValType>::operator[](size_t index) const
	{
		return val_[index];
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::front() const
	{
		return val_[0];
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::back() const
	{
		return val_[length_-1];
	}

	template <typename ArrValType>
	iterator<ArrValType> array<ArrValType>::begin()
	{
		return iterator<ArrValType>(this->val_, 0);
	}

	template <typename ArrValType>
	iterator<ArrValType> array<ArrValType>::end()
	{
		return iterator<ArrValType>(this->val_, length_);
	}

	template <typename ArrValType>
	bool array<ArrValType>::empty() const
	{
		if (length_ == 0) return true;
		else return false;
	}

	template <typename ArrValType>
	size_t array<ArrValType>::size() const
	{
		return length_;
	}

	template <typename ArrValType>
	void array<ArrValType>::export_vector(std::vector<ArrValType> &b) const
	{
		if (!b.empty())
		{
			b.clear();
		}

		b.resize(length_);
		for (size_t i = 0; i < length_; i++)
		{
			b[i] = val_[i];
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::import_vector(const std::vector<ArrValType> &b)
	{
		resize(b.size());
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = b[i];
		}
		return;
	}

#ifdef GCTL_EIGEN
	template <typename ArrValType>
	void array<ArrValType>::export_eigen_vector(Eigen::VectorXd &b) const
	{
		b.resize(length_);
		for (size_t i = 0; i < length_; i++)
		{
			b[i] = val_[i];
		}
		return;
	}

	
	template <typename ArrValType>
	void array<ArrValType>::import_eigen_vector(const Eigen::VectorXd &b)
	{
		resize(b.size());
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = b[i];
		}
		return;
	}
#endif // GCTL_EIGEN

	template <typename ArrValType>
	void array<ArrValType>::for_each(foreach_a_ptr func)
	{
		for (size_t i = 0; i < length_; i++)
		{
			func(&val_[i], i);
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::show(std::ostream &os, char sep)
	{
		if (length_ != 0)
		{
			os << val_[0];
			for (size_t i = 1; i < length_; i++)
			{
				os << sep << val_[i];
			}
			os << std::endl;
		}
		return;
	}

	template <typename VecValType>
	std::ostream &operator <<(std::ostream &os, const array<VecValType> &a)
	{
		if (!a.empty())
		{
			os << a[0];
			for (size_t i = 1; i < a.size(); i++)
			{
				os << " " << a[i];
			}
		}
		return os;
	}

	template <typename VecValType>
	std::istream &operator >>(std::istream &os, array<VecValType> &a)
	{
		for (size_t i = 0; i < a.size(); i++)
		{
			os >> a[i];
		}
		return os;
	}

	template <typename ArrValType>
	void array<ArrValType>::random_float(ArrValType np1, ArrValType np2, 
		random_type_e mode, unsigned int seed)
	{
		static_assert(std::is_floating_point<ArrValType>::value, 
			"gctl::array<T>::random_float(...) could only be used with a float point type.");

		if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
		std::default_random_engine generator(seed);

		if (mode == RdNormal)
		{
			//添加高斯分布的随机值
			std::normal_distribution<ArrValType> dist(np1, np2);
			for (size_t i = 0; i < length_; i++)
			{
				val_[i] = dist(generator);
			}
			return;
		}

		//添加均匀分布的随机值
		std::uniform_real_distribution<ArrValType> dist(np1, np2);
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = dist(generator);
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::random_int(ArrValType np1, ArrValType np2, unsigned int seed)
	{
		static_assert(std::is_integral<ArrValType>::value, 
			"gctl::array<T>::random_int(...) could only be used with an integral type.");

		if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
		std::default_random_engine generator(seed);

		//添加均匀分布的随机值
		std::uniform_int_distribution<ArrValType> dist(np1, np2);
		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = dist(generator);
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::sequence(ArrValType st_val, ArrValType inc)
	{
		//static_assert(std::is_arithmetic<ArrValType>::value, 
		//	"gctl::array<T>::sequence(...) could only be used with an arithmetic type.");

		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = st_val + i*inc;
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::sequence2d(ArrValType rs, ArrValType rinc, 
		ArrValType cs, ArrValType cinc, size_t rn, size_t cn)
	{
		//static_assert(std::is_arithmetic<ArrValType>::value, 
		//	"gctl::array<T>::sequence2d(...) could only be used with an arithmetic type.");

		if (rn*cn != length_)
		{
			throw invalid_argument("[gctl::array<T>::sequence2d] Invalid sequence sizes.");
		}

		array<ArrValType> out_x(cn), out_y(rn);
		out_x.sequence(cs, cinc);
        out_y.sequence(rs, rinc);

		for (size_t i = 0; i < rn; i++)
        {
            for (size_t j = 0; j < cn; j++)
            {
                val_[j + i*cn] = out_x[j] + out_y[i];
            }
        }
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::sequence3d(ArrValType ls, ArrValType linc, 
		ArrValType rs, ArrValType rinc, ArrValType cs, ArrValType cinc, 
		size_t ln, size_t rn, size_t cn)
	{
		//static_assert(std::is_arithmetic<ArrValType>::value, 
		//	"gctl::array<T>::sequence3d(...) could only be used with an arithmetic type.");

		if (ln*rn*cn != length_)
		{
			throw invalid_argument("[gctl::array<T>::sequence3d] Invalid sequence sizes.");
		}

		array<ArrValType> out_x(cn), out_y(rn), out_z(ln);
		out_x.sequence(cs, cinc);
        out_y.sequence(rs, rinc);
		out_z.sequence(ls, linc);

		for (size_t k = 0; k < ln; k++)
        {
			for (size_t i = 0; i < rn; i++)
			{
				for (size_t j = 0; j < cn; j++)
				{
					val_[j + i*cn + k*rn*cn] = out_x[j] + out_y[i] + out_z[k];
				}
			}
		}
		return;
	}

	template <typename ArrValType>
	void array<ArrValType>::scale(ArrValType in_val)
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::scale(...) could only be used with an arithmetic type.");

		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = val_[i]*in_val;
		}
		return;
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::dot(const array<ArrValType> &a)
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::dot(...) could only be used with an arithmetic type.");

		if (length_ != a.size())
		{
			throw std::runtime_error("Incompatible array sizes. gctl::array<T>::dot(...)");
		}

		ArrValType sum = 0;
		for (size_t i = 0; i < length_; i++)
		{
			sum += val_[i]*a[i];
		}
		return sum;
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::mean() const
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::mean(...) could only be used with an arithmetic type.");

		if (length_ == 0) return 0;

		ArrValType m = 0;
		for (size_t i = 0; i < length_; i++)
		{
			m = m + val_[i];
		}
		return m/length_;
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::std() const
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::std(...) could only be used with an arithmetic type.");

		if (length_ == 0 || length_ == 1) return 0;

		ArrValType m = mean(), s = 0;
		for (size_t i = 0; i < length_; i++)
		{
			s += (val_[i] - m)*(val_[i] - m);
		}
		return sqrt(s/length_);
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::rms() const
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::rms(...) could only be used with an arithmetic type.");

		if (length_ == 0) return 0;

		ArrValType m = 0;
		for (size_t i = 0; i < length_; i++)
		{
			m += val_[i]*val_[i];
		}
		return sqrt(m/length_);
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::max() const
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::max(...) could only be used with an arithmetic type.");

		if (length_ == 0) return 0;

		ArrValType m = val_[0];
		for (size_t i = 1; i < length_; i++)
		{
			m = std::max(val_[i], m);
		}
		return m;
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::min() const
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::min(...) could only be used with an arithmetic type.");

		if (length_ == 0) return 0;

		ArrValType m = val_[0];
		for (size_t i = 1; i < length_; i++)
		{
			m = std::min(val_[i], m);
		}
		return m;
	}

	template <typename ArrValType>
	ArrValType array<ArrValType>::module(norm_type_e n_type)
    {
        static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::module(...) could only be used with an arithmetic type.");

        ArrValType norm_sum = 0.0;

		if (n_type == L0)
		{
			int n = 0;
			for (size_t i = 0; i < length_; i++)
			{
				if (val_[i] != 0.0)
				{
					n++;
				}
			}
			norm_sum = (ArrValType) n;
		}

		if (n_type == L1)
		{
			for (size_t i = 0; i < length_; i++)
			{
				norm_sum += GCTL_FABS(val_[i]);
			}
		}

		if (n_type == L2)
		{
			for (size_t i = 0; i < length_; i++)
			{
				norm_sum += val_[i] * val_[i];
			}
			norm_sum = sqrt(norm_sum);
		}

		if (n_type == Linf)
		{
			norm_sum = val_[0];
			for (int i = 1; i < length_; i++)
			{
				norm_sum = GCTL_MAX(norm_sum, val_[i]);
			}
		}

		return norm_sum;
    }

	template <typename ArrValType>
	ArrValType array<ArrValType>::normalize(ArrValType norm, norm_type_e n_type)
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::normalize(...) could only be used with an arithmetic type.");

		ArrValType norm_sum = module(n_type);

		if (norm_sum <= GCTL_ZERO) return 0;

		for (size_t i = 0; i < length_; i++)
		{
			val_[i] = val_[i]*norm/norm_sum;
		}
		return norm_sum;
	}

	template <typename ArrValType>
	void array<ArrValType>::orth(array<ArrValType> &b)
    {
        static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::orth(...) could only be used with an arithmetic type.");

        if (length_ != b.size())
        {
            throw runtime_error("[gctl::array<T>::orth] Incompatible array size.");
        }

        ArrValType product = dot(b);
        for (int i = 0; i < length_; i++)
        {
            b[i] = b[i] - product*val_[i];
        }
        return;
    }

    template <typename ArrValType>
	void array<ArrValType>::log2linear(const ArrValType &base)
    {
        static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::logspace(...) could only be used with an arithmetic type.");

        for (size_t i = 0; i < length_; i++)
        {
            val_[i] = pow(base, val_[i]);
        }
        return;
    }

    template <typename ArrValType>
	void array<ArrValType>::linear2log(ArrValType base)
    {
        static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::linear2log(...) could only be used with an arithmetic type.");

        for (size_t i = 0; i < length_; ++i)
        {
            val_[i] = log(val_[i])/log(base);
        }
        return;
    }

    template <typename ArrValType>
	ArrValType array<ArrValType>::variance() const
    {
        static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::variance(...) could only be used with an arithmetic type.");

		if (length_ == 0) return 0;

        ArrValType mn = mean();
        ArrValType d = 0;
        for (size_t i = 0; i < length_; i++)
        {
            d = d + (val_[i] - mn)*(val_[i] - mn);
        }
        return d/length_;
    }

    template <typename ArrValType>
	ArrValType array<ArrValType>::std_unbiased() const
    {
        static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::std_unbiased(...) could only be used with an arithmetic type.");

		if (length_ < 2) return 0;

        ArrValType mn = mean();
        ArrValType d = 0;
        for (size_t i = 0; i < length_; i++)
        {
            d = d + (val_[i] - mn)*(val_[i] - mn);
        }
        return sqrt(d/(length_ - 1));
    }

    template <typename ArrValType>
	void array<ArrValType>::set2range(ArrValType min, ArrValType max, range_type_e rt)
	{
		static_assert(std::is_arithmetic<ArrValType>::value, 
			"gctl::array<T>::set2range(...) could only be used with an arithmetic type.");

        ArrValType amin = val_[0], amax = val_[0];
        for (size_t i = 1; i < length_; i++)
        {
            amin = amin<val_[i]?amin:val_[i];
            amax = amax>val_[i]?amax:val_[i];
        }

        if (rt == HardScale)
        {
            for (size_t i = 0; i < length_; i++)
            {
                val_[i] = (max - min)*(val_[i] - amin)/(amax - amin) + min;
            }
            return;
        }

        if (amin >= min && amax <= max) return; // aleardy in the range, nothing to be done

        if (rt == CutOff)
        {
            for (size_t i = 0; i < length_; i++)
            {
                if (val_[i] > max) val_[i] = max;
                if (val_[i] < min) val_[i] = min;
            }
            return;
        }

        // Soft Scale
        ArrValType min2 = amin>min?amin:min;
        ArrValType max2 = amax<max?amax:max;
        for (size_t i = 0; i < length_; i++)
        {
            val_[i] = (max2 - min2)*(val_[i] - amin)/(amax - amin) + min2;
        }
        return;
	}
}

#endif // _GCTL_ARRAY_H