Finished

.gitignore

@@ -0,0 +1,4 @@
+.venv/
+.vscode/
+bin/
+obj/

Calculator.cs

@@ -1,80 +1,185 @@
 namespace omp
 {
-    public class OverlappingPoints : Exception { }
+    using CalculationResults = Dictionary<string, double[,]>;
 
-    public static class Calculator
+    public class CalculationException : Exception;
+
+    public class OverlappingPoints : CalculationException;
+    public class SingularMatrix : CalculationException;
+
+    public interface ICalculator
     {
-        public const double StationPDOP = -1000;
+        CalculationResults Calculate();
+    }
 
-        public static double[,] Jacobian1T(List<Point> stations, Point position)
+    public class Calculator2D(List<Point2D> stations, int xSize, int ySize) : ICalculator
+    {
+        private readonly List<Point2D> stations = stations;
+        private readonly int xSize = xSize, ySize = ySize;
+
+        public CalculationResults Calculate()
+        {
+            double[,] xdop = new double[xSize, ySize];
+            double[,] ydop = new double[xSize, ySize];
+            double[,] hdop = new double[xSize, ySize];
+
+            Parallel.For(0, xSize, i =>
+            {
+                for (int j = 0; j < ySize; ++j)
+                {
+                    try
+                    {
+                        double[,] JR = Jacobian(stations, new Point2D(i, j));
+
+                        double[,] mult = JR.MultTrans();
+
+                        double det = mult.Det();
+
+                        if (det == 0)
+                        {
+                            throw new SingularMatrix();
+                        }
+
+                        double sx2 = mult.Adj(0, 0) / det;
+                        double sy2 = mult.Adj(1, 1) / det;
+
+                        xdop[i, j] = Math.Sqrt(sx2);
+                        ydop[i, j] = Math.Sqrt(sy2);
+                        hdop[i, j] = Math.Sqrt(sx2 + sy2);
+                    }
+                    catch (CalculationException)
+                    {
+                        xdop[i, j] = double.NaN;
+                        ydop[i, j] = double.NaN;
+                        hdop[i, j] = double.NaN;
+                    }
+                }
+            });
+
+            return new CalculationResults
+            {
+                ["xdop"] = xdop,
+                ["ydop"] = ydop,
+                ["hdop"] = hdop,
+            };
+        }
+
+        public static double[,] Jacobian(List<Point2D> stations, Point2D position)
         {
             int n = stations.Count;
-            double[,] J = new double[4, n];
+            double[,] J = new double[n, 2];
 
             for (int i = 0; i < n; ++i)
             {
-                double distance = position.Distance(stations[i]);
+                double distance = position.DistanceTo(stations[i]);
 
                 if (distance == 0)
                 {
                     throw new OverlappingPoints();
                 }
 
-                J[0, i] = -(stations[i].x - position.x) / distance;
+                J[i, 0] = -(stations[i].x - position.x) / distance;
-                J[1, i] = -(stations[i].y - position.y) / distance;
+                J[i, 1] = -(stations[i].y - position.y) / distance;
-                J[2, i] = -(stations[i].z - position.z) / distance;
-                J[3, i] = 1; // Clock bias
             }
 
             return J;
         }
+    }
 
-        public static double PDOP(List<Point> stations, Point position)
+    public class Calculator3D(List<Point3D> stations, int xSize, int ySize, int zCoordinate) : ICalculator
+    {
+        private readonly List<Point3D> stations = stations;
+        private readonly int xSize = xSize, ySize = ySize;
+        private readonly int z = zCoordinate;
+
+        public CalculationResults Calculate()
         {
-            double[,] JR = Jacobian1T(stations, position);
+            double[,] xdop = new double[xSize, ySize];
+            double[,] ydop = new double[xSize, ySize];
+            double[,] hdop = new double[xSize, ySize];
+            double[,] vdop = new double[xSize, ySize];
+            double[,] pdop = new double[xSize, ySize];
 
-            double[,] mult = MatrixOps.MultTrans(JR);
+            Parallel.For(0, xSize, i =>
-
-            double det = MatrixOps.Det(mult);
-
-            double sx2 = MatrixOps.Adj(mult, 0, 0) / det;
-            double sy2 = MatrixOps.Adj(mult, 1, 1) / det;
-            double sz2 = MatrixOps.Adj(mult, 2, 2) / det;
-
-            return Math.Sqrt(sx2 + sy2 + sz2);
-        }
-
-        public static void PDOPs(double[,] pdops, Configurator conf, out double min, out double max)
-        {
-            min = double.MaxValue;
-            max = double.MinValue;
-
-            for (int i = 0; i < pdops.GetLength(0); ++i)
             {
-                for (int j = 0; j < pdops.GetLength(1); ++j)
+                for (int j = 0; j < ySize; ++j)
                 {
                     try
                     {
-                        pdops[i, j] = PDOP(conf.stations, new Point(i, j, conf.zCoordinate));
+                        double[,] JR = Jacobian(stations, new Point3D(i, j, z));
 
-                        if (pdops[i, j] > max)
+                        double[,] mult = JR.MultTrans();
-                        {
-                            max = pdops[i, j];
-                        }
 
-                        if (pdops[i, j] < min)
+                        double det = mult.Det();
-                        {
+
-                            min = pdops[i, j];
+                        double sx2 = mult.Adj(0, 0) / det;
-                        }
+                        double sy2 = mult.Adj(1, 1) / det;
+                        double sz2 = mult.Adj(2, 2) / det;
+
+                        xdop[i, j] = Math.Sqrt(sx2);
+                        ydop[i, j] = Math.Sqrt(sy2);
+                        hdop[i, j] = Math.Sqrt(sx2 + sy2);
+                        vdop[i, j] = Math.Sqrt(sz2);
+                        pdop[i, j] = Math.Sqrt(sx2 + sy2 + sz2);
                     }
                     catch (OverlappingPoints)
                     {
-                        pdops[i, j] = StationPDOP;
+                        xdop[i, j] = double.NaN;
+                        ydop[i, j] = double.NaN;
+                        hdop[i, j] = double.NaN;
+                        vdop[i, j] = double.NaN;
+                        pdop[i, j] = double.NaN;
                     }
                 }
+            });
+
+            return new CalculationResults
+            {
+                ["xdop"] = xdop,
+                ["ydop"] = ydop,
+                ["hdop"] = hdop,
+                ["vdop"] = vdop,
+                ["pdop"] = pdop,
+            };
+        }
+
+        public static double[,] Jacobian(List<Point3D> stations, Point3D position)
+        {
+            int n = stations.Count;
+            double[,] J = new double[n, 3];
+
+            for (int i = 0; i < n; ++i)
+            {
+                double distance = position.DistanceTo(stations[i]);
+
+                if (distance == 0)
+                {
+                    throw new OverlappingPoints();
+                }
+
+                J[i, 0] = -(stations[i].x - position.x) / distance;
+                J[i, 1] = -(stations[i].y - position.y) / distance;
+                J[i, 2] = -(stations[i].z - position.z) / distance;
             }
 
-            Console.WriteLine("{0}...{1}", min, max);
+            return J;
         }
     }
-}
+
+    public static class CalculatorFactory
+    {
+        public static ICalculator Create(IConfigurator configurator)
+        {
+            var stations = configurator.Stations;
+
+            return configurator.PointType switch
+            {
+                _ when configurator.PointType == typeof(Point2D) =>
+                    new Calculator2D(stations.Cast<Point2D>().ToList(), configurator.XSize, configurator.YSize),
+                _ when configurator.PointType == typeof(Point3D) =>
+                    new Calculator3D(stations.Cast<Point3D>().ToList(), configurator.XSize, configurator.YSize, configurator.ZCoordinate),
+                _ => throw new InvalidOperationException("Unsupported point type")
+            };
+        }
+    }
+}

Configurator.cs

@@ -1,13 +1,47 @@
 namespace omp
 {
-    public class Configurator
+    class UndergroundException : Exception;
+
+    public interface IConfigurator
     {
-        public int xsize;
+        Type PointType { get; }
-        public int ysize;
+
-        public int zCoordinate;
+        public int XSize { get; }
-        public List<Point> stations = [];
+        public int YSize { get; }
+        public int ZCoordinate { get; }
+        public List<Point> Stations { get; }
+    }
+
+    public class Configurator : IConfigurator
+    {
+        public Type PointType { get; }
+
+        public int XSize { get; }
+        public int YSize { get; }
+        public int ZCoordinate { get; }
+        public List<Point> Stations { get; } = [];
 
         public Configurator(string[] args)
+        {
+            List<string> lines = ReadAllLines(args);
+
+            bool hasZ = CheckHasZ(lines[1]);
+
+            if (hasZ)
+            {
+                PointType = typeof(Point3D);
+            }
+            else
+            {
+                PointType = typeof(Point2D);
+            }
+
+            (XSize, YSize) = ParseSize(lines[0]);
+            ZCoordinate = ParseZCoord(lines[1], hasZ);
+            ParseStations(lines, hasZ);
+        }
+
+        private static List<string> ReadAllLines(string[] args)
         {
             List<string> lines;
 
@@ -32,16 +66,7 @@ namespace omp
                 throw new ArgumentException("You must either provide path to parameters file as an only cli argument, or insert it into STDIN");
             }
 
-            ParseParameters(lines);
+            return lines;
-        }
-
-        private void ParseParameters(IReadOnlyList<string> lines)
-        {
-            bool hasZ = CheckHasZ(lines[1]);
-
-            ParseSize(lines[0]);
-            ParseZCoord(lines[1], hasZ);
-            ParseStations(lines, hasZ);
         }
 
         private static bool CheckHasZ(string line)
@@ -49,48 +74,65 @@ namespace omp
             return line.Split(' ').Length == 1; // line contains a single number, z coordinate
         }
 
-        private void ParseSize(string line)
+        private static (int, int) ParseSize(string line)
         {
             string[] size_parts = line.Split(' ', 2);
-            xsize = int.Parse(size_parts[0]);
+            int xSize = int.Parse(size_parts[0]);
-            ysize = int.Parse(size_parts[1]);
+            int ySize = int.Parse(size_parts[1]);
+
+            return (xSize, ySize);
         }
 
-        private void ParseZCoord(string line, bool hasZ)
+        private int ParseZCoord(string line, bool hasZ)
         {
+            int zCoordinate;
             if (hasZ)
             {
                 zCoordinate = int.Parse(line);
+                if (ZCoordinate < 0)
+                {
+                    throw new UndergroundException();
+                }
             }
             else
             {
-                zCoordinate = 0;
+                zCoordinate = -1;
             }
+
+            return zCoordinate;
         }
 
         private void ParseStations(IReadOnlyList<string> lines, bool hasZ)
         {
-            int firstIndex = 1 + (hasZ ? 1 : 0);
+            if (hasZ) Parse3DStations(lines);
-            for (int i = firstIndex; i < lines.Count; ++i)
+            else Parse2DStations(lines);
+        }
+
+        private void Parse2DStations(IReadOnlyList<string> lines)
+        {
+            for (int i = 1; i < lines.Count; ++i)
             {
-                string[] point_parts = lines[i].Split(' ', 2 + (hasZ ? 1 : 0));
+                string[] point_parts = lines[i].Split(' ', 2);
 
-                int x, y, z;
+                int x = int.Parse(point_parts[0]);
+                int y = int.Parse(point_parts[1]);
 
-                x = int.Parse(point_parts[0]) - 1;
+                Stations.Add(new Point2D(x, y));
-                y = int.Parse(point_parts[1]) - 1;
+            }
+        }
 
-                if (hasZ)
+        private void Parse3DStations(IReadOnlyList<string> lines)
-                {
+        {
-                    z = int.Parse(point_parts[2]);
+            for (int i = 2; i < lines.Count; ++i)
-                }
+            {
-                else
+                string[] point_parts = lines[i].Split(' ', 3);
-                {
-                    z = 0;
-                }
 
-                stations.Add(new Point(x, y, z));
+                int x = int.Parse(point_parts[0]);
+                int y = int.Parse(point_parts[1]);
+                int z = int.Parse(point_parts[2]);
+
+                Stations.Add(new Point3D(x, y, z));
             }
         }
     }
-}
+}

Drawer.cs

@@ -2,25 +2,39 @@ using SkiaSharp;
 
 namespace omp
 {
-    class Drawer
+    using CalculationResults = Dictionary<string, double[,]>;
-    {
-        public static void Draw(Configurator conf, double[,] pdops, double min, double max)
-        {
-            var bitmap = new SKBitmap(conf.xsize, conf.ysize);
 
-            for (int i = 0; i < conf.xsize; ++i)
+    public class Drawer
+    {
+        public static void DrawAll(Configurator conf, CalculationResults results)
+        {
+            Directory.CreateDirectory("output");
+            foreach ((string name, double[,] matrix) in results)
             {
-                for (int j = 0; j < conf.ysize; ++j)
+                Draw(conf, matrix, Path.Join("output", name + ".png"));
+            }
+        }
+
+        public static void Draw(Configurator conf, double[,] matrix, string fileName)
+        {
+            (double min, double max) = matrix.MinMax();
+            Console.WriteLine("{0}: {1}...{2}", fileName, min, max);
+
+            var bitmap = new SKBitmap(conf.XSize, conf.YSize);
+
+            for (int i = 0; i < conf.XSize; ++i)
+            {
+                for (int j = 0; j < conf.YSize; ++j)
                 {
                     SKColor color;
 
-                    if (pdops[i, j] == Calculator.StationPDOP)
+                    if (double.IsNaN(matrix[i, j]))
                     {
                         color = SKColors.White;
                     }
                     else
                     {
-                        double t = (pdops[i, j] - min) / (max - min);
+                        double t = (matrix[i, j] - min) / (max - min);
                         color = Map2Color(t);
                     }
 
@@ -30,7 +44,7 @@ namespace omp
 
             using (var image = SKImage.FromBitmap(bitmap))
             using (var data = image.Encode(SKEncodedImageFormat.Png, 100))
-            using (var stream = File.OpenWrite("out.png"))
+            using (var stream = File.OpenWrite(fileName))
             {
                 data.SaveTo(stream);
             }
@@ -64,4 +78,4 @@ namespace omp
             }
         }
     }
-}
+}

MatrixExtensions.cs

@@ -1,16 +1,18 @@
+using System.ComponentModel;
+
 namespace omp
 {
-    public static class MatrixOps
+    public static class MatrixExtensions
     {
-        public static void Print(double[,] a)
+        public static void Print(this double[,] a)
         {
             int rows = a.GetLength(0);
             int columns = a.GetLength(1);
 
             Console.WriteLine("Matrix of size {0}x{1}:", rows, columns);
-            for (int i = 0; i < rows; i++)
+            for (int i = 0; i < rows; ++i)
             {
-                for (int j = 0; j < columns; j++)
+                for (int j = 0; j < columns; ++j)
                 {
                     Console.Write("\t{0}", a[i, j]);
                 }
@@ -18,14 +20,35 @@ namespace omp
             }
         }
 
-        public static double Adj(double[,] a, int i, int j)
+        public static void SaveToFile(this double[,] a, string filePath)
+        {
+            int rows = a.GetLength(0);
+            int cols = a.GetLength(0);
+
+            using (var stream = new FileStream(filePath, FileMode.Create))
+            using (var writer = new BinaryWriter(stream))
+            {
+                unsafe
+                {
+                    fixed (double* ptr = &a[0, 0])
+                    {
+                        byte* bytePtr = (byte*)ptr;
+                        int byteCount = rows * cols * sizeof(double);
+
+                        writer.Write(new ReadOnlySpan<byte>(bytePtr, byteCount));
+                    }
+                }
+            }
+        }
+
+        public static double Adj(this double[,] a, int i, int j)
         {
             double sgn = ((i + j) % 2 == 0) ? 1 : -1;
 
             return sgn * Det(Minor(a, j, i));
         }
 
-        public static double[,] Minor(double[,] a, int i, int j)
+        public static double[,] Minor(this double[,] a, int i, int j)
         {
             int rows = a.GetLength(0);
             int columns = a.GetLength(1);
@@ -62,12 +85,14 @@ namespace omp
             return minor;
         }
 
-        public static double Det(double[,] a)
+        public static double Det(this double[,] a)
         {
             int n = a.GetLength(0);
 
             switch (n)
             {
+                case 1:
+                    return a[0, 0];
                 case 2:
                     return Det2(a);
                 case 3:
@@ -86,11 +111,12 @@ namespace omp
             }
         }
 
-        public static double Det2(double[,] a)
+        public static double Det2(this double[,] a)
         {
             return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0];
         }
 
+
         public static double Det3(double[,] a)
         {
             if (a.GetLength(0) != 3 || a.GetLength(1) != 3)
@@ -104,19 +130,18 @@ namespace omp
                 - a[0, 0] * a[1, 2] * a[2, 1]
                 - a[2, 2] * a[0, 1] * a[1, 0];
         }
-
+        public static double[,] MultTrans(this double[,] a)
-        public static double[,] MultTrans(double[,] a)
         {
             int rows = a.GetLength(0);
             int columns = a.GetLength(1);
 
             double[,] r = new double[rows, rows];
 
-            for (int i = 0; i < rows; i++)
+            for (int i = 0; i < rows; ++i)
             {
-                for (int j = 0; j < rows; j++)
+                for (int j = 0; j < rows; ++j)
                 {
-                    for (int k = 0; k < columns; k++)
+                    for (int k = 0; k < columns; ++k)
                     {
                         r[i, j] += a[i, k] * a[j, k];
                     }
@@ -125,5 +150,47 @@ namespace omp
 
             return r;
         }
+
+        public static (double Min, double Max) MinMax(this double[,] a)
+        {
+            double min = double.MaxValue;
+            double max = double.MinValue;
+
+            object MMLock = new();
+
+            Parallel.For(0, a.GetLength(0), i =>
+            {
+                double localMin = double.MaxValue;
+                double localMax = double.MinValue;
+
+                for (int j = 0; j < a.GetLength(1); ++j)
+                {
+                    if (a[i, j] < localMin)
+                    {
+                        localMin = a[i, j];
+                    }
+                    if (a[i, j] > localMax)
+                    {
+                        localMax = a[i, j];
+                    }
+                }
+
+                lock (MMLock)
+                {
+                    if (localMin < min)
+                    {
+                        min = localMin;
+                    }
+
+                    if (localMax > max)
+                    {
+                        max = localMax;
+                    }
+
+                }
+            });
+
+            return (min, max);
+        }
     }
-}
+}

Point.cs

@@ -1,10 +1,12 @@
 namespace omp
 {
-    public class Point(int x, int y, int z)
+    public class Point { }
+
+    public class Point3D(int x, int y, int z) : Point
     {
         public int x = x, y = y, z = z;
 
-        public double Distance(Point a)
+        public double DistanceTo(Point3D a)
         {
             double dx = x - a.x;
             double dy = y - a.y;
@@ -13,9 +15,27 @@ namespace omp
             return Math.Sqrt(dx * dx + dy * dy + dz * dz);
         }
 
-        public bool Equals(Point a)
+        public bool Equals(Point3D a)
         {
             return (x == a.x) && (y == a.y) && (z == a.z);
         }
     }
+
+    public class Point2D(int x, int y) : Point
+    {
+        public int x = x, y = y;
+
+        public double DistanceTo(Point2D a)
+        {
+            double dx = x - a.x;
+            double dy = y - a.y;
+
+            return Math.Sqrt(dx * dx + dy * dy);
+        }
+
+        public bool Equals(Point2D a)
+        {
+            return (x == a.x) && (x == a.y);
+        }
+    }
 }

Program.cs

@@ -4,14 +4,15 @@ namespace omp
     {
         public static void Main(string[] args)
         {
-            Configurator conf = new(args);
+            var configurator = new Configurator(args);
+            var calculator = CalculatorFactory.Create(configurator);
+            var results = calculator.Calculate();
+            Drawer.DrawAll(configurator, results);
 
-            double[,] pdops = new double[conf.xsize, conf.ysize];
+            foreach ((string name, double[,] arr) in results)
-            double min, max;
+            {
-
+                arr.SaveToFile(Path.Join("output", name + ".dat"));
-            Calculator.PDOPs(pdops, conf, out min, out max);
+            }
-
-            Drawer.Draw(conf, pdops, min, max);
         }
     }
 }

test2.txt

@@ -1,4 +1,3 @@
-500 500
+1000 1000
-1
+340 200
-100 499 0
+270 690
-400 499 0

test3.txt

@@ -1,5 +1,5 @@
 1000 1000
-0
+50
-236 134 0
+340 200 100
-147 863 0
+270 690 150
-915 591 0
+650 350 120

test30.txt

@@ -0,0 +1,5 @@
+100 100
+5
+34 20 10
+27 69 15
+65 35 12