diff --git a/f.go b/f.go index c8e946f..6bceac8 100644 --- a/f.go +++ b/f.go @@ -4,6 +4,7 @@ import ( "bufio" "errors" "fmt" + "git.nunosempere.com/NunoSempere/fermi/sample" "math" "os" "strconv" @@ -26,14 +27,16 @@ type Lognormal struct { } func (l Lognormal) Samples() []float64 { - + sampler := func(r sample.Src) float64 { return sample.Sample_to(l.low, l.high, r) } + return sample.Sample_parallel(sampler, 1_000_000) } // Actually, I should look up how do do a) enums in go, b) union types -type Lognormal struct { +/*type Lognormal struct { low float64 high float64 } +*/ type Dist struct { Type string @@ -176,6 +179,11 @@ func prettyPrintDist(dist Dist) { /* Main event loop */ func main() { + + sample_0 := func(r sample.Src) float64 { return 0 } + x := sample.Sample_parallel(sample_0, 10) + fmt.Printf("%v\n", x) + reader := bufio.NewReader(os.Stdin) init_dist := Dist{Type: "Lognormal", Lognormal: Lognormal{low: 1, high: 1}, Samples: nil} // Could also just be a scalar old_dist := init_dist diff --git a/go.mod b/go.mod index 9232354..4ef7fee 100644 --- a/go.mod +++ b/go.mod @@ -1,3 +1,3 @@ -module git.nunosempere.com/NunoSempere/fermi.git +module git.nunosempere.com/NunoSempere/fermi go 1.22.1 diff --git a/sample/sample.go b/sample/sample.go index 4ceb225..b82ca64 100644 --- a/sample/sample.go +++ b/sample/sample.go @@ -1,55 +1,53 @@ -package squiggle +package sample -import "fmt" 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 +type Src = *rand.Rand +type func64 = func(Src) float64 -func Sample_unit_uniform(r src) float64 { +func Sample_unit_uniform(r Src) float64 { return r.Float64() } -func Sample_unit_normal(r src) 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_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_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_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 { +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) - + return Sample_normal(mean, std, r) } -func Sample_to(low float64, high float64, r src) float64 { +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 + // 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)) + return math.Exp(Sample_normal_from_90_ci(loglow, loghigh, r)) } -func Sample_mixture(fs []func64, weights []float64, r src) float64 { +func Sample_mixture(fs []func64, weights []float64, r Src) float64 { // fmt.Println("weights initially: ", weights) var sum_weights float64 = 0 @@ -104,6 +102,7 @@ func Sample_parallel(f func64, n_samples int) []float64 { return xs } + /* func main() { @@ -112,15 +111,15 @@ func main() { 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} + 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) } + model := func(r Src) float64 { return Sample_mixture(fs[0:], ws[0:], r) } n_samples := 1_000_000 - xs := sample_parallel(model, n_samples) + xs := Sample_parallel(model, n_samples) var avg float64 = 0 for _, x := range xs { avg += x @@ -133,10 +132,9 @@ func main() { 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 += Sample_mixture(fs[0:], ws[0:], r) } avg = avg / float64(n_samples) fmt.Printf("Average: %v\n", avg) - */ } */