/// <summary>
		/// Constructor. Takes an read only vector and evaluates the spaces between
		/// the vector elements.
		/// </summary>
		/// <param name="vec">The vector.</param>
		public VectorSpacingEvaluator(IROVector vec)
		{
			_numtotalsteps = vec.Length - 1;
			for (int i = 0; i < _numtotalsteps; i++)
			{
				double step = vec[i + 1] - vec[i];

				if (!double.IsNaN(step))
				{
					_numvalidsteps++;

					if (step > _stepmax)
						_stepmax = step;
					if (step < _stepmin)
						_stepmin = step;

					_sumsteps += step;
				}
			}

			// if all steps are valid, we calculate sumsteps from the boundaries
			// to enhance the accuracy.
			if (_numvalidsteps > 0 && _numtotalsteps == _numvalidsteps)
				_sumsteps = vec[_numtotalsteps] - vec[0];
		}

		/// <summary>
		/// Constructor of a bivariate linear spline. The vectors and the data matrix are not cloned, so make sure that they don't change during usage of this instance.
		/// </summary>
		/// <param name="x">Vector of x values corresponding to the rows of the data matrix. Must be strongly increasing or decreasing.</param>
		/// <param name="y">Vector of y values corresponding to the columns of the data matrix. Must be strongly increasing or decreasing.</param>
		/// <param name="datamatrix"></param>
		public BivariateLinearSpline(IROVector x, IROVector y, IROMatrix datamatrix)
		{
			_x = x;
			_y = y;
			_vmatrix = datamatrix;

			// check the arguments
			if (_x.Length < 2)
				throw new ArgumentException("x.Length is less or equal 1 (you can use univariate interpolation instead)");
			if (_y.Length < 2)
				throw new ArgumentException("y.Length is less or equal 1 (you can use univariate interpolation instead)");
			if (_x.Length != _vmatrix.Rows)
				throw new ArgumentException("Length of vector x is not equal to datamatrix.Rows");
			if (_y.Length != _vmatrix.Columns)
				throw new ArgumentException("Length of vector y is not equal to datamatrix.Columns");

			if (!VectorMath.IsStrictlyIncreasingOrDecreasing(_x, out _isXDecreasing))
				throw new ArgumentException("Vector x is not strictly increasing or decreasing");

			if (!VectorMath.IsStrictlyIncreasingOrDecreasing(_y, out _isYDecreasing))
				throw new ArgumentException("Vector y is not strictly increasing or decreasing");

			_lastIX = 0;
			_lastIY = 0;
		}

    /// <summary>
    /// Constructor. Takes an read only vector and evaluates the spaces between
    /// the vector elements.
    /// </summary>
    /// <param name="vec">The vector.</param>
    public VectorSpacingEvaluator(IROVector vec)
    {
      int lower = vec.LowerBound;
      int upper = vec.UpperBound;

      _numtotalsteps = upper-lower;
      for(int i=lower;i<upper;i++)
      {
        double step = vec[i+1]-vec[i];
        
        if(!double.IsNaN(step))
        {
          _numvalidsteps++;

          if(step>_stepmax)
            _stepmax = step;
          if(step<_stepmin)
            _stepmin = step;

          _sumsteps += step;
        }
      }

      // if all steps are valid, we calculate sumsteps from the boundaries
      // to enhance the accuracy.
      if(_numvalidsteps>0 && _numtotalsteps == _numvalidsteps)
        _sumsteps = vec[upper] - vec[lower];
    }

		public static void SetContentFromMatrix(DataTable destinationTable, IROMatrix matrix, string columnBaseName, IROVector rowHeaderColumn, string rowHeaderColumnName, IROVector colHeaderColumn, string colHeaderColumnName)
		{
			var c = new MatrixToDataTableConverter(matrix, destinationTable);
			c.ColumnBaseName = columnBaseName;
			c.AddMatrixColumnHeaderData(rowHeaderColumn, rowHeaderColumnName);
			c.AddMatrixColumnHeaderData(colHeaderColumn, colHeaderColumnName);
			c.Execute();
		}

		private static int FindIndex(IROVector v, bool isDecreasing, int lastIdx, double x)
		{
			if (isDecreasing) // strictly decreasing
			{
				if (x > v[lastIdx])
				{
					if (lastIdx == 0)
						return -1;
					if (x <= v[lastIdx - 1])
						return lastIdx - 1;
					return BinarySearchForIndex(v, isDecreasing, x);
				}
				else if (x < v[lastIdx + 1])
				{
					if (lastIdx + 2 <= v.Length)
						return -1;
					if (x >= v[lastIdx + 2])
						return lastIdx + 1;
					return BinarySearchForIndex(v, isDecreasing, x);
				}
				else
				{
					return lastIdx;
				}
			}
			else // strictly increasing
			{
				if (x < v[lastIdx])
				{
					if (lastIdx == 0)
						return -1;
					if (x >= v[lastIdx - 1])
						return lastIdx - 1;
					return BinarySearchForIndex(v, isDecreasing, x);
				}
				else if (x > v[lastIdx + 1])
				{
					if (lastIdx + 2 >= v.Length)
						return -1;
					if (x <= v[lastIdx + 2])
						return lastIdx + 1;
					return BinarySearchForIndex(v, isDecreasing, x);
				}
				else
				{
					return lastIdx;
				}
			}
		}

		public static double Nrd0(IROVector x)
		{
			if (x.Length < 2) throw new ArgumentException("need at least 2 data points");

			double hi = Statistics.StandardDeviation(x);
			double lo = Math.Min(hi, Statistics.InterQuartileRange(x) / 1.34);  // qnorm(.75) - qnorm(.25) = 1.34898
			if (lo.IsNaN())
			{
				lo = hi;
				if (lo.IsNaN())
				{
					lo = Math.Abs(x[0]);
					if (lo.IsNaN())
						lo = 1;
				}
			}

			return 0.9 * lo * Math.Pow(x.Length, (-0.2));
		}

		/// <summary>
		/// Constructs an Akima bivariate spline.
		/// </summary>
		/// <param name="x">ARRAY OF DIMENSION LX STORING THE X COORDINATES OF INPUT GRID POINTS (IN ASCENDING ORDER)</param>
		/// <param name="y">ARRAY OF DIMENSION LY STORING THE Y COORDINATES OF INPUT GRID POINTS (IN ASCENDING ORDER)</param>
		/// <param name="z">DOUBLY-DIMENSIONED ARRAY OF DIMENSION (LX,LY) STORING THE VALUES OF THE FUNCTION (Z VALUES) AT INPUT GRID POINTS</param>
		/// <param name="copyDataLocally">If true, the data where cloned before stored here in this instance. If false, the data
		/// are stored directly. Make sure then, that the data are not changed outside.</param>
		public BivariateAkimaSpline(IROVector x, IROVector y, IROMatrix z, bool copyDataLocally)
		{
			if (copyDataLocally)
			{
				_myX = VectorMath.ToVector(new double[x.Length]);
				VectorMath.Copy(x, (IVector)_myX);

				_myY = VectorMath.ToVector(new double[y.Length]);
				VectorMath.Copy(y, (IVector)_myY);

				_myZ = new MatrixMath.BEMatrix(_myZ.Rows, _myZ.Columns);
				MatrixMath.Copy(z, (IMatrix)_myZ);
			}
			else
			{
				_myX = x;
				_myY = y;
				_myZ = z;
			}
		}

      public void Append(IROVector a)
      {
        if(_length+a.Length>=_arr.Length)
          Redim((int)(32+1.3*(_length+a.Length)));


        for(int i=0;i<a.Length;i++)
          _arr[i+_length] = a[i+a.LowerBound];
        _length += a.Length;
      }

 /// <summary>
 /// Wraps a section of a original vector <c>x</c> into a new vector.
 /// </summary>
 /// <param name="x">Original vector.</param>
 /// <param name="start">Index of the start of the section to wrap.</param>
 /// <param name="len">Length (=number of elements) of the section to wrap.</param>
 /// <returns>A IROVector that contains the section from <c>start</c> to <c>start+len-1</c> of the original vector.</returns>
 public static IROVector ToROVector(IROVector x, int start, int len)
 {
   return new ROVectorSectionWrapper(x, start, len);
 }

 /// <summary>
 /// Processes the spectra in matrix xMatrix.
 /// </summary>
 /// <param name="xMatrix">The matrix of spectra. Each spectrum is a row of the matrix.</param>
 /// <param name="xMean">Output: On return, contains the ensemble mean of the spectra.</param>
 /// <param name="xScale">Not used.</param>
 /// <param name="regions">Vector of spectal regions. Each element is the index of the start of a new region.</param>
 public virtual void ProcessForPrediction(IMatrix xMatrix, IROVector xMean, IROVector xScale, int[] regions)
 {
 }

 /// <summary>
 /// Constructor, takes a double array for wrapping.
 /// </summary>
 /// <param name="x"></param>
 /// <param name="start">Start index of the section to wrap.</param>
 /// <param name="len">Length of the section to wrap.</param>
 public ROVectorSectionWrapper(IROVector x, int start, int len)
 {
   if(start>=x.Length)
     throw new ArgumentException("Start of the section is beyond length of the vector");
  if (start+len>=x.Length)
     throw new ArgumentException("End of the section is beyond length of the vector");
  
   _x = x;
   _start = start;
   _length = len;
 }

    /// <summary>
    /// Returns the sum of squared differences of the elements of xarray and yarray.
    /// </summary>
    /// <param name="xarray">The first array.</param>
    /// <param name="yarray">The other array.</param>
    /// <returns>The sum of squared differences all elements of xarray and yarray.</returns>
    public static double SumOfSquaredDifferences(IROVector xarray, IROVector yarray)
    {
      if(xarray.Length!=yarray.Length)
        throw new ArgumentException("Length of xarray is unequal length of yarray");

      double sum = 0;
      for(int i=0;i<xarray.Length;i++)
        sum += Square(xarray[i]-yarray[i]);

      return sum;
    }

		public virtual void StoreXOfX(IROVector xOfX, DataTable table)
		{
			DoubleColumn xColOfX = new DoubleColumn();
			VectorMath.Copy(xOfX, DataColumnWrapper.ToVector(xColOfX, xOfX.Length));
			table.DataColumns.Add(xColOfX, _XOfX_ColumnName, Altaxo.Data.ColumnKind.X, 0);
		}

		public void SetErrorVariance(IROVector dyy, double errvar)
		{
			dy.CopyFrom(dyy);
			var = errvar;
		}

    /// <summary>
    /// Adds (elementwise) two vectors a and b and stores the result in c. All vectors must have the same length.
    /// </summary>
    /// <param name="a">First summand.</param>
    /// <param name="b">Second summand.</param>
    /// <param name="c">The resulting vector.</param>
    public static void Add(IROVector a, IROVector b, IVector c)
    {
      if(a.Length != b.Length)
        throw new ArgumentException("Length of vectors a and b unequal");
      if(c.Length != b.Length)
        throw new ArgumentException("Length of vectors a and c unequal");
      if(a.LowerBound != b.LowerBound || a.LowerBound != c.LowerBound)
        throw new ArgumentException("Vectors a, b, and c have not the same LowerBound property");

      int end = c.UpperBound;
      for(int i=c.LowerBound;i<=end;i++)
        c[i]=a[i]+b[i];
    }

		/// <summary>
		///
		/// </summary>
		/// <param name="mcalib"></param>
		/// <param name="groupingStrategy"></param>
		/// <param name="preprocessOptions"></param>
		/// <param name="xOfX"></param>
		/// <param name="matrixX">Matrix of horizontal spectra, centered and preprocessed.</param>
		/// <param name="matrixY">Matrix of concentrations, centered.</param>
		/// <param name="numberOfFactors"></param>
		/// <param name="predictedY"></param>
		/// <param name="spectralResiduals"></param>
		public virtual void CalculateCrossPredictedY(
			IMultivariateCalibrationModel mcalib,
			ICrossValidationGroupingStrategy groupingStrategy,
			SpectralPreprocessingOptions preprocessOptions,
			IROVector xOfX,
			IMatrix matrixX,
			IMatrix matrixY,
			int numberOfFactors,
			IMatrix predictedY,
			IMatrix spectralResiduals)
		{
			MultivariateRegression.GetCrossYPredicted(xOfX,
				matrixX, matrixY, numberOfFactors, groupingStrategy, preprocessOptions,
				this.CreateNewRegressionObject(),
				predictedY);
		}

		/// <summary>
		/// Get the matrix of x and y values (raw data).
		/// </summary>
		/// <param name="srctable">The table where the data come from.</param>
		/// <param name="selectedColumns">The selected columns.</param>
		/// <param name="selectedRows">The selected rows.</param>
		/// <param name="selectedPropertyColumns">The selected property column(s).</param>
		/// <param name="bHorizontalOrientedSpectrum">True if a spectrum is a single row, False if a spectrum is a single column.</param>
		/// <param name="matrixX">On return, gives the matrix of spectra (each spectra is a row in the matrix).</param>
		/// <param name="matrixY">On return, gives the matrix of y-values (each measurement is a row in the matrix).</param>
		/// <param name="plsContent">Holds information about the analysis results.</param>
		/// <param name="xOfX">On return, this is the vector of values corresponding to each spectral bin, i.e. wavelength values, frequencies etc.</param>
		/// <returns></returns>
		public static string GetXYMatrices(
			Altaxo.Data.DataTable srctable,
			IAscendingIntegerCollection selectedColumns,
			IAscendingIntegerCollection selectedRows,
			IAscendingIntegerCollection selectedPropertyColumns,
			bool bHorizontalOrientedSpectrum,
			MultivariateContentMemento plsContent,
			out IMatrix matrixX,
			out IMatrix matrixY,
			out IROVector xOfX
			)
		{
			matrixX = null;
			matrixY = null;
			xOfX = null;
			plsContent.SpectrumIsRow = bHorizontalOrientedSpectrum;

			Altaxo.Data.DataColumn xColumnOfX = null;
			Altaxo.Data.DataColumn labelColumnOfX = new Altaxo.Data.DoubleColumn();

			Altaxo.Data.DataColumnCollection concentration = bHorizontalOrientedSpectrum ? srctable.DataColumns : srctable.PropertyColumns;

			// we presume for now that the spectrum is horizontally,
			// if not we exchange the collections later

			AscendingIntegerCollection numericDataCols = new AscendingIntegerCollection();
			AscendingIntegerCollection numericDataRows = new AscendingIntegerCollection();
			AscendingIntegerCollection concentrationIndices = new AscendingIntegerCollection();

			AscendingIntegerCollection spectralIndices = bHorizontalOrientedSpectrum ? numericDataCols : numericDataRows;
			AscendingIntegerCollection measurementIndices = bHorizontalOrientedSpectrum ? numericDataRows : numericDataCols;

			plsContent.ConcentrationIndices = concentrationIndices;
			plsContent.MeasurementIndices = measurementIndices;
			plsContent.SpectralIndices = spectralIndices;
			plsContent.SpectrumIsRow = bHorizontalOrientedSpectrum;
			plsContent.OriginalDataTableName = srctable.Name;

			bool bUseSelectedColumns = (null != selectedColumns && 0 != selectedColumns.Count);
			// this is the number of columns (for now), but it can be less than this in case
			// not all columns are numeric
			int prenumcols = bUseSelectedColumns ? selectedColumns.Count : srctable.DataColumns.ColumnCount;
			// check for the number of numeric columns
			int numcols = 0;
			for (int i = 0; i < prenumcols; i++)
			{
				int idx = bUseSelectedColumns ? selectedColumns[i] : i;
				if (srctable[idx] is Altaxo.Data.INumericColumn)
				{
					numericDataCols.Add(idx);
					numcols++;
				}
			}

			// check the number of rows
			bool bUseSelectedRows = (null != selectedRows && 0 != selectedRows.Count);
			int numrows;
			if (bUseSelectedRows)
			{
				numrows = selectedRows.Count;
				numericDataRows.Add(selectedRows);
			}
			else
			{
				numrows = 0;
				for (int i = 0; i < numcols; i++)
				{
					int idx = bUseSelectedColumns ? selectedColumns[i] : i;
					numrows = Math.Max(numrows, srctable[idx].Count);
				}
				numericDataRows.Add(ContiguousIntegerRange.FromStartAndCount(0, numrows));
			}

			if (bHorizontalOrientedSpectrum)
			{
				if (numcols < 2)
					return "At least two numeric columns are neccessary to do Partial Least Squares (PLS) analysis!";

				// check that the selected columns are in exactly two groups
				// the group which has more columns is then considered to have
				// the spectrum, the other group is the y-values
				int group0 = -1;
				int group1 = -1;
				int groupcount0 = 0;
				int groupcount1 = 0;

				for (int i = 0; i < numcols; i++)
//.........这里部分代码省略.........

		/// <summary>
		/// Calculate the cross PRESS values and stores the results in the provided table.
		/// </summary>
		/// <param name="xOfX">Vector of spectral wavelengths. Necessary to divide the spectras in different regions.</param>
		/// <param name="matrixX">Matrix of spectra (horizontal oriented).</param>
		/// <param name="matrixY">Matrix of concentrations.</param>
		/// <param name="plsOptions">Analysis options.</param>
		/// <param name="plsContent">Information about this analysis.</param>
		/// <param name="table">Table to store the results.</param>
		public virtual void CalculateCrossPRESS(
			IROVector xOfX,
			IMatrix matrixX,
			IMatrix matrixY,
			MultivariateAnalysisOptions plsOptions,
			MultivariateContentMemento plsContent,
			DataTable table
			)
		{
			IROVector crossPRESSMatrix;

			Altaxo.Data.DoubleColumn crosspresscol = new Altaxo.Data.DoubleColumn();

			double meanNumberOfExcludedSpectra = 0;
			if (plsOptions.CrossPRESSCalculation != CrossPRESSCalculationType.None)
			{
				// now a cross validation - this can take a long time for bigger matrices

				MultivariateRegression.GetCrossPRESS(
					xOfX, matrixX, matrixY, plsOptions.MaxNumberOfFactors, GetGroupingStrategy(plsOptions),
					plsContent.SpectralPreprocessing,
					this.CreateNewRegressionObject(),
					out crossPRESSMatrix);

				VectorMath.Copy(crossPRESSMatrix, DataColumnWrapper.ToVector(crosspresscol, crossPRESSMatrix.Length));

				table.DataColumns.Add(crosspresscol, GetCrossPRESSValue_ColumnName(), Altaxo.Data.ColumnKind.V, 4);

				plsContent.MeanNumberOfMeasurementsInCrossPRESSCalculation = plsContent.NumberOfMeasurements - meanNumberOfExcludedSpectra;
			}
			else
			{
				table.DataColumns.Add(crosspresscol, GetCrossPRESSValue_ColumnName(), Altaxo.Data.ColumnKind.V, 4);
			}
		}

		/// <summary>
		/// Execute an analysis and stores the result in the provided table.
		/// </summary>
		/// <param name="matrixX">The matrix of spectra (horizontal oriented), centered and preprocessed.</param>
		/// <param name="matrixY">The matrix of concentrations, centered.</param>
		/// <param name="plsOptions">Information how to perform the analysis.</param>
		/// <param name="plsContent">A structure to store information about the results of the analysis.</param>
		/// <param name="table">The table where to store the results to.</param>
		/// <param name="press">On return, gives a vector holding the PRESS values of the analysis.</param>
		public virtual void ExecuteAnalysis(
			IMatrix matrixX,
			IMatrix matrixY,
			MultivariateAnalysisOptions plsOptions,
			MultivariateContentMemento plsContent,
			DataTable table,
			out IROVector press
			)
		{
			int numFactors = Math.Min(matrixX.Columns, plsOptions.MaxNumberOfFactors);
			MultivariateRegression regress = this.CreateNewRegressionObject();
			regress.AnalyzeFromPreprocessed(matrixX, matrixY, numFactors);
			plsContent.NumberOfFactors = regress.NumberOfFactors;
			plsContent.CrossValidationType = plsOptions.CrossPRESSCalculation;
			press = regress.GetPRESSFromPreprocessed(matrixX);

			Import(regress.CalibrationModel, table);
		}

		public virtual void StorePRESSData(
			IROVector PRESS,
			DataTable table)
		{
			StoreNumberOfFactors(PRESS.Length, table);

			Altaxo.Data.DoubleColumn presscol = new Altaxo.Data.DoubleColumn();
			for (int i = 0; i < PRESS.Length; i++)
				presscol[i] = PRESS[i];
			table.DataColumns.Add(presscol, GetPRESSValue_ColumnName(), Altaxo.Data.ColumnKind.V, 4);
		}