savepoint

makefile

@@ -1,2 +1,9 @@
 dev:
 	go run probppl.go
+
+prod:
+	go build -o probppl
+	./probppl
+
+record:
+	make prod > record.txt

probppl.go

@@ -5,36 +5,36 @@ import (
 	"git.nunosempere.com/NunoSempere/probppl/choose"
 	"math"
 	rand "math/rand/v2"
+	"sync"
 )
 
 type src = *rand.Rand
-type pplKnownDistrib = map[int64]float64
+type IntProb struct {
+	N int64
+	p float64
+}
+type IntProbs = []IntProb
+type IntProbsWeights struct {
+	IntProbs IntProbs
+	w        int64
+}
 
-func generatePeopleKnownDistribution(r src) map[int64]float64 {
-	mapping := make(map[int64]float64)
+func generatePeopleKnownDistribution(r src) IntProbs {
+
+	var probabilities IntProbs
 	sum := 0.0
 
-	// Consider zero case separately
-	/*p0 := r.Float64()
-	mapping[0.0] = p0
-	sum += p0
-	*/
-
-	// Consider successive exponents of 1.5
-	num := 16.0
-	base := 2.0
-	for i := 1; i < 8; i++ {
-		num = num * base
+	for i := 16; i <= 2048; i *= 2 {
 		p := r.Float64()
-		mapping[int64(num)] = p
+		probabilities = append(probabilities, IntProb{N: int64(i), p: p})
 		sum += p
 	}
 
-	for key, value := range mapping {
-		mapping[key] = value / sum
+	for i := range probabilities {
+		probabilities[i].p = probabilities[i].p / sum
 	}
 
-	return mapping
+	return probabilities
 }
 
 func chooseWrapper(n int64, k int64) int64 {
@@ -75,13 +75,13 @@ func getMatchesDrawGivenNPeopleKnown(n int64, r src) int64 {
 	≥3: 9.5% | 14
 */
 
-func drawFromDistributionWithReplacement(d pplKnownDistrib, r src) int64 {
+func drawFromDistributionWithReplacement(d IntProbs, r src) int64 {
 	pp := r.Float64()
 	sum := 0.0
-	for i, p := range d {
-		sum += p
+	for i := range d {
+		sum += d[i].p
 		if pp <= sum {
-			return int64(i)
+			return int64(d[i].N) // this introduces some non-determinism, as order of maps in go isn't guaranteed
 		}
 
 	}
@@ -95,7 +95,7 @@ func aboutEq(a int64, b int64) bool {
 	return ((-h) <= (a - b)) && ((a - b) <= h)
 }
 
-func draw148PplFromDistributionAndCheck(d pplKnownDistrib, r src, show bool) int64 {
+func draw148PplFromDistributionAndCheck(d IntProbs, r src, show bool) int64 {
 
 	count := make(map[int64]int64)
 	count[0] = 0
@@ -118,9 +118,9 @@ func draw148PplFromDistributionAndCheck(d pplKnownDistrib, r src, show bool) int
 	}
 }
 
-func getUnnormalizedBayesianUpdateForDistribution(d pplKnownDistrib, r src) int64 {
+func getUnnormalizedBayesianUpdateForDistribution(d IntProbs, r src) int64 {
 	var sum int64 = 0
-	n := 100
+	n := 30_000
 	for i := 0; i < n; i++ {
 		/* if i%1000 == 0 {
 			fmt.Println(i)
@@ -134,22 +134,52 @@ func getUnnormalizedBayesianUpdateForDistribution(d pplKnownDistrib, r src) int6
 
 func main() {
 
-	var r = rand.New(rand.NewPCG(uint64(1), uint64(2)))
+	n_dists := 30_000
+	var dists = make([]IntProbsWeights, n_dists)
+
+	// Prepare for concurrency
+	num_threads := 32
+	var wg sync.WaitGroup
+	wg.Add(num_threads)
+
+	for i := range num_threads {
+		go func() {
+			defer wg.Done()
+			var r = rand.New(rand.NewPCG(uint64(i), uint64(i+1)))
+			for j := i * (n_dists / num_threads); j < (i+1)*(n_dists/num_threads); j++ {
+				people_known_distribution := generatePeopleKnownDistribution(r)
+				result := getUnnormalizedBayesianUpdateForDistribution(people_known_distribution, r)
+				/*
+					if i%10 == 0 {
+						fmt.Printf("%d/%d\n", i, n_dists)
+					}
+				*/
+				if result > 0 {
+					dists[j] = IntProbsWeights{IntProbs: people_known_distribution, w: result}
+				}
+			}
+
+		}()
 
-	sum := int64(0)
-	for i := 0; i < 1000; i++ {
+	}
 
-		people_known_distribution := generatePeopleKnownDistribution(r)
-		// fmt.Println(people_known_distribution)
-		result := getUnnormalizedBayesianUpdateForDistribution(people_known_distribution, r)
-		fmt.Println(i)
-		if result > 0 {
-			fmt.Println(people_known_distribution)
-			fmt.Println(result)
-		}
-		sum += result
-		// fmt.Println(result)
+	wg.Wait()
+	// Now calculate the posterior
+	sum_weights := int64(0)
+	for _, dist := range dists {
+		sum_weights += dist.w
 	}
-	fmt.Println(sum)
 
+	for i := int64(16); i <= 2048; i *= 2 {
+		p := 0.0
+		for _, dist := range dists {
+			for _, int_prob := range dist.IntProbs {
+				if int_prob.N == i {
+					p += float64(dist.w) * int_prob.p
+				}
+			}
+		}
+		p = p / float64(sum_weights)
+		fmt.Printf("%d: %f\n", i, p)
+	}
 }