C# DenseMatrix.SubMatrix方法代码示例

//.........这里部分代码省略.........
                {
                    c.SetColumn(f, c.Column(f) * generator.NextDouble());
                }
            }

            // todo:
            // generator.shuffle(C)

            // Loop over all clusters
            int pos = 0;
            int posEnd = 0;

            for (int k = 0; k < nClusters; k++)
            {
                // Number of samples in cluster k
                int nSamplesK = nSamplesPerCluster[k];

                // Define the range of samples
                pos = posEnd;
                posEnd = pos + nSamplesK;

                // Assign labels
                for (int l = pos; l < posEnd; l++)
                {
                    y[l] = k % nClasses;
                }

                // Draw features at random
                var subMatrix = DenseMatrix.CreateRandom(
                    nSamplesK,
                    nInformative,
                    new Normal { RandomSource = generator });

                x.SetSubMatrix(pos, nSamplesK, 0, nInformative, subMatrix);

                // Multiply by a random matrix to create co-variance of the features
                var uniform = new ContinuousUniform(-1, 1) { RandomSource = generator };
                Matrix a = DenseMatrix.CreateRandom(nInformative, nInformative, uniform);

                x.SetSubMatrix(
                    pos,
                    nSamplesK,
                    0,
                    nInformative,
                    x.SubMatrix(pos, nSamplesK, 0, nInformative) * a);

                // Shift the cluster to a vertice
                var v = x.SubMatrix(pos, nSamplesK, 0, nInformative).AddRowVector(c.Row(k));
                x.SetSubMatrix(pos, nSamplesK, 0, nInformative, v);
            }

            // Create redundant features
            if (nRedundant > 0)
            {
                var uniform = new ContinuousUniform(-1, 1) { RandomSource = generator };
                Matrix b = DenseMatrix.CreateRandom(nInformative, nRedundant, uniform);
                x.SetSubMatrix(
                    0,
                    x.RowCount,
                    nInformative,
                    nRedundant,
                    x.SubMatrix(0, x.RowCount, 0, nInformative) * b);
            }

            // Repeat some features
            if (nRepeated > 0)

        /// <summary>
        /// Run example
        /// </summary>
        public void Run()
        {
            // Format vector output to console
            var formatProvider = (CultureInfo)CultureInfo.InvariantCulture.Clone();
            formatProvider.TextInfo.ListSeparator = " ";

            // Create new empty square matrix
            var matrix = new DenseMatrix(10);
            Console.WriteLine(@"Empty matrix");
            Console.WriteLine(matrix.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 1. Fill matrix by data using indexer []
            var k = 0;
            for (var i = 0; i < matrix.RowCount; i++)
            {
                for (var j = 0; j < matrix.ColumnCount; j++)
                {
                    matrix[i, j] = k++;
                }
            }

            Console.WriteLine(@"1. Fill matrix by data using indexer []");
            Console.WriteLine(matrix.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 2. Fill matrix by data using At. The element is set without range checking.
            for (var i = 0; i < matrix.RowCount; i++)
            {
                for (var j = 0; j < matrix.ColumnCount; j++)
                {
                    matrix.At(i, j, k--);
                }
            }

            Console.WriteLine(@"2. Fill matrix by data using At");
            Console.WriteLine(matrix.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 3. Clone matrix
            var clone = matrix.Clone();
            Console.WriteLine(@"3. Clone matrix");
            Console.WriteLine(clone.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 4. Clear matrix
            clone.Clear();
            Console.WriteLine(@"4. Clear matrix");
            Console.WriteLine(clone.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 5. Copy matrix into another matrix
            matrix.CopyTo(clone);
            Console.WriteLine(@"5. Copy matrix into another matrix");
            Console.WriteLine(clone.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 6. Get submatrix into another matrix
            var submatrix = matrix.SubMatrix(2, 2, 3, 3);
            Console.WriteLine(@"6. Copy submatrix into another matrix");
            Console.WriteLine(submatrix.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 7. Get part of the row as vector. In this example: get 4 elements from row 5 starting from column 3
            var row = matrix.Row(5, 3, 4);
            Console.WriteLine(@"7. Get part of the row as vector");
            Console.WriteLine(row.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 8. Get part of the column as vector. In this example: get 3 elements from column 2 starting from row 6
            var column = matrix.Column(2, 6, 3);
            Console.WriteLine(@"8. Get part of the column as vector");
            Console.WriteLine(column.ToString("#0.00\t", formatProvider));
            Console.WriteLine();

            // 9. Get columns using column enumerator. If you need all columns you may use ColumnEnumerator without parameters
            Console.WriteLine(@"9. Get columns using column enumerator");
            foreach (var keyValuePair in matrix.ColumnEnumerator(2, 4))
            {
                Console.WriteLine(@"Column {0}: {1}", keyValuePair.Item1, keyValuePair.Item2.ToString("#0.00\t", formatProvider));
            }

            Console.WriteLine();

            // 10. Get rows using row enumerator. If you need all rows you may use RowEnumerator without parameters
            Console.WriteLine(@"10. Get rows using row enumerator");
            foreach (var keyValuePair in matrix.RowEnumerator(4, 3))
            {
                Console.WriteLine(@"Row {0}: {1}", keyValuePair.Item1, keyValuePair.Item2.ToString("#0.00\t", formatProvider));
            }

            Console.WriteLine();

            // 11. Convert matrix into multidimensional array
            var data = matrix.ToArray();
            Console.WriteLine(@"11. Convert matrix into multidimensional array");
            for (var i = 0; i < data.GetLongLength(0); i++)
//.........这里部分代码省略.........

        public DenseMatrix createMatrix(Model model)
        {
            int numConstraints = model.Constraints.Count;
            int numDecisionVars = model.Goal.Coefficients.Length;
            int varCounter = numDecisionVars;
            //  matrix(rows, columns)
            DenseMatrix coefficients = new DenseMatrix(numConstraints, numDecisionVars);
            DenseMatrix artificialVars = new DenseMatrix(numConstraints, 1);
            var constraintCounter = 0;
            this.rhsValues = new DenseVector(numConstraints);
            this.basics = new List<int>();
            this.artificialRows = new List<int>();
            foreach (var constraint in model.Constraints) {
                rhsValues[constraintCounter] = constraint.Value;

                // if the constraint RHS is negative, invert the coefficients and flip the inequality sign
                if (constraint.Value < 0)
                {
                    for (int i = 0; i< model.Goal.Coefficients.Length; i++) {
                        model.Goal.Coefficients[i] = model.Goal.Coefficients[i] * -1;
                    }
                    if (constraint.Relationship == Relationship.LessThanOrEquals)
                    {
                        constraint.Relationship = Relationship.GreaterThanOrEquals;
                    }
                    else if (constraint.Relationship == Relationship.GreaterThanOrEquals)
                    {
                        constraint.Relationship = Relationship.LessThanOrEquals;
                    }
                    // also flip the rhs value which we already put in the array for the simplex setup
                    rhsValues[constraintCounter] = rhsValues[constraintCounter] * -1;
                }

                coefficients.SetRow(constraintCounter, 0, constraint.Coefficients.Length, new DenseVector(constraint.Coefficients));
                // if it's a less than, add a slack column to the coefs matrix
                if (constraint.Relationship == Relationship.LessThanOrEquals)
                {
                    DenseVector slack = DenseVector.Create(model.Constraints.Count, delegate(int s) { return 0; });
                    slack.At(constraintCounter, 1);
                    coefficients = (DenseMatrix)coefficients.Append(slack.ToColumnMatrix());

                    this.basics.Add(varCounter);
                }
                else
                {
                    // Need to add an artificial variable for >= and = constraints

                    DenseVector surplus = DenseVector.Create(model.Constraints.Count, delegate(int s) { return 0; });
                    surplus.At(constraintCounter, -1);
                    coefficients = (DenseMatrix)coefficients.Append(surplus.ToColumnMatrix());

                    DenseVector artificial = DenseVector.Create(model.Constraints.Count, delegate(int s) { return 0; });
                    artificial.At(constraintCounter, 1);
                    artificialVars = (DenseMatrix)artificialVars.Append(artificial.ToColumnMatrix());

                    // Keeps track of the rows with artificial variable, for setting w
                    artificialRows.Add(constraintCounter);
                }
                varCounter++;
                constraintCounter++;
            }

            // put the constraints and stuff into the matrix
            if (artificialVars.ColumnCount > 1)
            {
                artificialVars = (DenseMatrix)artificialVars.SubMatrix(0, artificialVars.RowCount, 1, artificialVars.ColumnCount - 1);

                for (int i = coefficients.ColumnCount; i < coefficients.ColumnCount + artificialVars.ColumnCount; i++)
                {
                    this.basics.Add(i);
                }

                coefficients = (DenseMatrix)coefficients.Append(artificialVars);

                numArtificial = artificialVars.ColumnCount;
            }
            else
            {
                numArtificial = 0;
            }

            return coefficients;
        }

        private void optimize(DenseMatrix coefficients, DenseVector objFunValues, bool artifical)
        {
            //for calculations on the optimal solution row
            int cCounter,
                width = coefficients.ColumnCount;
            DenseVector cBVect = new DenseVector(basics.Count);

            //Sets up the b matrix
            DenseMatrix b = new DenseMatrix(basics.Count, 1);

            //basics will have values greater than coefficients.ColumnCount - 1 if there are still artificial variables
            //or if Nathan is bad and didn't get rid of them correctly
            foreach (int index in basics)
            {
                b = (DenseMatrix)b.Append(DenseVector.OfVector(coefficients.Column(index)).ToColumnMatrix());
            }
            // removes the first column
            b = (DenseMatrix)b.SubMatrix(0, b.RowCount, 1, b.ColumnCount - 1);

            double[] cPrimes = new double[width];
            double[] rhsOverPPrime;
            DenseMatrix[] pPrimes = new DenseMatrix[width];
            DenseMatrix bInverse;

            int newEntering, exitingRow;

            bool optimal = false;

            if(artifical)
            {
                rhsOverPPrime = new double[numConstraints + 1];
            }
            else
            {
                rhsOverPPrime = new double[numConstraints];
            }

            while (!optimal)
            {
                //calculates the inverse of b for this iteration
                bInverse = (DenseMatrix)b.Inverse();

                //updates the C vector with the most recent basic variables
                cCounter = 0;
                foreach (int index in basics)
                {
                    cBVect[cCounter++] = objFunValues.At(index);
                }

                //calculates the pPrimes and cPrimes
                for (int i = 0; i < coefficients.ColumnCount; i++)
                {
                    if (!basics.Contains(i))
                    {
                        pPrimes[i] = (DenseMatrix)bInverse.Multiply((DenseMatrix)coefficients.Column(i).ToColumnMatrix());

                        //c' = objFunVals - cB * P'n
                        //At(0) to turn it into a double
                        cPrimes[i] = objFunValues.At(i) - (pPrimes[i].LeftMultiply(cBVect)).At(0);
                    }
                    else
                    {
                        pPrimes[i] = null;
                    }
                }

                //RHS'
                xPrime = (DenseMatrix)bInverse.Multiply((DenseMatrix)rhsValues.ToColumnMatrix());

                //Starts newEntering as the first nonbasic
                newEntering = -1;
                int iter = 0;
                while(newEntering == -1)
                {
                    if(!basics.Contains(iter))
                    {
                        newEntering = iter;
                    }

                    iter++;
                }

                //new entering becomes the small cPrime that corresponds to a non-basic value
                for (int i = 0; i < cPrimes.Length; i++)
                {
                    if (cPrimes[i] < cPrimes[newEntering] && !basics.Contains(i))
                    {
                        newEntering = i;
                    }
                }

                //if the smallest cPrime is >= 0, ie they are all positive
                if (cPrimes[newEntering] >= 0)
                {
                    optimal = true;
                }
                else
                {
                    //fix me to deal with if all these values are negative
                    exitingRow = 0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........

            //First column of U
            U = new DenseMatrix(Rjk.RowCount, 1);
            U.SetColumn(0, Rjk.Column(0));

            // Norm of (approximate) Z, to test error
            double d1 = U.L2Norm();
            double d2 = V.L2Norm();
            double normZ = d1 * d1 * d2 * d2;

            //Find 2nd row index I(2)
            int row = 0;
            max = -1.0;

            foreach (int ind in i)
            {
                if (Math.Abs(Rjk[ind, 0]) > max)
                {
                    max = Math.Abs(Rjk[ind, 0]);
                    row = ind;
                }
            }

            //I[1] = i[row];
            I[1] = row;
            i.Remove(I[1]);

            //Iteration
            for (int k = 1; k < Math.Min(M, N); k++)
            {
                //Update (Ik)th row of the approximate error matrix:
                List<int> t1 = new List<int>();
                t1.Add(m[I[k]]);
                Rik = (Matrix)(UserImpedance(t1, n) - U.SubMatrix(I[k], 1, 0, U.ColumnCount).Multiply(V));
                //CHECKED OK all code before works fine
                //Find kth column index Jk
                max = -1.0;
                col = 0;

                foreach (int ind in j)
                {
                    if (Math.Abs(Rik[0, ind]) > max)
                    {
                        max = Math.Abs(Rik[0, ind]);
                        col = ind;
                    }
                }

                J[k] = col;
                j.Remove(J[k]);

                //Terminate if R(I(k),J(k)) == 0
                if (Rik[0, J[k]] == 0)
                {
                    break;
                }

                //Set k-th row of V equal to normalized error
                Matrix Vk = (Matrix)Rik.Divide(Rik[0, J[k]]);

                //Update (Jk)th column of the approximate error matrix
                List<int> n1 = new List<int>();
                n1.Add(n[J[k]]);
                Rjk = (Matrix)(UserImpedance(m, n1) - U.Multiply(V.SubMatrix(0, V.RowCount, J[k], 1)));

                // Set k-th column of U equal to updated error

        static Tuple<int, double, double> RunValidation(double[][] trainingData, double[][] validationData, double[][] outSampleData, int k)
        {
            int N = trainingData.Length;
              var ZZ = new DenseMatrix(N, 8);
              var Y = new DenseVector(N);
              for (int j = 0; j < N; j++)
              {
            double x1 = trainingData[j][0], x2 = trainingData[j][1];
            ZZ[j, 0] = 1;
            ZZ[j, 1] = x1;
            ZZ[j, 2] = x2;
            ZZ[j, 3] = x1 * x1;
            ZZ[j, 4] = x2 * x2;
            ZZ[j, 5] = x1 * x2;
            ZZ[j, 6] = Math.Abs(x1 - x2);
            ZZ[j, 7] = Math.Abs(x1 + x2);

            Y[j] = trainingData[j][2];
              }

              var Z = ZZ.SubMatrix(0, N, 0, k + 1);
              var WW = Z.TransposeThisAndMultiply(Z).Inverse().TransposeAndMultiply(Z).Multiply(Y);

              var W = new DenseVector(8);
              WW.CopySubVectorTo(W, 0, 0, k + 1);

              Func<double, double, double> h = (x1, x2) =>
            W[0] + W[1] * x1 + W[2] * x2 + W[3] * x1 * x1 + W[4] * x2 * x2 + W[5] * x1 * x2
            + W[6] * Math.Abs(x1 - x2) + W[7] * Math.Abs(x1 + x2) >= 0 ? 1.0 : -1.0;

              return Tuple.Create(
            k,
            (validationData.Count(v => Math.Sign(h(v[0], v[1])) != Math.Sign(v[2])) + 0.0) / validationData.Length,
            (outSampleData.Count(v => Math.Sign(h(v[0], v[1])) != Math.Sign(v[2])) + 0.0) / outSampleData.Length);
        }

        public static ResultFEM Solve(Matrix<double> K, Matrix<double> f, int[] bc)
        {
            // Create dof array:
            int ndof = f.RowCount;

            // Initialize displacement vector
            Matrix u = new DenseMatrix(ndof, 1);

            int[] dof = new int[ndof];
            for (int i = 0; i < ndof; i++)
            {
                dof[i] = i;
            }
            // Create the free dofs:
            int[] fdof = dof.Except(bc).ToArray();

            int nfdof = fdof.Length; // Number of free dofs.
            int nbdof = ndof - nfdof; //Constrained dofs

            // Create the permutation array which puts the constrained dofs last:
            int[] permute = bc.Union(fdof).ToArray();

            Permutation perm = new Permutation(permute);
            Permutation invPerm = perm.Inverse();
            Matrix<double> KPermuted = DenseMatrix.OfMatrix(K);
            KPermuted.PermuteRows(invPerm);
            KPermuted.PermuteColumns(invPerm);

            // Split K:::
            Matrix<double> Kff = KPermuted.SubMatrix(nbdof, nfdof, nbdof, nfdof);
            System.Console.WriteLine(Kff);  //Down right corner of matrix
            Matrix<double> Kfp = KPermuted.SubMatrix(nbdof, nfdof, 0, nbdof); //Down left corner of matrix

            // Split F:::
            Matrix<double> fPermuted = DenseMatrix.OfMatrix(f);
            fPermuted.PermuteRows(invPerm);
            Matrix<double> ff = fPermuted.SubMatrix(nbdof, nfdof, 0, 1);

            Matrix<double> up = u.SubMatrix(0, nbdof, 0, 1); // Must set up to constrained values in bc.

            // Solve for the unknown displacements:::
            Matrix<double> s = Kff.Solve(ff.Subtract(Kfp.Multiply(up)));
            u.SetSubMatrix(nbdof, 0, s); // Set displacements back.

            System.Console.WriteLine(u);
            // Permute back u
            u.PermuteRows(perm);

            System.Console.WriteLine("U after permut");
            System.Console.WriteLine(K);
            System.Console.WriteLine(u);
            System.Console.WriteLine(f);

            // Get reaction forces:
            Matrix<double> Q = K.Multiply(u).Subtract(f);

            ResultFEM result = new ResultFEM(u, Q);

            return result;
        }