Merge pull request #15 from rickyrombo/moar_tests

Moar tests

RaikesSimplexService/Implementation/Extensions/DataModelEqualValues.cs

@@ -10,14 +10,14 @@ public static class DataModelEqualValues
     {
         public static bool EqualValues(this Goal self, Goal other)
         {
-            bool sameCoefficients = self.Coefficients.SequenceEqual(other.Coefficients);
+            bool sameCoefficients = self.Coefficients.EqualValues(other.Coefficients);
             bool sameConstantTerm = self.ConstantTerm == other.ConstantTerm;
             return sameCoefficients && sameConstantTerm;
         }
 
         public static bool EqualValues(this LinearConstraint self, LinearConstraint other)
         {
-            bool sameCoefficients = self.Coefficients.SequenceEqual(other.Coefficients);
+            bool sameCoefficients = self.Coefficients.EqualValues(other.Coefficients);
             bool sameRelationship = self.Relationship == other.Relationship;
             bool sameValue = self.Value == other.Value;
             return sameCoefficients && sameRelationship && sameValue;
@@ -40,10 +40,18 @@ public static bool EqualValues(this Model self, Model other)
 
         public static bool EqualValues(this Solution self, Solution other)
         {
-            bool sameDecisions = self.Decisions.SequenceEqual(other.Decisions);
+            bool sameOptimal = self.OptimalValue.NearlyEqual(other.OptimalValue);
+            bool sameDecisions = self.Decisions.EqualValues(other.Decisions);
             bool sameQuality = self.Quality == other.Quality;
-            bool sameAltSols = self.AlternateSolutionsExist == other.AlternateSolutionsExist;
-            return sameDecisions && sameQuality && sameAltSols;
+            bool sameAltSols = self.AlternateSolutionsExist == other.AlternateSolutionsExist
+            return sameDecisions && sameQuality && sameOptimal && sameAltSols;
+        }
+
+        public static bool EqualValues(this double[] self, double[] other)
+        {
+            var pairs = self.Zip(other, (a, b) => new { First = a, Second = b });
+            bool allPairsEqual = pairs.All(pair => pair.First.NearlyEqual(pair.Second));
+            return allPairsEqual;
         }
 
     }

RaikesSimplexService/Implementation/Extensions/Delta.cs

@@ -0,0 +1,22 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace RaikesSimplexService.Implementation.Extensions
+{
+    public static class Delta
+    {
+        private static readonly double TOLERANCE = 0.00000001;
+
+        public static bool NearlyZero(this double d)
+        {
+            return NearlyEqual(d, 0.0);
+        }
+
+        public static bool NearlyEqual(this double self, double other)
+        {
+            return self >= other - TOLERANCE && self <= other + TOLERANCE;
+        }
+    }
+}

RaikesSimplexService/Implementation/Solver.cs

@@ -8,30 +8,46 @@
 using RaikesSimplexService.DataModel;
 using MathNet.Numerics.LinearAlgebra;
 using MathNet.Numerics.LinearAlgebra.Double;
+using RaikesSimplexService.Implementation.Extensions;
 
 namespace RaikesSimplexService.Implementation
 {
     public class Solver : ISolver
     {
-        private static readonly double ZERO_TOLERANCE = 0.0000001;
 
         public Solution Solve(Model m)
         {
             var model = StandardModel.FromModel(m);
             return SolveStandardModel(model);
         }
+
+        public Tuple<int, double> CnPrime(StandardModel model, Matrix<double> cb, Tuple<int, Vector<double>> p)
+        {
+            double cn = model.ObjectiveRow.At(0, p.Item1);
+            double cnPrime = cn - cb.Multiply(p.Item2).At(0);
+            return new Tuple<int, double>(p.Item1, cnPrime);
+        }
+
+        public Tuple<int, double> GetMinCnPrime(IEnumerable<Tuple<int, Vector<double>>> primes, StandardModel model, Matrix<double> cb, List<int> basicColumnIndices)
+        {
+            var cnPrimes = primes.Where(p => !basicColumnIndices.Contains(p.Item1)).Select(
+                    p => CnPrime(model, cb, p)
+                );
+            Tuple<int, double> minCnPrime = cnPrimes.OrderBy(s => s.Item2).FirstOrDefault(); // Get min
+            return minCnPrime;
+        }
         public Solution SolveStandardModel(StandardModel model)
         {
             //Get initial basic columns
-            var basicColumns = model.LHS.EnumerateColumnsIndexed().Where(v => v.Item2.Count(s => s != 0) == 1 && v.Item2.Any(s => s == 1)).ToList();
+            var basicColumns = model.LHS.EnumerateColumnsIndexed().Where(v => v.Item2.Count(s => !s.NearlyZero()) == 1 && v.Item2.Any(s => s.NearlyEqual(1))).ToList();
             var sol = new Solution();
             while (true)
             {
-                var basicColumnIndices = basicColumns.Select(s => s.Item1).ToList();
+                List<int> basicColumnIndices = basicColumns.Select(s => s.Item1).ToList();
                 var nonbasicColumns = model.LHS.EnumerateColumnsIndexed().Where(v => !basicColumnIndices.Contains(v.Item1)).ToList();
                 var bInv = Matrix<double>.Build.DenseOfColumnVectors(basicColumns.Select(s => s.Item2)).Inverse();
                 //Get the P1' P2' etc from the nonbasic columns * inverse basic matrix
-                var primes = model.LHS.EnumerateColumnsIndexed().Select(
+                List<Tuple<int, Vector<double>>> primes = model.LHS.EnumerateColumnsIndexed().Select(
                     s => new Tuple<int, Vector<double>>(s.Item1, bInv.Multiply(s.Item2))
                 ).ToList();
                 var xb = bInv.Multiply(model.RHS);
@@ -41,13 +57,9 @@ public Solution SolveStandardModel(StandardModel model)
                     .Select(s => s.Item2)
                 );
                 //Calculate C1' C2' etc and select the minimum - that's our entering basic variable
-                var enteringCol = primes.Select(
-                    p => (model.ObjectiveRow.At(0, p.Item1) - cb.Multiply(p.Item2)).Select(
-                        s => new Tuple<int, double>(p.Item1, (double)s)
-                    ).First() //There's only ever one element because the width of cb is equal to the height of the primes
-                ).OrderBy(s => s.Item2).FirstOrDefault();
+                var minCnPrime = GetMinCnPrime(primes, model, cb, basicColumnIndices);
                 //If all the C1' C2' etc are positive, then we're done - we've optimized the solution
-                if (enteringCol.Item2 >= 0)
+                if (minCnPrime.Item2 >= 0 || minCnPrime.Item2.NearlyZero())
                 {
                     if (model.ArtificialVariables > 0)
                     {
@@ -62,7 +74,7 @@ public Solution SolveStandardModel(StandardModel model)
                         };
                         var zRowIndex = phase2.LHS.Column(0)
                             .EnumerateIndexed()
-                            .Where(pair => pair.Item2 == 1)
+                            .Where(pair => pair.Item2.NearlyEqual(1.0))
                             .Select(pair => pair.Item1)
                             .FirstOrDefault();
                         phase2.LHS = phase2.LHS.RemoveRow(zRowIndex);
@@ -74,7 +86,7 @@ public Solution SolveStandardModel(StandardModel model)
                         for (int i = 0; i < model.ArtificialVariables; i++)
                         {
                             var artificialCol = phase2.LHS.Column(phase2.LHS.ColumnCount - 1).Enumerate();
-                            if (artificialCol.Count(s => s != 0) == 1 && artificialCol.Any(s => s == 1))
+                            if (artificialCol.Count(s => !s.NearlyZero()) == 1 && artificialCol.Any(s => s.NearlyEqual(1.0)))
                             {
                                 sol.Quality = SolutionQuality.Infeasible;
                                 return sol;
@@ -89,22 +101,30 @@ public Solution SolveStandardModel(StandardModel model)
                         sol.Decisions = new double[model.DecisionVariables];
                         _mapDecisionVariables(sol.Decisions, basicColumnIndices, xb);
                         sol.OptimalValue = _calculateGoalValue(sol.Decisions, model.OriginalModel.Goal.Coefficients);
-                        sol.AlternateSolutionsExist = sol.Decisions.Any(s => NearlyZero(s));
+                        sol.AlternateSolutionsExist = sol.Decisions.Any(s => s.NearlyZero());
+                        Console.WriteLine(model);
                     }
                     break;
                 }
                 //else, get the divisor from the Pn' we selected
-                var divisor = primes.Where(s => s.Item1 == enteringCol.Item1).FirstOrDefault().Item2;
-                var ratios = xb.PointwiseDivide(divisor);
+                var enteringVar = minCnPrime;
+                var divisor = primes.Where(s => s.Item1 == enteringVar.Item1).FirstOrDefault().Item2;
+                Vector<double> ratios = xb.PointwiseDivide(divisor);
+                List<Tuple<int, double>> columnsWithRatios = new List<Tuple<int, double>>();
+                for (int i = 0; i < basicColumns.Count; i++)
+                {
+                    columnsWithRatios.Add(new Tuple<int, double>(basicColumnIndices[i], ratios[i]));
+                }
                 //Get the minimum ratio that's > 0 - that's our exiting basic variable
-                var exitingCol = ratios.EnumerateIndexed().Where(s => s.Item2 > 0).OrderBy(s => s.Item2).FirstOrDefault();
-                if (exitingCol == null)
+                var exitCol = ratios.EnumerateIndexed().Where(s => s.Item2 > 0 && !s.Item2.NearlyZero() && model.ArtificialVariables == 0 || IndexArtificial(basicColumns[s.Item1].Item1, model)).OrderBy(s => s.Item2).FirstOrDefault();
+                if (exitCol == null)
                 {
                     sol.Quality = SolutionQuality.Unbounded;
                     break;
                 }
-                var newCol = nonbasicColumns.FirstOrDefault(s => s.Item1 == enteringCol.Item1);
-                basicColumns.RemoveAt(exitingCol.Item1);
+                var newCol = nonbasicColumns.FirstOrDefault(s => s.Item1 == enteringVar.Item1);
+                //basicColumns[exitCol.Item1] = newCol;
+                basicColumns.RemoveAt(exitCol.Item1);
                 int insertHere = 0;
                 foreach (var col in basicColumnIndices)
                 {
@@ -115,10 +135,16 @@ public Solution SolveStandardModel(StandardModel model)
                     insertHere++;
                 }
                 basicColumns.Insert(insertHere, newCol);
+
             }
             return sol;
         }
 
+        public bool IndexArtificial(int i, StandardModel model)
+        {
+            return i >= model.DecisionVariables + model.SlackVariables;
+        }
+
         /// <summary>
         /// Runs the goal equation on the optimized decision variables
         /// </summary>
@@ -152,10 +178,5 @@ private static void _mapDecisionVariables(double[] decisionVariables, List<int>
                 }
             }
         }
-
-        public static bool NearlyZero(double d)
-        {
-            return d >= -ZERO_TOLERANCE && d <= ZERO_TOLERANCE;
-        }
     }
 }

RaikesSimplexService/Implementation/StandardModel.cs

@@ -25,6 +25,7 @@ public enum OutputFormat
         public int DecisionVariables { get; private set; }
         public int SlackVariables { get; private set; }
         public int ArtificialVariables { get; private set; }
+        public double ConstantTerm { get; set; }
         #endregion
 
         public StandardModel(int constraintCount, int decisionCount, int slackCount, int artificialCount, Model originalModel)
@@ -37,6 +38,7 @@ public StandardModel(int constraintCount, int decisionCount, int slackCount, int
             this.SlackVariables = slackCount;
             this.ArtificialVariables = artificialCount;
             this.OriginalModel = originalModel;
+            this.ConstantTerm = originalModel.Goal.ConstantTerm;
         }
 
         /// <summary>
@@ -88,7 +90,7 @@ public static StandardModel FromModel(Model model)
                 standardModel.LHS.SetRow(i++, coeffs.ToArray<double>());
             }
             //Catch for artificial column (z)
-            var goalCoeffsList = (artificialCount == 0) ? new List<double>() : new List<double>{0};
+            var goalCoeffsList = (artificialCount == 0) ? new List<double>() : new List<double> { 0 };
             goalCoeffsList.AddRange(model.Goal.Coefficients);
             var goalCoeffs = goalCoeffsList.ToArray<double>();
             Array.Resize<double>(ref goalCoeffs, goalCoeffs.Length + slackCount + artificialCount);
@@ -121,9 +123,10 @@ private void AddWRow()
             this.LHS.At(this.LHS.RowCount - 1, 0, 1);
             this.ObjectiveRow.SetRow(0, new double[this.LHS.ColumnCount]);
             //Sum all variables for the w row
-            foreach(var row in this.LHS.EnumerateRows().Where(r => ContainsArtificalVariable(r)))
+            foreach (var rowAndIndex in this.LHS.EnumerateRowsIndexed().Where(rowAndIndex => ContainsArtificalVariable(rowAndIndex.Item2)))
             {
-                foreach (var pair in row.EnumerateIndexed())
+                this.ConstantTerm = this.ConstantTerm - this.RHS.ElementAt(rowAndIndex.Item1);
+                foreach (var pair in rowAndIndex.Item2.EnumerateIndexed())
                 {
                     this.ObjectiveRow.At(0, pair.Item1, this.ObjectiveRow.At(0, pair.Item1) - pair.Item2);
                 }
@@ -163,11 +166,12 @@ public string ToString(OutputFormat format)
         {
             if (format == OutputFormat.Expression)
             {
-                return string.Format("Constraints:\n{0}\nObjective: Maximize\nZ\t+ {1} \t= 0",
+                return string.Format("Constraints:\n{0}\nObjective: Maximize\nZ\t+ {1} \t= {2}",
                     string.Join("\n", LHS.EnumerateRowsIndexed().Select(
                         r => string.Format("{0}\t= {1}", _stringExpression(r.Item2), this.RHS.At(r.Item1))
                     )),
-                    _stringExpression(ObjectiveRow.Row(0))
+                    _stringExpression(ObjectiveRow.Row(0)),
+                    ConstantTerm
                 );
             }
             if (format == OutputFormat.Original)

RaikesSimplexService/RaikesSimplexService.csproj

@@ -58,6 +58,7 @@
     <Compile Include="DataModel\Solution.cs" />
     <Compile Include="DataModel\SolutionQuality.cs" />
     <Compile Include="Implementation\Extensions\DataModelEqualValues.cs" />
+    <Compile Include="Implementation\Extensions\Delta.cs" />
     <Compile Include="Implementation\StandardModel.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
     <Compile Include="Implementation\Solver.cs" />

UnitTests/Helpers/ModelGenerator.cs

@@ -114,11 +114,51 @@ public static Model GetTwoPhaseModel()
                 },
                 Goal = new Goal
                 {
-                    Coefficients = new double[] {6, 3},
+                    Coefficients = new double[] { 6, 3 },
                     ConstantTerm = 0
                 },
                 GoalKind = GoalKind.Maximize
             };
         }
+
+        public static Model GetAshuModel()
+        {
+            return new Model
+            {
+                Constraints = new List<LinearConstraint>
+                {
+                    new LinearConstraint 
+                    {
+                        Coefficients = new double[2] { 8, 12},
+                        Relationship = Relationship.GreaterThanOrEquals,
+                        Value = 24
+                    },
+                    new LinearConstraint 
+                    {
+                        Coefficients = new double[2] {12, 12},
+                        Relationship = Relationship.GreaterThanOrEquals,
+                        Value = 36
+                    },
+                    new LinearConstraint 
+                    {
+                        Coefficients = new double[2] {2, 1},
+                        Relationship = Relationship.GreaterThanOrEquals,
+                        Value = 4
+                    },
+                    new LinearConstraint 
+                    {
+                        Coefficients = new double[2] {1, 1},
+                        Relationship = Relationship.LessThanOrEquals,
+                        Value = 5
+                    },
+                },
+                Goal = new Goal
+                {
+                    Coefficients = new double[2] { .2, .3 },
+                    ConstantTerm = 0
+                },
+                GoalKind = GoalKind.Minimize
+            };
+        }
     }
 }

UnitTests/Helpers/SolutionGenerator.cs

@@ -15,7 +15,8 @@ public static Solution GetSimpleSolution()
             {
                 AlternateSolutionsExist = true,
                 OptimalValue = 20,
-                Decisions = new double[] { 0, 3.333333333333333, 0 }
+                Decisions = new double[] { 0, 3.333333333333333, 0 },
+                Quality = SolutionQuality.Optimal
             };
         }
 
@@ -26,6 +27,18 @@ public static Solution GetTwoPhaseSolution()
                 AlternateSolutionsExist = true,
                 OptimalValue = 6,
                 Decisions = new double[] { 1, 0 },
+                Quality = SolutionQuality.Optimal
+            };
+        }
+
+        public static Solution GetAshuSolution()
+        {
+            return new Solution
+            {
+                AlternateSolutionsExist = true,
+                OptimalValue = .6,
+                Quality = SolutionQuality.Optimal,
+                Decisions = new double[] { 3, 0 }
             };
         }
     }

UnitTests/Helpers/StandardModelGenerator.cs

@@ -28,5 +28,10 @@ public static StandardModel GetTwoPhaseStandardModel()
         {
             return StandardModel.FromModel(ModelGenerator.GetTwoPhaseModel());
         }
+
+        public static StandardModel GetAshuStandardModel()
+        {
+            return StandardModel.FromModel(ModelGenerator.GetAshuModel());
+        }
     }
 }