```﻿using System;
using System.Linq;
using System.Numerics;
using System.Text;
using Clustering;
using HilbertTransformation;
using HilbertTransformation.Random;
using HilbertTransformationTests.Data;
using NUnit.Framework;

namespace HilbertTransformationTests
{
/// <summary>
/// Test the Hilbert Transformation.
///
///    The transformation:
///
///       A. Hilbert Index to HilbertPoint to N-Dimensional coordinates
///
///          int bits = ???;       // Pick so that 2^bits exceeds the larges value in any coordinate.
///          int dimensions = ???; // Number of dimensions for the point.
///          var index1 = new BigInteger(...);
///          var hPoint1 = new HilbertPoint(index1, dimensions, bits);
/// 	     uint[] coordinates = hPoint.Coordinates;
///
///       B. Coordinates to Hilbert Index
///
///          var hPoint2 = new HilbertPoint(coordinates, bits);
///          BigInteger index2 = hPoint2.Index;
/// </summary>
[TestFixture]
public class HilbertPointTests
{
[Test]
{
}

[Test]
{
}

[Test]
{
}

/// <summary>
/// Verify that the square Cartesian distance calculation is correct.
/// </summary>
[Test]
public void MeasureDistanceSquared()
{
const int dims = 100;
var p1 = new uint[dims];
var p2 = new uint[dims];
var expectedSquareDistance = 0L;
for (var i = 0; i < dims; i++)
{
p1[i] = (uint)(i % 37) * 10;
p2[i] = (uint)(i % 18) * 17;
long delta = (long)p1[i] - (long)p2[i];
expectedSquareDistance += delta * delta;
}
var up1 = new UnsignedPoint(p1);
var up2 = new UnsignedPoint(p2);
var actualSquareDistance = up1.Measure(up2);
Assert.AreEqual(expectedSquareDistance, actualSquareDistance, "Distances do not match");
}

/// <summary>
/// The proof of this test is studying the Console output by eye.
/// </summary>
[Test]
public void CartesianToHilbert_Dim2Bits2()
{
var bits = 2;
var size = 1 << bits;
var sb = new StringBuilder();
for (var row = 0; row < size; row++)
{
for (var column = 0; column < size; column++)
{
var cartesianPoint = new int[] { row, column };
var hilbertPoint = new HilbertPoint(cartesianPoint, bits);
var hilbertIndex = hilbertPoint.HilbertIndex;
sb.Append("Cart = [")
.Append(string.Join(",", cartesianPoint))
.Append("] Hilbert = ")
.Append(hilbertIndex.ToString())
.AppendLine();
}
}
var diagnostic = sb.ToString();
Console.WriteLine(diagnostic);
}

/// <summary>
/// Verify the transformation in both directions, from 1-Dimensional index to N-dimensional point and back.
///
/// This evaluates 2^(dims*bits) points, so be careful or the test will run for a long time and consume a lot of memory.
/// </summary>
/// <param name="dims">Dimensions for each point.</param>
/// <param name="bits">Bits per dimension.</param>
public void AdjacentPointsCase(int dims, int bits)
{
var points = new HilbertPoint[1 << (bits*dims)];
for (var i = 0; i < points.Length; i++)
{
var hilbertIndex = new BigInteger(i);
points[i] = new HilbertPoint(hilbertIndex, dims, bits);
if (i > 0)
{
var p1 = points[i - 1];
var p2 = points[i];
string.Format("Points {0} and {1}",
FormatPoint(p1), FormatPoint(p2)));
}
AssertPointMapsToHilbertIndex(points[i].Coordinates, hilbertIndex, dims, bits);
}
}

static void AssertPointMapsToHilbertIndex(uint[] coordinates, BigInteger hilbertIndex, int dims, int bits)
{
var hPoint = new HilbertPoint(coordinates, bits);
Assert.AreEqual(hilbertIndex, hPoint.HilbertIndex , \$"Coordinates {UintsToString(coordinates)} do not map back to the expected Hilbert Index {hilbertIndex.ToString()}, but instead {hPoint.HilbertIndex}");
}

/// <summary>
/// Test if two points are adjacent, meaning that only a single coordiante differs between them and
/// the difference in coordinate value is exactly one.
/// </summary>
/// <returns><c>true, if points are adjacent, false otherwise.
/// <param name="p1">First point.</param>
/// <param name="p2">Second point.</param>
static bool ArePointsAdjacent(HilbertPoint p1, HilbertPoint p2)
{
var maxCoordinateDistance = 0;
var differentDimensionsCount = 0;
for (var dim = 0; dim < p1.Dimensions; dim++)
{
var diff = Math.Abs(p1[dim] - p2[dim]);
if (diff != 0)
{
differentDimensionsCount++;
maxCoordinateDistance = Math.Max(diff, maxCoordinateDistance);
}
}
return maxCoordinateDistance == 1 && differentDimensionsCount == 1;
}

/// <summary>
/// Pretty print a HilbertPoint.
/// </summary>
/// <returns>Formatted point.</returns>
/// <param name="p">Point to pretty print.</param>
static string FormatPoint(HilbertPoint p)
{
return string.Format("Index: {0} Coords: [{1}]", p.HilbertIndex, string.Join(",", p.Coordinates));
}

/// <summary>
/// Verify that two vectors have the same values, except that one coordinate differs form another by one.
/// </summary>
/// <returns><c>true, if only one coordinate differs, false otherwise.
/// <param name="u">First vector to compare.</param>
/// <param name="v">Second vector to compare.</param>
private static bool DifferByOne(uint[] u, uint[] v)
{
var dims = u.Length;
var sum = 0U;
for (var dim = 0; dim < dims; dim++)
{
sum += u[dim] > v[dim] ? u[dim] - v[dim] : v[dim] - u[dim];
}
return sum == 1U;
}

/// <summary>
/// Pretty print a uint array.
/// </summary>
/// <returns>Formatted string.</returns>
/// <param name="uVec">Vector of points.</param>
private static string UintsToString(uint[] uVec)
{
var sb = new StringBuilder();
foreach (var u in uVec)
{
if (sb.Length == 0)
sb.Append("[");
else
sb.Append(",");
sb.Append(u);
}
sb.Append("]");
return sb.ToString();
}

/// <summary>
/// Make sure that SquareDistanceCompare returns the correct result, and if every test succeeds,
/// report on the Average percentage of the time that the SquareDistanceCompare optimization could be used.
/// </summary>
[Test]
public void SquareDistanceCompareValidation()
{
var repeats = 25;
var avg = Enumerable.Range(1, repeats).Select(i => SquareDistanceCompareValidationCase(10)).Average();
Console.WriteLine(\$"For {repeats} tries, Average usage of optimization is {avg} %");
}

/// <summary>
/// Vary the number of triangulation points used when comparing distances
/// to see where the point of diminishing returns is.
/// </summary>
[Test]
public void SquareDistanceCompareVaryTriangulation()
{
/*
Sample output:
For 25 tries, Triangulating with 2 points, Average usage of optimization is 37.1748 %
For 25 tries, Triangulating with 4 points, Average usage of optimization is 40.8068 %
For 25 tries, Triangulating with 6 points, Average usage of optimization is 44.0816 %
For 25 tries, Triangulating with 8 points, Average usage of optimization is 45.6396 %
For 25 tries, Triangulating with 10 points, Average usage of optimization is 47.7976 %
For 25 tries, Triangulating with 12 points, Average usage of optimization is 48.8452 %
For 25 tries, Triangulating with 14 points, Average usage of optimization is 50.2412 %
For 25 tries, Triangulating with 16 points, Average usage of optimization is 50.6484 %
For 25 tries, Triangulating with 18 points, Average usage of optimization is 50.52 %
*/

var repeats = 25;
var triangulationValues = new[] { 2, 4, 6, 8, 10, 12, 14, 16, 18 };
var report = "";
foreach (var numTriangulationPoints in triangulationValues)
{
var avg = Enumerable.Range(1, repeats).Select(i => SquareDistanceCompareValidationCase(numTriangulationPoints)).Average();
var record = \$"    For {repeats} tries, Triangulating with {numTriangulationPoints} points, Average usage of optimization is {avg} %\n";
report += record;
Console.WriteLine(report);
}
Console.WriteLine("\n\nFinal report:\n");
Console.WriteLine(report);
}

public double SquareDistanceCompareValidationCase(int numTriangulationPoints)
{
var correctResult = 0;
var wrongResult = 0;
var totalComparisons = 10000;
var extraShortTrianagulatable = 0;
var extraShortNotTrianagulatable = 0;
var shortTrianagulatable = 0;
var shortNotTrianagulatable = 0;
var longTrianagulatable = 0;
var longNotTrianagulatable = 0;

// 1. Make test data.
var bitsPerDimension = 10;
var data = new GaussianClustering
{
ClusterCount = 100,
Dimensions = 100,
MaxCoordinate = (1 << bitsPerDimension) - 1,
MinClusterSize = 50,
MaxClusterSize = 150
};
var clusters = data.MakeClusters();

// 2. Create HilbertIndex for points.
var hIndex = new HilbertIndex(clusters, bitsPerDimension);
hIndex.SetTriangulation(numTriangulationPoints);

// 3. Deduce the characteristic distance.
var counter = new ClusterCounter
{
OutlierSize = 5,
NoiseSkipBy = 10
};
var count = counter.Count(hIndex.SortedPoints);
var mergeDistance = count.MaximumSquareDistance;
var longDistance = 5 * mergeDistance;

// 4. Select random pairs of the HilbertPoints points and see how many distance comparisons yield the correct result.
var rng = new FastRandom();
var points = hIndex.SortedPoints.ToList();

for (var i = 0; i < totalComparisons; i++)
{
var p1 = points[rng.Next(points.Count)];
var p2 = points[rng.Next(points.Count)];
var d = p1.Measure(p2);
if (d.CompareTo(mergeDistance) == p1.SquareDistanceCompare(p2, mergeDistance))
correctResult++;
else
wrongResult++;

if (d.CompareTo(longDistance) == p1.SquareDistanceCompare(p2, longDistance))
correctResult++;
else
wrongResult++;

if (p1.Triangulatable(p2, mergeDistance/2))
extraShortTrianagulatable++;
else
extraShortNotTrianagulatable++;

if (p1.Triangulatable(p2, mergeDistance))
shortTrianagulatable++;
else
shortNotTrianagulatable++;

if (p1.Triangulatable(p2, longDistance))
longTrianagulatable++;
else
longNotTrianagulatable++;
}
var extraShortPct = 100.0 * extraShortTrianagulatable / (extraShortTrianagulatable + extraShortNotTrianagulatable);
var shortPct = 100.0 * shortTrianagulatable / (shortTrianagulatable + shortNotTrianagulatable);
var longPct = 100.0 * longTrianagulatable / (longTrianagulatable + longNotTrianagulatable);
Console.WriteLine(\$"Triangulatable? \n    XS: {extraShortPct} % \n    Short: {shortPct} % Yes {shortTrianagulatable}, No {shortNotTrianagulatable}\n    Long: {longPct} % Yes {longTrianagulatable}, No {longNotTrianagulatable}");
Assert.AreEqual(wrongResult, 0, \$"{correctResult} correct, {wrongResult} wrong");

return shortPct;
}

[Test]
public void SquareDistanceCompareAverage()
{
// Example output:
//    After 200 trials, Optimizations were possible on Average 29.2514 %, with Min 20.14 % and Max 39.26 %
var sumPercent = 0.0;
var count = 200;
var maxPercent = 0.0;
var minPercent = 100.0;
for (var i = 0; i < count; i++)
{
var percent = SquareDistanceCompareOptimizableCase(10000, false);
maxPercent = Math.Max(maxPercent, percent);
minPercent = Math.Min(minPercent, percent);
sumPercent += percent;
}
var avgPercent = sumPercent / count;
var message = \$"After {count} trials, Optimizations were possible on Average {avgPercent} %, with Min {minPercent} % and Max {maxPercent} %";
Console.WriteLine(message);
Assert.GreaterOrEqual(avgPercent, 25.0, message);
}

[Test]
public void SquareDistanceExtendedCompareAverage()
{
// Example output:
//    After 200 trials, Optimizations were possible on Average 37.959 %, with Min 26.17 % and Max 46.72 %
var sumPercent = 0.0;
var count = 200;
var maxPercent = 0.0;
var minPercent = 100.0;
for (var i = 0; i < count; i++)
{
var percent = SquareDistanceCompareOptimizableCase(10000, true);
maxPercent = Math.Max(maxPercent, percent);
minPercent = Math.Min(minPercent, percent);
sumPercent += percent;
}
var avgPercent = sumPercent / count;
var message = \$"After {count} trials, Optimizations were possible on Average {avgPercent} %, with Min {minPercent} % and Max {maxPercent} %";
Console.WriteLine(message);
Assert.GreaterOrEqual(avgPercent, 25.0, message);
}

/// <summary>
/// UnsignedPoint.SquareDistanceCompare has an optimization. This tests how often this optimization
/// can be exploited in a realistic test. The comparison will be against an estimated characteristic distance
/// between points. This distance is assumed to be close enough to trigger two points to be merged into a single cluster.
/// </summary>
private double SquareDistanceCompareOptimizableCase(int totalComparisons, bool useExtendedOptimization = false)
{
// 1. Make test data.
var bitsPerDimension = 10;
var data = new GaussianClustering
{
ClusterCount = 100,
Dimensions = 100,
MaxCoordinate = (1 << bitsPerDimension) - 1,
MinClusterSize = 50,
MaxClusterSize = 150
};
var clusters = data.MakeClusters();

// 2. Create HilbertIndex for points.
var hIndex = new HilbertIndex(clusters, bitsPerDimension);

// 3. Deduce the characteristic distance.
var counter = new ClusterCounter
{
OutlierSize = 5,
NoiseSkipBy = 10
};
var count = counter.Count(hIndex.SortedPoints);
var mergeDistance = count.MaximumSquareDistance;
var longDistance = 5 * mergeDistance;

// 4. Select random pairs of points and see how many distance comparisons can exploit the optimization.
var rng = new FastRandom();
var points = clusters.Points().ToList();
var ableToUseOptimizationsAtShortDistance = 0;
var ableToUseOptimizationsAtLongDistance = 0;

for (var i = 0; i < totalComparisons; i++)
{
var p1 = points[rng.Next(points.Count)];
var p2 = points[rng.Next(points.Count)];
if (useExtendedOptimization)
{
if (IsExtendedDistanceOptimizationUsable(p1, p2, mergeDistance, bitsPerDimension))
ableToUseOptimizationsAtShortDistance++;
if (IsExtendedDistanceOptimizationUsable(p1, p2, longDistance, bitsPerDimension))
ableToUseOptimizationsAtLongDistance++;
}
else
{
if (IsDistanceOptimizationUsable(p1, p2, mergeDistance))
ableToUseOptimizationsAtShortDistance++;
if (IsDistanceOptimizationUsable(p1, p2, longDistance))
ableToUseOptimizationsAtLongDistance++;
}
}
var percentOptimizable = 100.0 * ableToUseOptimizationsAtShortDistance / totalComparisons;
var percentOptimizableLongDistance = 100.0 * ableToUseOptimizationsAtLongDistance / totalComparisons;
var message = \$"Comparisons were {percentOptimizable} % Optimizable at short distance, {percentOptimizableLongDistance} % at long distance";
Console.WriteLine(message);
return percentOptimizable;
}

/// <summary>
/// UnsignedPoint.SquareDistanceCompare has an optimization.
/// This tests if that optimization can be used in a given case.
/// </summary>
/// <returns><c>true, if distance optimization is usable, false otherwise.
/// <param name="p1">First point to compare.</param>
/// <param name="p2">Second point to compare.</param>
/// <param name="squareDistance">Test if the distance between the points is less than, equal to or greater than this given distance.</param>
private bool IsDistanceOptimizationUsable(UnsignedPoint p1, UnsignedPoint p2, long squareDistance)
{
var delta = p1.Magnitude - p2.Magnitude;
var low = (long)Math.Floor(delta * delta);
if (squareDistance < low) return true;

var high = p1.SquareMagnitude + p2.SquareMagnitude;
return (squareDistance > high);
}

/// <summary>
/// UnsignedPoint.SquareDistanceCompare has an optimization.
/// This tests if an extension of that optimization can be used in a given case.
/// </summary>
/// <returns><c>true, if distance optimization is usable, false otherwise.
/// <param name="p1">First point to compare.</param>
/// <param name="p2">Second point to compare.</param>
/// <param name="squareDistance">Test if the distance between the points is less than, equal to or greater than this given distance.</param>
/// <param name="bitsPerDimension">Number of bits needed to represent the larges value of any coordinate of any point.</param>
private bool IsExtendedDistanceOptimizationUsable(UnsignedPoint p1, UnsignedPoint p2, long squareDistance, int bitsPerDimension)
{
if (IsDistanceOptimizationUsable(p1, p2, squareDistance))
return true;

var maxCoordinate = (1 << bitsPerDimension) - 1;
var cornerSqDistance1 = 0L;
var cornerSqDistance2 = 0L;
for (var d = 0; d < p1.Dimensions; d++)
{
long delta;
delta = maxCoordinate - p1[d];
cornerSqDistance1 += delta * delta;
delta = maxCoordinate - p2[d];
cornerSqDistance2 += delta * delta;
}
var cornerDistance1 = Math.Sqrt(cornerSqDistance1);
var cornerDistance2 = Math.Sqrt(cornerSqDistance2);

var delta2 = cornerDistance1 - cornerDistance2;
var low = (long)Math.Floor(delta2 * delta2);
if (squareDistance < low) return true;

var high = cornerSqDistance1 + cornerSqDistance2;
if (squareDistance > high) return true;

return false;
}
}
}

```