/// <summary>
		/// Execution of the fast nonnegative least squares algorithm. The algorithm finds a vector x with all elements xi&gt;=0 which minimizes |X*x-y|.
		/// </summary>
		/// <param name="XtX">X transposed multiplied by X, thus a square matrix.</param>
		/// <param name="Xty">X transposed multiplied by Y, thus a matrix with one column and same number of rows as X.</param>
		/// <param name="tolerance">Used to decide if a solution element is less than or equal to zero. If this is null, a default tolerance of tolerance = MAX(SIZE(XtX)) * NORM(XtX,1) * EPS is used.</param>
		/// <param name="x">Output: solution vector (matrix with one column and number of rows according to dimension of X.</param>
		/// <param name="w">Output: Lagrange vector. Elements which take place in the fit are set to 0. Elements fixed to zero contain a negative number.</param>
		/// <remarks>
		/// <para>
		/// Literature: Rasmus Bro and Sijmen De Jong, 'A fast non-negativity-constrained least squares algorithm', Journal of Chemometrics, Vol. 11, 393-401 (1997)
		/// </para>
		/// <para>
		/// Algorithm modified by Dirk Lellinger 2015 to allow a mixture of restricted and unrestricted parameters.
		/// </para>
		/// </remarks>
		public static void Execution(IROMatrix XtX, IROMatrix Xty, double? tolerance, out IMatrix x, out IMatrix w)
		{
			Execution(XtX, Xty, (i) => true, tolerance, out x, out w);
		}

    ///<summary>Solves a system on linear equations, AX=B, where A is the factored matrixed.</summary>
    ///<param name="B">RHS side of the system.</param>
    ///<returns>the solution matrix, X.</returns>  
    ///<exception cref="ArgumentNullException">B is null.</exception>
    ///<exception cref="NotPositiveDefiniteException">A is not positive definite.</exception>
    ///<exception cref="ArgumentException">The number of rows of A and B must be the same.</exception>
    public DoubleMatrix Solve (IROMatrix B) 
    {
      if ( B == null ) 
      {
        throw new System.ArgumentNullException("B cannot be null.");
      }
      Compute();
      if ( !ispd ) 
      {
        throw new NotPositiveDefiniteException();
      } 
      else 
      {
        if ( B.Rows != order ) 
        {
          throw new System.ArgumentException("Matrix row dimensions must agree." );
        }
#if MANAGED
        // Copy right hand side.
        int cols = B.Columns;
        DoubleMatrix X = new DoubleMatrix(B);
        for (int c = 0; c < cols; c++ ) 
        {
          // Solve L*Y = B;
          for (int i = 0; i < order; i++) 
          {
            double sum = B[i,c];
            for (int k = i-1; k >= 0; k--) 
            {
              sum -= l.data[i][k] * X.data[k][c];
            }
            X.data[i][c] = sum / l.data[i][i];
          }

          // Solve L'*X = Y;
          for (int i =order-1; i >= 0; i--) 
          {
            double sum = X.data[i][c];
            for (int k = i+1; k < order; k++) 
            {
              sum -= l.data[k][i] * X.data[k][c];
            }
            X.data[i][c] = sum / l.data[i][i];
          }
        }

        return X;
#else
                double[] rhs = DoubleMatrix.ToLinearArray(B);
                Lapack.Potrs.Compute(Lapack.UpLo.Lower,order,B.Columns,l.data,order,rhs,B.Rows);
                DoubleMatrix ret = new DoubleMatrix(order,B.Columns);
                ret.data = rhs;
                return ret;
#endif
      }
    }

		/// <summary>
		/// Creates an analyis from preprocessed spectra and preprocessed concentrations.
		/// </summary>
		/// <param name="matrixX">The spectral matrix (each spectrum is a row in the matrix). They must at least be centered.</param>
		/// <param name="matrixY">The matrix of concentrations (each experiment is a row in the matrix). They must at least be centered.</param>
		/// <param name="maxFactors">Maximum number of factors for analysis.</param>
		/// <returns>A regression object, which holds all the loads and weights neccessary for further calculations.</returns>
		protected override void AnalyzeFromPreprocessedWithoutReset(IROMatrix matrixX, IROMatrix matrixY, int maxFactors)
		{
			int numberOfFactors = _calib.NumberOfFactors = Math.Min(matrixX.Columns, maxFactors);
			IMatrix helperY = new MatrixMath.BEMatrix(matrixY.Rows, 1);

			_PRESS = null;

			for (int i = 0; i < matrixY.Columns; i++)
			{
				MatrixMath.Submatrix(matrixY, helperY, 0, i);

				PLS2Regression r = PLS2Regression.CreateFromPreprocessed(matrixX, helperY, maxFactors);

				IPLS2CalibrationModel cal = r.CalibrationModel;
				_calib.NumberOfFactors = Math.Min(_calib.NumberOfFactors, cal.NumberOfFactors);
				_calib.XLoads[i] = cal.XLoads;
				_calib.YLoads[i] = cal.YLoads;
				_calib.XWeights[i] = cal.XWeights;
				_calib.CrossProduct[i] = cal.CrossProduct;

				if (_PRESS == null)
					_PRESS = VectorMath.CreateExtensibleVector(r.PRESS.Length);
				VectorMath.Add(_PRESS, r.PRESS, _PRESS);
			}
		}

		/// <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;
		}

		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();
		}

		public void Decompose(IROMatrix A)
		{
			// Initialize.
			if (m == A.Rows && n == A.Columns)
			{
				MatrixMath.Copy(A, new JaggedArrayMatrix(QR, m, n));
				//JaggedArrayMath.Copy(A, QR);
			}
			else
			{
				QR = JaggedArrayMath.GetMatrixCopy(A);
				m = A.Rows;
				n = A.Columns;
				Rdiag = new double[n];
			}

			// Main loop.
			for (int k = 0; k < n; k++)
			{
				// Compute 2-norm of k-th column without under/overflow.
				double nrm = 0;
				for (int i = k; i < m; i++)
				{
					nrm = RMath.Hypot(nrm, QR[i][k]);
				}

				if (nrm != 0.0)
				{
					// Form k-th Householder vector.
					if (QR[k][k] < 0)
					{
						nrm = -nrm;
					}
					for (int i = k; i < m; i++)
					{
						QR[i][k] /= nrm;
					}
					QR[k][k] += 1.0;

					// Apply transformation to remaining columns.
					for (int j = k + 1; j < n; j++)
					{
						double s = 0.0;
						for (int i = k; i < m; i++)
						{
							s += QR[i][k] * QR[i][j];
						}
						s = -s / QR[k][k];
						for (int i = k; i < m; i++)
						{
							QR[i][j] += s * QR[i][k];
						}
					}
				}
				Rdiag[k] = -nrm;
			}
		}

		public MatrixToDataTableConverter(IROMatrix sourceMatrix, DataTable destinationTable)
		{
			if (null == sourceMatrix)
				throw new ArgumentNullException("sourceMatrix");
			if (null == destinationTable)
				throw new ArgumentNullException("destinationTable");

			_sourceMatrix = sourceMatrix;
			_destinationTable = destinationTable;
		}

		/// <summary>
		/// Allocates an array of the same dimensions than the provided matrix a, and copies the element of a to the new array.
		/// </summary>
		/// <param name="a">The matrix to copy.</param>
		/// <returns>New jagged array with the same element values than matrix a.</returns>
		public static double[][] GetMatrixCopy(IROMatrix a)
		{
			int rows = a.Rows;
			int cols = a.Columns;
			double[][] result = new double[rows][];
			double[] row;
			for (int i = 0; i < rows; i++)
			{
				result[i] = row = new double[cols];
				for (int j = 0; j < cols; j++)
					row[j] = a[i, j];
			}
			return result;
		}

		public override IROVector GetPRESSFromPreprocessed(IROMatrix matrixX)
		{
			IROVector result;
			CalculatePRESS(
				matrixX,
				_calib.XLoads,
				_calib.YLoads,
				_calib.XScores,
				_calib.CrossProduct,
				_calib.NumberOfFactors,
				out result);

			return result;
		}

    ///<summary>Constructor for Cholesky decomposition class. The constructor performs the factorization of a symmetric positive
    ///definite matrax and the Cholesky factored matrix is accessible by the <c>Factor</c> property. The factor is the lower 
    ///triangular factor.</summary>
    ///<param name="matrix">The matrix to factor.</param>
    ///<exception cref="ArgumentNullException">matrix is null.</exception>
    ///<exception cref="NotSquareMatrixException">matrix is not square.</exception>
    ///<remarks>This class only uses the lower triangle of the input matrix. It ignores the
    ///upper triangle.</remarks>
    public DoubleCholeskyDecomp(IROMatrix matrix)
    {
      if ( matrix == null ) 
      {
        throw new System.ArgumentNullException("matrix cannot be null.");
      }

      if ( matrix.Rows != matrix.Columns ) 
      {
        throw new NotSquareMatrixException("Matrix must be square.");
      }

      order = matrix.Columns;
      this.matrix = new DoubleMatrix(matrix);
    }

    /// <summary>
    /// <see cref="ICrossValidationGroupingStrategy.Group" />
    /// </summary>
    /// <param name="Y"></param>
    /// <returns></returns>
    public int[][] Group(IROMatrix Y)
    {
      System.Collections.ArrayList groups = new System.Collections.ArrayList();

      // add the first y-row to the first group
      System.Collections.ArrayList newcoll = new System.Collections.ArrayList();
      newcoll.Add(0);
      groups.Add(newcoll);
      // now test all other rows of the y-matrix against the existing groups
      for(int i=1;i<Y.Rows;i++)
      {
        bool bNewGroup=true;
        for(int gr=0;gr<groups.Count;gr++)
        {
          int refrow = (int)(((System.Collections.ArrayList)groups[gr])[0]);
          bool match = true;
          for(int j=0;j<Y.Columns;j++)
          {
            if(Y[i,j]!= Y[refrow,j])
            {
              match=false;
              break;
            }
          }
            
          if(match)
          {
            bNewGroup=false;
            ((System.Collections.ArrayList)groups[gr]).Add(i);
            break;
          }
        }
        if(bNewGroup)
        {
          newcoll = new System.Collections.ArrayList();
          newcoll.Add(i);
          groups.Add(newcoll);
        }
      }

      int[][] result = new int[groups.Count][];
      for(int i=0;i<result.Length;i++)
        result[i] = (int[])((System.Collections.ArrayList)groups[i]).ToArray(typeof(int));
      return result;
    }

		/// <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;
			}
		}

    /// <summary>
    /// Calculates the spectral residuals obtained from cross validation.
    /// </summary>
    /// <param name="spectralRegions">Array of ascending indices representing the starting indices of spectral regions.</param>
    /// <param name="X">Matrix of spectra (a spectrum = a row in the matrix).</param>
    /// <param name="Y">Matrix of y values (e.g. concentrations).</param>
    /// <param name="numFactors">Number of factors used for calculation.</param>
    /// <param name="groupingStrategy">The strategy how to group the spectra for cross prediction.</param>
    /// <param name="preprocessOptions">Information how to preprocess the data.</param>
    /// <param name="regress">The type of regression (e.g. PCR, PLS1, PLS2) provided as an empty regression object.</param>
    /// <param name="crossXResiduals">Returns the matrix of spectral residuals</param>
    /// <returns>Mean number of spectra used for prediction.</returns>
    public static double GetCrossXResiduals(
      int[] spectralRegions,
      IROMatrix X, // matrix of spectra (a spectra is a row of this matrix)
      IROMatrix Y, // matrix of concentrations (a mixture is a row of this matrix)
      int numFactors,
      ICrossValidationGroupingStrategy groupingStrategy,
      SpectralPreprocessingOptions preprocessOptions,
      MultivariateRegression regress,

      out IROMatrix crossXResiduals
      )
    {
      CrossPredictedXResidualsEvaluator worker = new CrossPredictedXResidualsEvaluator(X.Rows,spectralRegions,numFactors,groupingStrategy,preprocessOptions,regress);
      double result = CrossValidationIteration(X,Y,groupingStrategy,new CrossValidationIterationFunction(worker.EhCrossValidationWorker));
      crossXResiduals = worker.XCrossResiduals;
      return result;
    }

 protected override void InternalGetXLeverageFromPreprocessed(IROMatrix matrixX, int numFactors, IMatrix xLeverage)
 {
   PLS2Regression.CalculateXLeverageFromPreprocessed(matrixX,_calib.XWeights,numFactors,xLeverage,0);
 }

    /// <summary>
    /// Fits a data set linear to a given x base.
    /// </summary>
    /// <param name="xbase">The matrix of x values of the data set. Dimensions: numberOfData x numberOfParameters. The matrix is changed during calculation!</param>
    /// <param name="yarr">The array of y values of the data set.</param>
    /// <param name="stddev">The array of y standard deviations of the data set. Can be null if the standard deviation is unkown.</param>
    /// <param name="numberOfData">The number of data points (may be smaller than the array sizes of the data arrays).</param>
    /// <param name="numberOfParameter">The number of parameters to fit == size of the function base.</param>
    /// <param name="threshold">A treshold value (usually 1E-5) used to chop the unimportant singular values away.</param>
    public LinearFitBySvd Calculate(
      IROMatrix xbase, // NumberOfData, NumberOfParameters 
      double[] yarr,
      double[] stddev,
      int numberOfData,
      int numberOfParameter,
      double threshold)
    {
      _numberOfParameter = numberOfParameter;
      _numberOfFreeParameter = numberOfParameter;
      _numberOfData      = numberOfData;
      _parameter = new double[numberOfParameter];
      _residual = new double[numberOfData];
      _predicted = new double[numberOfData];
      _reducedPredictionVariance = new double[numberOfData];

      double[] scaledY      = new double[numberOfData];
      
      // Calculated some useful values
      _yMean = Mean(yarr,0,_numberOfData);
      _yCorrectedSumOfSquares = CorrectedSumOfSquares(yarr,_yMean,0,_numberOfData);

      MatrixMath.BEMatrix u = new MatrixMath.BEMatrix(numberOfData,numberOfParameter);
      // Fill the function base matrix (rows: numberOfData, columns: numberOfParameter)
      // and scale also y
      if (null == stddev)
      {
        for (int i = 0; i < numberOfData; i++)
        {
          for (int j = 0; j < numberOfParameter; j++)
            u[i, j] = xbase[i, j];

          scaledY[i] = yarr[i];
        }
      }
      else
      {
        for (int i = 0; i < numberOfData; i++)
        {
          double scale = 1 / stddev[i];

          for (int j = 0; j < numberOfParameter; j++)
            u[i, j] = scale * xbase[i, j];

          scaledY[i] = scale * yarr[i];
        }
      }
      _decomposition = MatrixMath.GetSingularValueDecomposition(u);

      // set singular values < thresholdLevel to zero
      // ChopSingularValues makes only sense if all columns of the x matrix have the same variance
      //decomposition.ChopSingularValues(1E-5);
      // recalculate the parameters with the chopped singular values
      _decomposition.Backsubstitution(scaledY,_parameter);

      _chiSquare = 0;
      for(int i=0;i<numberOfData;i++)
      {
        double ypredicted=0;
        for(int j=0;j<numberOfParameter;j++)
          ypredicted += _parameter[j]*xbase[i,j];
        double deviation = yarr[i]-ypredicted;
        _predicted[i] = ypredicted;
        _residual[i] = deviation;
        _chiSquare += deviation*deviation;
      }
    
      _covarianceMatrix = _decomposition.GetCovariances();


      //calculate the reduced prediction variance x'(X'X)^(-1)x
      for(int i=0;i<numberOfData;i++)
      {
        double total = 0;
        for(int j=0;j<numberOfParameter;j++)
        {
          double sum=0;
          for(int k=0;k<numberOfParameter;k++)
            sum += _covarianceMatrix[j][k]*u[i,k];

          total += u[i,j]*sum;
        }
        _reducedPredictionVariance[i] = total;
      }
 
      return this;
    }

 /// <summary>
 /// Fits a data set linear to a given x base.
 /// </summary>
 /// <param name="xbase">The matrix of x values of the data set. Dimensions: numberOfData x numberOfParameters. The matrix is changed during calculation!</param>
 /// <param name="yarr">The array of y values of the data set.</param>
 /// <param name="stddev">The array of y standard deviations of the data set.</param>
 /// <param name="numberOfData">The number of data points (may be smaller than the array sizes of the data arrays).</param>
 /// <param name="numberOfParameter">The number of parameters to fit == size of the function base.</param>
 /// <param name="threshold">A treshold value (usually 1E-5) used to chop the unimportant singular values away.</param>
 public LinearFitBySvd(
   IROMatrix xbase, // NumberOfData, NumberOfParameters 
   double[] yarr,
   double[] stddev,
   int numberOfData,
   int numberOfParameter,
   double threshold)
 {
   Calculate(xbase,yarr,stddev,numberOfData,numberOfParameter,threshold);
 }

 /// <summary>
 /// Creates an analyis from preprocessed spectra and preprocessed concentrations.
 /// </summary>
 /// <param name="matrixX">The spectral matrix (each spectrum is a row in the matrix). They must at least be centered.</param>
 /// <param name="matrixY">The matrix of concentrations (each experiment is a row in the matrix). They must at least be centered.</param>
 /// <param name="maxFactors">Maximum number of factors for analysis.</param>
 /// <returns>A regression object, which holds all the loads and weights neccessary for further calculations.</returns>
 public static PLS2Regression CreateFromPreprocessed(IROMatrix matrixX, IROMatrix matrixY, int maxFactors)
 {
   PLS2Regression result = new PLS2Regression();
   result.AnalyzeFromPreprocessed(matrixX,matrixY,maxFactors);
   return result;
 }

    /// <summary>
    /// Creates an analyis from preprocessed spectra and preprocessed concentrations.
    /// </summary>
    /// <param name="matrixX">The spectral matrix (each spectrum is a row in the matrix). They must at least be centered.</param>
    /// <param name="matrixY">The matrix of concentrations (each experiment is a row in the matrix). They must at least be centered.</param>
    /// <param name="maxFactors">Maximum number of factors for analysis.</param>
    /// <returns>A regression object, which holds all the loads and weights neccessary for further calculations.</returns>
    protected override void AnalyzeFromPreprocessedWithoutReset(IROMatrix matrixX, IROMatrix matrixY, int maxFactors)
    {
      int numberOfFactors = _calib.NumberOfFactors = Math.Min(matrixX.Columns, maxFactors);

      MatrixMath.BEMatrix _xLoads   = new MatrixMath.BEMatrix(0,0);
      MatrixMath.BEMatrix _yLoads   = new MatrixMath.BEMatrix(0,0);
      MatrixMath.BEMatrix _W       = new MatrixMath.BEMatrix(0,0);
      MatrixMath.REMatrix _V       = new MatrixMath.REMatrix(0,0);
      _PRESS   = VectorMath.CreateExtensibleVector(0);

      ExecuteAnalysis(matrixX, matrixY, ref numberOfFactors, _xLoads, _yLoads, _W, _V, _PRESS);
      _calib.NumberOfFactors = Math.Min(_calib.NumberOfFactors,numberOfFactors);
      _calib.XLoads = _xLoads;
      _calib.YLoads = _yLoads;
      _calib.XWeights = _W;
      _calib.CrossProduct = _V;
    }

    } // end partial-least-squares-predict


    /// <summary>
    /// Get the prediction score matrix. To get the predictions, you have to multiply
    /// the spectras to predict by this prediction score matrix (in case of a single y-variable), this is
    /// simply the dot product.
    /// </summary>
    /// <param name="xLoads">Matrix of spectral loads [factors,spectral bins].</param>
    /// <param name="yLoads">Matrix of concentration loads[factors, number of concentrations].</param>
    /// <param name="W">Matrix of spectral weightings [factors,spectral bins].</param>
    /// <param name="V">Cross product matrix[1,factors].</param>
    /// <param name="numFactors">Number of factors to use to calculate the score matrix.</param>
    /// <param name="predictionScores">Output: the resulting score matrix[ spectral bins, numberOfConcentrations]</param>
    public static void GetPredictionScoreMatrix(
      IROMatrix xLoads, // x-loads matrix
      IROMatrix yLoads, // y-loads matrix
      IROMatrix W, // weighting matrix
      IROMatrix V,  // Cross product vector
      int numFactors, // number of factors to use for prediction
      IMatrix predictionScores 
      )
    {

      Matrix bidiag = new Matrix(numFactors,numFactors);
      IROMatrix subweights = MatrixMath.ToROSubMatrix(W,0,0,numFactors,W.Columns);
      IROMatrix subxloads  = MatrixMath.ToROSubMatrix(xLoads,0,0,numFactors,xLoads.Columns);
      MatrixMath.MultiplySecondTransposed(subxloads,subweights,bidiag);
      IMatrix invbidiag = bidiag.Inverse;

      
      Matrix subyloads = new Matrix(numFactors,yLoads.Columns);
      MatrixMath.Submatrix(yLoads,subyloads,0,0);
      // multiply each row of the subyloads matrix by the appropriate weight
      for(int i=0;i<subyloads.Rows;i++)
        for(int j=0;j<subyloads.Columns;j++)
          subyloads[i,j] *= V[0,i];

      Matrix helper = new Matrix(numFactors,yLoads.Columns);
      MatrixMath.Multiply(invbidiag,subyloads,helper);

      //Matrix scores = new Matrix(yLoads.Columns,xLoads.Columns);
      //MatrixMath.MultiplyFirstTransposed(helper,subweights,predictionScores); // we calculate the transpose of scores (i.e. scores are horizontal oriented here)
      MatrixMath.MultiplyFirstTransposed(subweights,helper,predictionScores); // we calculate the transpose of scores (i.e. scores are horizontal oriented here)

      // now calculate the ys from the scores and the spectra

      //MultiplySecondTransposed(XU,scores,predictedY);
    } // end partial-least-squares-predict

    } // end partial-least-squares-predict


   
    public static void CalculateXLeverageFromPreprocessed(
      IROMatrix matrixX,
      IROMatrix W, // weighting matrix
      int numFactors, // number of factors to use for prediction
      IMatrix leverage, // Matrix of predicted y-values, must be same number of rows as spectra
      int leverageColumn
      )
    {

      // get the score matrix
      MatrixMath.BEMatrix weights = new MatrixMath.BEMatrix(numFactors,W.Columns);
      MatrixMath.Submatrix(W,weights,0,0);
      MatrixMath.BEMatrix scoresMatrix = new MatrixMath.BEMatrix(matrixX.Rows,weights.Rows);
      MatrixMath.MultiplySecondTransposed(matrixX,weights,scoresMatrix);
    
      MatrixMath.SingularValueDecomposition decomposition = MatrixMath.GetSingularValueDecomposition(scoresMatrix);


      for(int i=0;i<matrixX.Rows;i++)
        leverage[i,leverageColumn] = decomposition.HatDiagonal[i];
    }