Spatial index (doesn't work yet)

Icfpc2023.Tests/Tests.fs

@@ -3,7 +3,6 @@ module Tests
 open System.IO
 
 open Icfpc2023
-open Icfpc2023.Scoring
 open Xunit
 
 [<Fact>]
@@ -91,3 +90,19 @@ let ``Stadium detects intersection when blocking musician is right on the line``
         }
     let blockingMusician = PointD(1100.0, 150.0)
     Assert.True(stadium.Contains blockingMusician)
+
+[<Fact>]
+let ``CalculateScore for problem 1``(): unit =
+    let problem = JsonDefs.ReadProblemFromFile(Path.Combine(DirectoryLookup.problemsDir, "1.json"))
+    let solution = JsonDefs.ReadSolutionFromFile(Path.Combine(DirectoryLookup.solutionsDir, "1.json"))
+    let score = Scoring.CalculateScore problem solution
+    let meta = JsonDefs.ReadSolutionMetadataFromFile(Path.Combine(DirectoryLookup.solutionsDir, "1.meta.json"))
+    Assert.Equal(meta.Score, score)
+
+[<Fact>]
+let ``CalculateScore for problem 30``(): unit =
+    let problem = JsonDefs.ReadProblemFromFile(Path.Combine(DirectoryLookup.problemsDir, "30.json"))
+    let solution = JsonDefs.ReadSolutionFromFile(Path.Combine(DirectoryLookup.solutionsDir, "30.json"))
+    let score = Scoring.CalculateScore problem solution
+    let meta = JsonDefs.ReadSolutionMetadataFromFile(Path.Combine(DirectoryLookup.solutionsDir, "30.meta.json"))
+    Assert.Equal(meta.Score, score)

Icfpc2023/Geometry.fs

@@ -1,5 +1,7 @@
 namespace Icfpc2023
 
+open System
+
 [<Struct>]
 type PointD =
     | PointD of double * double
@@ -31,6 +33,8 @@ type Line =
 
         abs((x2 - x1) * (y1 - y0) - (x1 - x0) * (y2 - y1)) / sqrt((x2 - x1) ** 2.0 + (y2 - y1) ** 2)
 
+#nowarn "3391"
+
 [<Struct>]
 type Rectangle =
     {
@@ -101,3 +105,66 @@ type Stadium =
         // doesn't exceed the length of the segment
         let segment_length = this.Center1.DistanceTo(this.Center2)
         closest.DistanceTo(this.Center1) < segment_length && closest.DistanceTo(this.Center2) < segment_length
+
+    static member NormalizeVector(v: PointD): PointD =
+        let (PointD(x, y)) = v
+        let length = sqrt(x ** 2.0 + y ** 2.0)
+        PointD(x / length, y / length)
+
+    static member RotateVector (radians: double) (v: PointD): PointD =
+        let (PointD(x, y)) = v
+        let cos = cos radians
+        let sin = sin radians
+        PointD(x * cos - y * sin, x * sin + y * cos)
+
+    member this.RectanglePoints = [|
+        this.Center1 + (Stadium.NormalizeVector(this.Center2 - this.Center1) |> Stadium.RotateVector(Math.PI/2.0)) * this.Radius
+        this.Center1 + (Stadium.NormalizeVector(this.Center2 - this.Center1) |> Stadium.RotateVector(-Math.PI/2.0)) * this.Radius
+        this.Center2 + (Stadium.NormalizeVector(this.Center2 - this.Center1) |> Stadium.RotateVector(Math.PI/2.0)) * this.Radius
+        this.Center2 + (Stadium.NormalizeVector(this.Center2 - this.Center1) |> Stadium.RotateVector(-Math.PI/2.0)) * this.Radius
+    |]
+
+    /// https://stackoverflow.com/a/10965077/2684760
+    static member DoPolygonsIntersect(a: PointD[], b: PointD[]): bool =
+        let mutable foundNoIntersection = false
+        for polygon in [|a; b|] do
+            for i1 in 0 .. polygon.Length - 1 do
+                if not foundNoIntersection then
+                    let i2 = (i1 + 1) % polygon.Length
+                    let p1 = polygon[i1]
+                    let p2 = polygon[i2]
+
+                    let normal = PointD(p2.Y - p1.Y, p1.X - p2.X)
+
+                    let mutable minA = Nullable<double>()
+                    let mutable maxA = Nullable<double>()
+                    for p in a do
+                        let projected = normal.X * p.X + normal.Y * p.Y
+                        if not minA.HasValue || projected < minA.Value then
+                            minA <- projected
+                        if not maxA.HasValue || projected > maxA.Value then
+                            maxA <- projected
+
+                    let mutable minB = Nullable<double>()
+                    let mutable maxB = Nullable<double>()
+                    for p in b do
+                        let projected = normal.X * p.X + normal.Y * p.Y
+                        if not minB.HasValue || projected < minB.Value then
+                            minB <- projected
+                        if not maxB.HasValue || projected > maxB.Value then
+                            maxB <- projected
+
+                    if Nullable.Compare(maxA, maxB) < 0 || Nullable.Compare(maxB, minA) < 0 then
+                        foundNoIntersection <- true
+
+        if foundNoIntersection then false else true
+
+    member this.RectangularPartIntersectsWith(rect: Rectangle): bool =
+        let p1Points = [|
+            rect.BottomLeft
+            rect.BottomLeft + PointD(rect.Width, 0.0)
+            rect.BottomLeft + PointD(rect.Width, rect.Height)
+            rect.BottomLeft + PointD(0.0, rect.Height)
+        |]
+        let p2Points = this.RectanglePoints
+        Stadium.DoPolygonsIntersect(p1Points, p2Points)

Icfpc2023/Scoring.fs

@@ -9,18 +9,58 @@ type private Musician = {
     Volume: double
 }
 
+type SectoredIndex(problem: Problem) =
+    let sectorSize = 50.0
+    let sectors =
+        let w = int <| ceil(problem.StageWidth / sectorSize)
+        let h = int <| ceil(problem.StageHeight / sectorSize)
+        Array2D.init w h (fun _ _ -> ResizeArray<PointD>())
+
+    let getSectorSquares(stadium: Stadium) =
+        let stadium' =
+            { stadium with
+                Center1 = stadium.Center1 - problem.StageBottomLeft
+                Center2 = stadium.Center2 - problem.StageBottomLeft
+            }
+
+        let w = Array2D.length1 sectors
+        let h = Array2D.length2 sectors
+        seq {
+            for x in 0..w - 1 do
+                for y in 0..h - 1 do
+                    let square = {
+                        BottomLeft = PointD(float x * sectorSize, float y * sectorSize)
+                        Width = sectorSize
+                        Height = sectorSize
+                    }
+                    if stadium'.RectangularPartIntersectsWith square then
+                        yield struct(x, y)
+        }
+
+    member _.Add(point: PointD): unit =
+        let sector = sectors[int ((point.X - problem.StageBottomLeft.X) / sectorSize), int ((point.Y- problem.StageBottomLeft.Y) / sectorSize)]
+        sector.Add point
+
+    member _.GetPointsIn(s: Stadium): PointD seq =
+        getSectorSquares s
+        |> Seq.collect(fun (struct(x, y)) -> sectors[x, y])
+
 let private CalculateAttendeeMusicianScore(attendee: Attendee, musician: Musician, closeness: ClosenessFactor): Score =
     let d_squared = (attendee.X - musician.Location.X) ** 2.0 + (attendee.Y - musician.Location.Y) ** 2.0
     musician.Volume * ceil(closeness * 1_000_000.0 * attendee.Tastes[musician.Instrument] / d_squared)
 
-let private AnyOtherMusicianBlocksSound (musicians: Musician[]) (attendee: Attendee) (mIndex: int): bool =
-    let musician = PointD(musicians[mIndex].Location.X, musicians[mIndex].Location.Y)
+let private AnyOtherMusicianBlocksSound(index: SectoredIndex,
+                                        attendee: Attendee,
+                                        musician: Musician): bool =
+    let musician = musician.Location
     let attendee = PointD(attendee.X, attendee.Y)
     let blockZone = { Center1 = musician; Center2 = attendee; Radius = 5.0 }
+    let candidates = index.GetPointsIn blockZone |> Seq.toArray // TODO REMOVE TOARRAY
 
-    Seq.indexed musicians
-    |> Seq.filter(fun (i, _) -> i <> mIndex)
-    |> Seq.exists(fun (_, m) -> blockZone.Contains(m.Location))
+    candidates
+    |> Seq.filter(fun p -> p <> musician)
+    |> Seq.tryHead
+    |> Option.isSome
 
 let private AnyPillarBlocksSound (pillars: Pillar[]) (musician: Musician) (attendee: Attendee): bool =
     let musician = PointD(musician.Location.X, musician.Location.Y)
@@ -34,10 +74,14 @@ let private AnyPillarBlocksSound (pillars: Pillar[]) (musician: Musician) (atten
             }
             blockZone.Contains(p.Center))
 
-let private CalculateAttendeeScore (pillars: Pillar[]) (musicians: Musician[]) (closenessFactors: ClosenessFactor[]) (attendee: Attendee): Score =
+let private CalculateAttendeeScore(pillars: Pillar[],
+                                   musicians: Musician[],
+                                   index: SectoredIndex,
+                                   closenessFactors: ClosenessFactor[],
+                                   attendee: Attendee): Score =
     Seq.indexed musicians
     |> Seq.sumBy(fun (i, musician) ->
-        if AnyOtherMusicianBlocksSound musicians attendee i then 0.0
+        if AnyOtherMusicianBlocksSound(index, attendee, musicians[i]) then 0.0
         else if AnyPillarBlocksSound pillars musicians.[i] attendee then 0.0
         else CalculateAttendeeMusicianScore(attendee, musician, closenessFactors.[i])
     )
@@ -105,7 +149,11 @@ let CalculateScore(problem: Problem) (solution: Solution): Score =
     if tooClose
     then BadScore
     else
-        problem.Attendees |> Array.sumBy(CalculateAttendeeScore problem.Pillars musicians closeness_factors)
+        let index = SectoredIndex problem
+        musicians |> Seq.iter(fun m -> index.Add m.Location)
+
+        problem.Attendees
+        |> Array.sumBy(fun a -> CalculateAttendeeScore(problem.Pillars, musicians, index, closeness_factors, a))
 
 let CalculateNoBlockingScore(problem: Problem) (solution: IPartialSolution): Score =
     let musicians =