C# Matrix.SetColumn方法代码示例

        public Molecule(double[,] inputAtoms, bool getDirections)
        {
            numAtoms = inputAtoms.GetLength(0);
            atoms = DenseMatrix.OfArray(inputAtoms);
            centroid = atoms.ColumnSums() / (double)numAtoms;

            if (getDirections)
            {
                for (int i = 0; i < 3; i++)
                {
                    atoms.SetColumn(i, atoms.Column(i) - centroid[i]);
                }
                direction = atoms.Svd(true).VT.Row(0);
                if ((atoms.Row(atoms.RowCount - 1) - atoms.Row(0)).DotProduct(direction) < 0)
                {
                    direction *= -1;
                }
            }
        }

 /// <summary>
 /// Sets the rows and columns of the stiffness corresponding 
 /// to the given DOF to 0, and the diagonal element to 1.
 /// </summary>
 private void EraseDOF(ref Matrix<double> k, int index)
 {
     Vector<double> zeroes = Vector<double>.Build.Dense(k.RowCount, 0);
     k.SetRow(index, zeroes);
     k.SetColumn(index, zeroes);
     k[index, index] = 1;
 }

        // Computes d(ei)/d(P) for ith line numericaly
        public void ComputeJacobian_Numerical(Matrix<double> J)
        {
            Vector<double> error_n = new DenseVector(_currentErrorVector.Count);
            Vector<double> error_p = new DenseVector(_currentErrorVector.Count);
            for(int k = 0; k < ResultsVector.Count; ++k)
            {
                double oldK = ResultsVector[k];
                double k_n = Math.Abs(oldK) > float.Epsilon ? oldK * (1 - NumericalDerivativeStep) : -NumericalDerivativeStep * 0.01;
                double k_p = Math.Abs(oldK) > float.Epsilon ? oldK * (1 + NumericalDerivativeStep) : NumericalDerivativeStep * 0.01;

                ResultsVector[k] = k_n;
                UpdateAll();
                ComputeErrorVector(error_n);

                ResultsVector[k] = k_p;
                UpdateAll();
                ComputeErrorVector(error_p);

                Vector<double> diff_e = 1.0 / (k_p - k_n) * (error_p - error_n);
                J.SetColumn(k, diff_e);

                ResultsVector[k] = oldK;

                // Throw if NaN found in jacobian
                bool nanInfFound = diff_e.Exists((e) => { return double.IsNaN(e) || double.IsInfinity(e); });
                if(nanInfFound)
                {
                    throw new NotFiniteNumberException("NaN or Infinity found on jacobian");
                }
            }

            UpdateAll();
        }

        private static Matrix<double> ntrp45(Matrix<double> tinterp, double t, Matrix<double> y, double h, Matrix<double> f)
        {
	        Matrix<double> BI = new Matrix<double>(new double[,] {{1.0,	-183.0/64.0,     37.0/12.0,    -145.0/128.0},
					                              {0.0,  0.0,            0.0,           0.0},
                                                  {0.0,  1500.0/371.0,  -1000.0/159.0,  1000.0/371.0},
                                                  {0.0, -125.0/32.0,     125.0/12.0,   -375.0/64.0}, 
                                                  {0.0,  9477.0/3392.0, -729.0/106.0,   25515.0/6784.0},
                                                  {0.0, -11.0/7.0,       11.0/3.0,     -55.0/28.0},
                                                  {0.0,  3.0/2.0,       -4.0,           5.0/2.0}});


            Matrix<double> s = (tinterp - t)/h;
            Matrix<double> S = new Matrix<double>(4, s.NumCols);

            for (int i = 1; i <= 4; i++)
                S.SetRow(i, s);

            Matrix<double> Y = new Matrix<double>(y.NumRows, tinterp.NumCols);

            for (int i = 1; i <= Y.NumCols; i++)
                Y.SetColumn(i, y);

            Matrix<double> yinterp = new Matrix<double>(y.NumRows, tinterp.NumCols);

            yinterp = Y + f*(h*BI)*Matrix<double>.Cumprod(S);

            return yinterp;
        }

        // Computes d(ei)/d(P) for ith line numericaly
        public void ComputeJacobianForLine_Numerical(Matrix<double> J, int l, int p0)
        {
            Vector<double> error_n = new DenseVector(_currentErrorVector.Count);
            Vector<double> error_p = new DenseVector(_currentErrorVector.Count);
            for(int k = 0; k < DistortionModel.ParametersCount; ++k)
            {
                double oldK = DistortionModel.Parameters[k];
                double k_n = Math.Abs(oldK) > float.Epsilon ? oldK * (1 - NumericalDerivativeStep) : -NumericalDerivativeStep * 0.01;
                double k_p = Math.Abs(oldK) > float.Epsilon ? oldK * (1 + NumericalDerivativeStep) : NumericalDerivativeStep * 0.01;

                DistortionModel.Parameters[k] = k_n;
                UpdateAll(l);
                ComputeErrorVector(error_n);
                Vector<double> error_n_line = error_n.SubVector(p0, LinePoints[l].Count);

                DistortionModel.Parameters[k] = k_p;
                UpdateAll(l);
                ComputeErrorVector(error_p);
                Vector<double> error_p_line = error_p.SubVector(p0, LinePoints[l].Count);

                Vector<double> diff_e = (1.0 / (k_p - k_n)) * (error_p_line - error_n_line);
                J.SetColumn(k, p0, LinePoints[l].Count, diff_e);

                DistortionModel.Parameters[k] = oldK;

                // Throw if NaN found in jacobian
                bool nanInfFound = diff_e.Exists((e) => { return double.IsNaN(e) || double.IsInfinity(e); });
                if(nanInfFound)
                {
                    throw new DivideByZeroException("NaN or Infinity found on jacobian");
                }
            }

            UpdateAll(l);
        }

        public static void PCA(DenseMatrix input, out Vector<double> latent, out Matrix<double> score, out Matrix<double> coeff)
        {
            int n = input.RowCount;
            int p = input.ColumnCount;

            //de-mean input
            var tmpInput = DenseMatrix.OfMatrix(input);
            for (int i = 0; i < tmpInput.ColumnCount; i++)
            {
                double avg = tmpInput.Column(i).Average();
                tmpInput.SetColumn(i, tmpInput.Column(i).Subtract(avg));
            }

            var svd = tmpInput.Svd(true);
            var sigma = svd.S;
            var tmpCoeff = svd.VT.Transpose();

            score = DenseMatrix.Create(n, p, (_, __) => 0);
            var U = svd.U.SubMatrix(0, n, 0, p);
            for (int i = 0; i < U.RowCount; i++)
            {
                score.SetRow(i, U.Row(i).PointwiseMultiply(sigma));
            }

            sigma = sigma.Divide(Math.Sqrt(n - 1));
            latent = sigma.PointwiseMultiply(sigma);

            //give the largest absolute value in each column a positive sign
            var maxIndices = tmpCoeff.EnumerateColumns().Select(x => x.AbsoluteMaximumIndex());
            var colSigns = maxIndices.Select((x, j) => Math.Sign(tmpCoeff[x, j])).ToList();
            for (int j = 0; j < tmpCoeff.ColumnCount; j++)
            {
                tmpCoeff.SetColumn(j, tmpCoeff.Column(j) * colSigns[j]);
                score.SetColumn(j, score.Column(j) * colSigns[j]);
            }

            coeff = tmpCoeff;
        }

 private void InsertOnesColumn(ref Matrix<double> data)
 {
     data = data.InsertColumn(0, Vector.Build.Dense(data.RowCount));
     data.SetColumn(0, Enumerable.Range(0, data.RowCount).Select(_ => 1.0).ToArray());
 }