C# Matrix.SetRow方法代码示例

        public override void ComputeRectificationMatrices()
        {
            //function[T1, T2, Pn1, Pn2] = rectify(Po1, Po2)
            //% RECTIFY: compute rectification matrices
            //% factorize old PPMs
            //[A1, R1, t1] = art(Po1);
            //[A2, R2, t2] = art(Po2);
            //% optical centers(unchanged)
            //c1 = - inv(Po1(:,1:3))*Po1(:,4);
            //        c2 = - inv(Po2(:,1:3))*Po2(:,4);
            //% new x axis(= direction of the baseline)
            //v1 = (c1-c2);
            //% new y axes(orthogonal to new x and old z)
            //v2 = cross(R1(3,:)',v1);
            //% new z axes(orthogonal to baseline and y)
            //v3 = cross(v1, v2);
            //% new extrinsic parameters
            //R = [v1'/norm(v1)
            //v2'/norm(v2)
            //v3'/norm(v3)];
            //% translation is left unchanged
            //% new intrinsic parameters(arbitrary)
            //A = (A1 + A2)./2;
            //A(1,2)=0; % no skew
            //% new projection matrices
            //Pn1 = A*[R - R * c1];
            //Pn2 = A*[R - R * c2];
            //% rectifying image transformation
            //T1 = Pn1(1:3, 1:3) * inv(Po1(1:3, 1:3));
            //T2 = Pn2(1:3,1:3)* inv(Po2(1:3,1:3));
            Vector<double> c1 = CalibrationData.Data.TranslationLeft;
            Vector<double> c2 = CalibrationData.Data.TranslationRight;

            Vector<double> v1 = c1 - c2;
            Vector<double> v2 = CalibrationData.Data.RotationLeft.Row(2).Cross(v1);
            Vector<double> v3 = v1.Cross(v2);

            _R = new DenseMatrix(3, 3);
            _R.SetRow(0, v1.Normalize(2));
            _R.SetRow(1, v2.Normalize(2));
            _R.SetRow(2, v3.Normalize(2));

            Matrix<double> halfRevolve = new DenseMatrix(3, 3);
            RotationConverter.EulerToMatrix(new double[] { 0.0, 0.0, Math.PI }, halfRevolve);
            _R = halfRevolve * _R;

            _K = (CalibrationData.Data.CalibrationLeft + CalibrationData.Data.CalibrationRight).Multiply(0.5);
            _K[0, 1] = 0.0;

            RectificationLeft = (_K * _R) * ((CalibrationData.Data.CalibrationLeft * CalibrationData.Data.RotationLeft).Inverse());
            RectificationRight = (_K * _R) * ((CalibrationData.Data.CalibrationRight * CalibrationData.Data.RotationRight).Inverse());
            ComputeScalingMatrices(ImageWidth, ImageHeight);

            RectificationLeft = _Ht_L * RectificationLeft;
            RectificationLeft = _Ht_R * RectificationRight;
            RectificationLeft_Inverse = RectificationLeft.Inverse();
            RectificationRight_Inverse = RectificationRight.Inverse();
        }

        public void InsertRowTest()
        {
            Matrix<double> A = new Matrix<double>(new double[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } });
            A.SetRow(2, new Matrix<double>(new double[,] { { 1, 2, 3 } }));
            Matrix<double> B = new Matrix<double>(new double[,] { { 1, 2, 3 }, { 1, 2, 3 }, { 7, 8, 9 } });

            Assert.AreEqual(B, A);
        }

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

//.........这里部分代码省略.........
                int seek;
                Vector<float> signal;
                //seek to the right position in the audio signal
                if (online)
                    //step back a complete windowSize after moving forward 1 hopSize
                    //so that the current position is at the stop of the window
                    seek = (int)((frame + 1) * HopSize - windowSize);
                else
                    //step back half of the windowSize so that the frame represents the centre of the window
                    seek = (int)(frame * HopSize - windowSize / 2);
                //read in the right portion of the audio
                if (seek >= _wav.Samples)
                    //stop of file reached
                    break;
                else if (seek + windowSize > _wav.Samples)
                {
                    //stop behind the actual audio stop, append zeros accordingly
                    int zeroAmount = seek + windowSize - _wav.Samples;
                    //var zeros = Vector<float>.Build.Dense(zeroAmount, 0);

                    var t = PythonUtilities.Slice<float>(cArray, seek, cArray.Length).ToArray();

                    //t.AddRange(zeros.ToList());

                    signal = _allocator.GetFloatVector(t.Length + zeroAmount);
                    for (int i = 0; i < t.Length; i++)
                    {
                        signal[i] = t[i];
                    }
                    //signal.SetValues(t);
                    //signal = Vector<float>.Build.DenseOfEnumerable(t);
                }
                else if (seek < 0)
                {
                    //start before actual audio start, pad with zeros accordingly
                    int zeroAmount = -seek;
                    var zeros = Vector<float>.Build.Dense(zeroAmount, 0).ToList();

                    var t = PythonUtilities.Slice<float>(cArray, 0, seek + windowSize).ToArray();
                    zeros.AddRange(t);

                    signal = _allocator.GetFloatVector(t.Length + zeroAmount);
                    signal.SetValues(zeros.ToArray());
                    //signal = Vector<float>.Build.DenseOfEnumerable(zeros);
                }
                else
                {
                    //normal read operation
                    var slice = PythonUtilities.Slice<float>(cArray, seek, seek + windowSize).ToArray();
                    signal = _allocator.GetFloatVector(slice.Length);
                    signal.SetValues(slice);

                    //signal = Vector<float>.Build.DenseOfEnumerable(PythonUtilities.Slice<float>(cArray, seek, seek + windowSize));
                }
                //multiply the signal with the window function
                signal = signal.PointwiseMultiply(Window);
                //only shift and perform complex DFT if needed
                if (phase)
                {
                    //circular shift the signal (needed for correct phase)
                    signal = NumpyCompatibility.FFTShift(signal);
                }
                //perform DFT
                //sanity check
                Debug.Assert(result.Length == signal.Count);
                for (int i = 0; i < result.Length; i++)
                {
                    result[i] = signal[i];
                }
                MathNet.Numerics.IntegralTransforms.Fourier.BluesteinForward(result, MathNet.Numerics.IntegralTransforms.FourierOptions.NoScaling);
                _STFT.SetRow(frame, result.Select(r => new Complex32((float)r.Real, (float)r.Imaginary)).Take(_ffts).ToArray());
                //var _newSTFTRow = result.Select(r => new Complex32((float)r.Real, (float)r.Imaginary)).Take(_ffts).ToArray();
                //_STFT.SetRow(frame, _newSTFTRow);
                //next frame
                _allocator.ReturnFloatVectorStorage((MathNet.Numerics.LinearAlgebra.Storage.DenseVectorStorage<float>)signal.Storage);
            }
            //magnitude spectrogram

            Spec = _allocator.GetFloatMatrix(_STFT.RowCount, _STFT.ColumnCount);
            if (phase)
                Phase = _allocator.GetFloatMatrix(_STFT.RowCount, _STFT.ColumnCount);
            for (int i = 0; i < Spec.RowCount; i++)
            {
                for (int j = 0; j < Spec.ColumnCount; j++)
                {
                    Spec.At(i, j, _STFT.At(i, j).Magnitude);
                    if (phase)
                        Phase.At(i, j, _STFT.At(i, j).Phase);
                }
            }
            //Spec = _STFT.Map(c => (float)c.Magnitude);

            //phase
            //if (phase)
            //{
            //    var imag = _STFT.Map(c => (float)c.Imaginary);
            //    var real = _STFT.Map(c => (float)c.Real);
            //    Phase = real.Map2((r, i) => (float)Math.Atan2(i,r), imag);
            //}
        }

        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
        public virtual void ComputeJacobianForLine(Matrix<double> J, int l, int p0)
        {
            double A = LineCoeffs[l, 0];
            double B = LineCoeffs[l, 1];
            double C = LineCoeffs[l, 2];
            double Ex = _sumX[l];
            double Ey = _sumY[l];
            double Exy = _sumXY[l];
            double Ex2 = _sumX2[l];
            double Ey2 = _sumY2[l];
            Vector<double> dEx = _dEx[l];
            Vector<double> dEy = _dEy[l];
            Vector<double> dExy = _dExy[l];
            Vector<double> dEx2 = _dEx2[l];
            Vector<double> dEy2 = _dEy2[l];
            double N = _n[l];

            //   A = (-b - sqrt(D)) / 2a
            //   D = b^2 + 4a^2
            //   a = Exy - Ex*Ey/N
            //   b = E(y^2) - (Ey)^2/N - (E(x^2) - (Ex)^2/N)
            //double a = Exy - Ex * Ey / N;
            //double b = Ey2 - Ex2 + (Ex * Ex - Ey * Ey) / N;
            //double D = b * b + 4 * a * a;

            double a = Get_a(l);
            double b = Get_b(l);
            double D = Get_D(l);

            // d(a) = d(Exy) - 1/N * (Ey*d(Ex) + *Ex*d(Ey))
            Vector<double> diff_a = dExy - (1.0 / N) * (dEx * Ey + dEy * Ex);

            // d(b) = d(Ey2) - 2/N * Ey * d(Ey) - d(Ex2) + 2/N * Ex * d(Ex)
            Vector<double> diff_b = dEy2 - dEx2 + (2.0 / N) * (dEx * Ex - dEy * Ey);

            // d(D) = 2b*d(b) + 8a*d(a)
            Vector<double> diff_D = (2.0 * b) * diff_b + (8.0 * a) * diff_a;

            // d(A) = -a(d(b) + 1/2sqrt(D) * d(D)) + d(a)(b+sqrt(D)) / 2a^2
            Vector<double> diff_A = (0.5 / (a * a)) *
                ((b + Math.Sqrt(D)) * diff_a - a * (diff_b + (0.5 / Math.Sqrt(D)) * diff_D));

            // d(C)/d(P) = -(1/N)(Ex*d(A) + A*d(Ex) + d(Ey))
            Vector<double> diff_C = (-1.0 / N) * (Ex * diff_A + A * dEx + dEy);

            double A212 = 1.0 / ((A * A + 1.0) * (A * A + 1.0));
            var line = CorrectedPoints[l];
            for(int p = 0; p < line.Count; ++p)
            {
                var point = line[p];
                double dist = A * point.Pf.X + point.Pf.Y + C;
                Vector<double> diff_dist = diff_A * point.Pf.X + A * point.Diff_Xf + point.Diff_Yf + diff_C;
                // J[point, k] = d(e)/d(P[k]) =
                //  2 * (Ax+y+C)/(A^2+1)^2 * d((Ax+y+C)) - A*d(A)*(Ax+y+C)
                //  d((Ax+y+C)) = a(A)*x + A*d(x) + d(y) + d(C)
                Vector<double> diff_e = (2.0 * dist * A212) * (
                    diff_dist * (A * A + 1.0) - (A * dist) * diff_A);

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

                J.SetRow(p0 + p, diff_e);
            }
        }

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

        public void CalculateDirectionDerivatives()
        {
            curls = DenseMatrix.Create(centroids.RowCount, 3, 0.0);
            totalDirectionDerivatives = new List<Matrix<double>>();

            for (int point = 0; point < centroids.RowCount; point++)
            {
                Matrix<double> tempMatrix = DenseMatrix.Create(3, 3, 0);
                List<int> neighbors = Calculator.NearestNeighbors(centroids.Row(point), this);
                Matrix<double> divideBy = DenseMatrix.Create(3, 3, 0);
                foreach (int neighbor in neighbors)
                {
                    if (neighbor != point)
                    {
                        for (int i = 0; i < 3; i++)
                        {
                            for (int j = 0; j < 3; j++)
                            {
                                if (Math.Abs(centroids.Row(neighbor)[j] - centroids.Row(point)[j]) > 0.001)
                                {
                                    tempMatrix[i, j] += (directions.Row(neighbor)[i] - directions.Row(point)[i]) / (centroids.Row(neighbor)[j] - centroids.Row(point)[j]);
                                    divideBy[i, j]++;
                                }
                            }
                        }
                    }
                }
                tempMatrix.PointwiseDivide(divideBy);
                curls.SetRow(point, DenseVector.OfArray(new double[] { tempMatrix[2, 1] - tempMatrix[1, 2], tempMatrix[0, 2] - tempMatrix[2, 0], tempMatrix[1, 0] - tempMatrix[0, 1] }));
                totalDirectionDerivatives.Add(tempMatrix);
            }
        }