fermi

sample.go (5332B)

---

 package sample
 
 import "math"
 import "sync"
 import rand "math/rand/v2"
 
 // https://pkg.go.dev/math/rand/v2
 
 type Src = *rand.Rand
 type func64 = func(Src) float64
 
 var global_r = rand.New(rand.NewPCG(uint64(1), uint64(2)))
 
 func Sample_unit_uniform(r Src) float64 {
 	return r.Float64()
 }
 
 func Sample_unit_normal(r Src) float64 {
 	return r.NormFloat64()
 }
 
 func Sample_uniform(start float64, end float64, r Src) float64 {
 	return Sample_unit_uniform(r)*(end-start) + start
 }
 
 func Sample_normal(mean float64, sigma float64, r Src) float64 {
 	return mean + Sample_unit_normal(r)*sigma
 }
 
 func Sample_lognormal(logmean float64, logstd float64, r Src) float64 {
 	return (math.Exp(Sample_normal(logmean, logstd, r)))
 }
 
 func Sample_normal_from_90_ci(low float64, high float64, r Src) float64 {
 	var normal90 float64 = 1.6448536269514727
 	var mean float64 = (high + low) / 2.0
 	var std float64 = (high - low) / (2.0 * normal90)
 	return Sample_normal(mean, std, r)
 }
 
 func Sample_to(low float64, high float64, r Src) float64 {
 	// Given a (positive) 90% confidence interval,
 	// returns a sample from a lognorma with a matching 90% c.i.
 	// Key idea: If we want a lognormal with 90% confidence interval [a, b]
 	// we need but get a normal with 90% confidence interval [log(a), log(b)].
 	// Then see code for Sample_normal_from_90_ci
 	var loglow float64 = math.Log(low)
 	var loghigh float64 = math.Log(high)
 	return math.Exp(Sample_normal_from_90_ci(loglow, loghigh, r))
 }
 
 func Sample_gamma(alpha float64, r Src) float64 {
 
 	// a simple method for generating gamma variables, marsaglia and wan tsang, 2001
 	// https://dl.acm.org/doi/pdf/10.1145/358407.358414
 	// see also the references/ folder
 	// note that the wikipedia page for the gamma distribution includes a scaling parameter
 	// k or beta
 	// https://en.wikipedia.org/wiki/gamma_distribution
 	// such that gamma_k(alpha, k) = k * gamma(alpha)
 	// or gamma_beta(alpha, beta) = gamma(alpha) / beta
 	// so far i have not needed to use this, and thus the second parameter is by default 1.
 
 	if alpha >= 1 {
 		var d, c, x, v, u float64
 		d = alpha - 1.0/3.0
 		c = 1.0 / math.Sqrt(9.0*d)
 
 		for {
 
 		InnerLoop:
 			for {
 				x = Sample_unit_normal(r)
 				v = 1.0 + c*x
 				if v > 0.0 {
 					break InnerLoop
 				}
 			}
 
 			v = v * v * v
 			u = Sample_unit_uniform(r)
 
 			if u < 1.0-0.0331*(x*x*x*x) { // Condition 1
 				// the 0.0331 doesn't inspire much confidence
 				// however, this isn't the whole story
 				// by knowing that Condition 1 implies condition 2
 				// we realize that this is just a way of making the algorithm faster
 				// i.e., of not using the logarithms
 				return d * v
 			}
 			if math.Log(u) < 0.5*(x*x)+d*(1.0-v+math.Log(v)) { // Condition 2
 				return d * v
 			}
 
 		}
 
 	} else {
 		return Sample_gamma(1.0+alpha, r) * math.Pow(Sample_unit_uniform(r), 1.0/alpha)
 	}
 }
 
 func Sample_beta(a float64, b float64, r Src) float64 {
 	gamma_a := Sample_gamma(a, r)
 	gamma_b := Sample_gamma(b, r)
 	return gamma_a / (gamma_a + gamma_b)
 }
 
 func Sample_mixture(fs []func64, weights []float64, r Src) float64 {
 
 	// fmt.Println("weights initially: ", weights)
 	var sum_weights float64 = 0
 	for _, weight := range weights {
 		sum_weights += weight
 	}
 
 	var total float64 = 0
 	var cumsummed_normalized_weights = append([]float64(nil), weights...)
 	for i, weight := range weights {
 		total += weight / sum_weights
 		cumsummed_normalized_weights[i] = total
 	}
 
 	var result float64
 	var flag int = 0
 	var p float64 = r.Float64()
 
 	for i, cnw := range cumsummed_normalized_weights {
 		if p < cnw {
 			result = fs[i](r)
 			flag = 1
 			break
 		}
 	}
 
 	if flag == 0 {
 		result = fs[len(fs)-1](r)
 	}
 	return result
 
 }
 
 func Sample_serially(f func64, n_samples int) []float64 {
 	xs := make([]float64, n_samples)
 	// var global_r = rand.New(rand.NewPCG(uint64(1), uint64(2)))
 	for i := 0; i < n_samples; i++ {
 		xs[i] = f(global_r)
 	}
 	return xs
 }
 
 func Sample_parallel(f func64, n_samples int) []float64 {
 	var num_threads = 16
 	var xs = make([]float64, n_samples)
 	var wg sync.WaitGroup
 	var h = n_samples / num_threads
 	wg.Add(num_threads)
 	for i := range num_threads {
 		var xs_i = xs[i*h : (i+1)*h]
 		go func(f func64) {
 			defer wg.Done()
 			var r = rand.New(rand.NewPCG(uint64(i), uint64(i+1)))
 			for i := range xs_i {
 				xs_i[i] = f(r)
 			}
 		}(f)
 	}
 
 	wg.Wait()
 	return xs
 
 }
 
 /*
 func main() {
 
 	var p_a float64 = 0.8
 	var p_b float64 = 0.5
 	var p_c float64 = p_a * p_b
 	ws := [4](float64){1 - p_c, p_c / 2, p_c / 4, p_c / 4}
 
 	Sample_0 := func(r Src) float64 { return 0 }
 	Sample_1 := func(r Src) float64 { return 1 }
 	Sample_few := func(r Src) float64 { return Sample_to(1, 3, r) }
 	Sample_many := func(r Src) float64 { return Sample_to(2, 10, r) }
 	fs := [4](func64){Sample_0, Sample_1, Sample_few, Sample_many}
 
 	model := func(r Src) float64 { return Sample_mixture(fs[0:], ws[0:], r) }
 	n_samples := 1_000_000
 	xs := Sample_parallel(model, n_samples)
 	var avg float64 = 0
 	for _, x := range xs {
 		avg += x
 	}
 	avg = avg / float64(n_samples)
 	fmt.Printf("Average: %v\n", avg)
 	/*
 		  // Without concurrency:
 			n_samples := 1_000_000
 			var r = rand.New(rand.NewPCG(uint64(1), uint64(2)))
 			var avg float64 = 0
 			for i := 0; i < n_samples; i++ {
 				avg += Sample_mixture(fs[0:], ws[0:], r)
 			}
 			avg = avg / float64(n_samples)
 			fmt.Printf("Average: %v\n", avg)
 }
 */