more feng shui

README.md

@@ -1,6 +1,6 @@
-# A minimalist calculator for fermi estimation
+# A minimalist calculator for f estimation
 
-This project is a minimalist, stack-based DSL for Fermi estimation. It can multiply and divide scalars, lognormals and beta distributions.
+This project is a minimalist, stack-based DSL for f estimation. It can multiply and divide scalars, lognormals and beta distributions.
 
 ## Motivation
 
@@ -11,7 +11,7 @@ Sometimes, [Squiggle](https://github.com/quantified-uncertainty/squiggle), [simp
 Here is an example
 
 ```
-$ go run fermi.go
+$ go run f.go
 5000000 12000000
 => 5.0M 12.0M
 * beta 1 200
@@ -31,7 +31,7 @@ $ go run fermi.go
 Perhaps this example is more understandable with comments and better units:
 
 ```
-$ sed -u "s|#.*||" | sed -u 's|M|000000|g' | go run fermi.go
+$ sed -u "s|#.*||" | sed -u 's|M|000000|g' | go run f.go
 5M 12M # number of people living in Chicago
 => 5.0M 12.0M
 * beta 1 200 # fraction of people that have a piano
@@ -74,12 +74,40 @@ x
 
 The difference between `=: x` and `=. y` is that `=.` clears the stack after the assignment.
 
+If you type "help", you can see a small grammar:
+
+```
+help
+   Operation | Variable assignment | Special
+      Operation:                             operator operand
+            operator:                        (empty) | * | / | + | -
+            operand:                         scalar | lognormal | beta | variable
+              lognormal:                     low high
+              beta:                          beta alpha beta
+      Variable assignment:                   =: variable_name
+      Clear stack:                           . | c | clear
+      Variable assignment and clear stack:   =. variable_name
+      Other special operations:              help | debug | exit
+    Examples: 
+      + 2
+      / 2.5
+      * 1 10 (interpreted as lognormal)
+      + 1 10
+      * beta 1 10
+      1 10 (multiplication taken as default operation)
+        =: x
+      .
+      1 100
+      + x
+      exit
+```
+
 ## Installation
 
 ```
 make build
 sudo make install
-f # rather than the previous go run fermi.go
+f # rather than the previous go run f.go
 ```
 
 Why use make instead of the built-in go commands? Because the point of make is to be able to share command-line recipes. 
@@ -93,8 +121,8 @@ sed -u "s|#.*||" | sed -u 's|M|000000|g' | f
 cat more/piano-tuners.f | f
 cat more/piano-tuners-commented.f | sed -u "s|#.*||" | sed -u 's|M|000000|g' | f
 
-tee -a input.log | go run fermi.go | tee -a output.log
-tee -a io.log | go run fermi.go | tee -a io.log
+tee -a input.log | go run f.go | tee -a output.log
+tee -a io.log | go run f.go | tee -a io.log
 
 function f(){
   sed -u "s|#.*||" | 
@@ -114,8 +142,17 @@ Note that these sed commands are just hacks, and won't parse e.g., `3.5K` correc
 - Sums and divisions now also supported
 - For things between 0 and 1, consider using a beta distribution
 
+## Different levels of complexity
+
+The top level f.go file (400 lines) has a bunch of complexity: variables, parenthesis, samples, beta distributions. In the simple/ folder:
+
+- f_simple.go (370 lines) strips variables and parenthesis, but keeps beta distributions, samples, and addition and substraction
+- f_minimal.go (140 lines) strips everything that isn't lognormal and scalar multiplication and addition, plus a few debug options.
+
 ## Roadmap 
 
+Done: 
+
 - [x] Write README
 - [x] Add division?
 - [x] Read from file?
@@ -124,13 +161,10 @@ Note that these sed commands are just hacks, and won't parse e.g., `3.5K` correc
   - [x] Use a sed filter? 
 - [x] Add show more info version
 - [x] Scalar multiplication and division
-- [ ] Program into a small device, like a calculator?
-- [-] Think of some way of calling bc
 - [x] Think how to integrate with squiggle.c to draw samples
   - [x] Copy the time to botec go code
   - [x] Define samplers
   - [x] Call those samplers when operating on distributions that can't be operted on algebraically
-- [ ] Think about how to draw a histogram from samples
 - [x] Display output more nicely, with K/M/B/T
 - [x] Consider the following: make this into a stack-based DSL, with:
   - [x] Variables that can be saved to and then displayed
@@ -140,8 +174,20 @@ Note that these sed commands are just hacks, and won't parse e.g., `3.5K` correc
   - Joint types
   - Enums
 - [x] Fix correlation problem, by spinning up a new randomness thing every time some serial computation is done.
+- [x] Clean up error code. Right now only needed for division
+- [x] Maintain *both* a more complex thing that's more featureful *and* the more simple multiplication of lognormals thing.
+
+To (possibly) do:
+
+- [ ] Document parenthesis syntax
+- [ ] Allow input with K/M/T
+- [ ] Add functions. Now easier to do with an explicit representation of the stakc
+- [ ] Think about how to draw a histogram from samples
 - [ ] Dump samples to file
 - [ ] Represent samples/statistics in some other way
 - [ ] Perhaps use qsort rather than full sorting
+- [ ] Program into a small device, like a calculator?
 
-Some possible syntax for a more expressive stack-based DSL (now implemented)
+Discarded: 
+
+- [ ] ~~Think of some way of calling bc~~

f.go

@@ -4,8 +4,8 @@ import (
 	"bufio"
 	"errors"
 	"fmt"
-	"git.nunosempere.com/NunoSempere/fermi/sample"
 	"git.nunosempere.com/NunoSempere/fermi/pretty"
+	"git.nunosempere.com/NunoSempere/fermi/sample"
 	"math"
 	"os"
 	"sort"
@@ -14,7 +14,6 @@ import (
 )
 
 /* Types and interfaces */
-
 type Stack struct {
 	old_dist Dist
 	vars     map[string]Dist
@@ -52,14 +51,12 @@ func (p Scalar) Samples() []float64 {
 
 func (ln Lognormal) Samples() []float64 {
 	sampler := func(r sample.Src) float64 { return sample.Sample_to(ln.low, ln.high, r) }
-	// return sample.Sample_parallel(sampler, N_SAMPLES)
 	// Can't do parallel because then I'd have to await throughout the code
 	return sample.Sample_serially(sampler, N_SAMPLES)
 }
 
 func (beta Beta) Samples() []float64 {
 	sampler := func(r sample.Src) float64 { return sample.Sample_beta(beta.a, beta.b, r) }
-	// return sample.Sample_parallel(sampler, N_SAMPLES)
 	return sample.Sample_serially(sampler, N_SAMPLES)
 }
 
@@ -68,7 +65,7 @@ func (fs FilledSamples) Samples() []float64 {
 }
 
 /* Constants */
-const GENERAL_ERR_MSG = "  Operation | Variable assignment | Special\n" +
+const HELP_MSG = "  Operation | Variable assignment | Special\n" +
 	"    Operation:                             operator operand\n" +
 	"          operator:                        (empty) | * | / | + | -\n" +
 	"          operand:                         scalar | lognormal | beta | variable\n" +
@@ -94,131 +91,55 @@ const NORMAL90CONFIDENCE = 1.6448536269514727
 const INIT_DIST Scalar = Scalar(1)
 const N_SAMPLES = 100_000
 
-/* Pretty print for distributions */
-// Needs types
+/* Printers */
 func prettyPrintDist(dist Dist) {
 	switch v := dist.(type) {
 	case Lognormal:
 		fmt.Printf("=> ")
 		pretty.PrettyPrint2Floats(v.low, v.high)
-	case FilledSamples:
-		tmp_xs := make([]float64, N_SAMPLES)
-		copy(tmp_xs, v.xs)
-		sort.Slice(tmp_xs, func(i, j int) bool {
-			return tmp_xs[i] < tmp_xs[j]
-		})
-		low_int := N_SAMPLES / 20
-		low := tmp_xs[low_int]
-		high_int := N_SAMPLES * 19 / 20
-		high := tmp_xs[high_int]
-		fmt.Printf("=> ")
-		pretty.PrettyPrintFloat(low)
-		fmt.Printf(" ")
-		pretty.PrettyPrintFloat(high)
-		fmt.Printf(" (")
-		pretty.PrettyPrintInt(N_SAMPLES)
-		fmt.Printf(" samples)\n")
+		fmt.Println()
 	case Beta:
 		fmt.Printf("=> beta ")
 		pretty.PrettyPrint2Floats(v.a, v.b)
+		fmt.Println()
 	case Scalar:
 		fmt.Printf("=> scalar ")
 		w := float64(v)
 		pretty.PrettyPrintFloat(w)
 		fmt.Println()
+	case FilledSamples:
+		sorted_xs := make([]float64, N_SAMPLES)
+		copy(sorted_xs, v.xs)
+		sort.Slice(sorted_xs, func(i, j int) bool {
+			return sorted_xs[i] < sorted_xs[j]
+		})
+
+		low := sorted_xs[N_SAMPLES/20]
+		high := sorted_xs[N_SAMPLES*19/20]
+		fmt.Printf("=> ")
+		pretty.PrettyPrint2Floats(low, high)
+
+		fmt.Printf(" (")
+		pretty.PrettyPrintInt(N_SAMPLES)
+		fmt.Printf(" samples)")
+		fmt.Println()
 	default:
-		fmt.Printf("%v", v)
+		fmt.Printf("%v\n", v)
 	}
 }
 
-// Parse line into Distribution
-func parseLineErr(err_msg string) (string, Dist, error) {
-	fmt.Println(GENERAL_ERR_MSG)
+func printAndReturnErr(err_msg string) error {
 	fmt.Println(err_msg)
-	var errorDist Dist
-	return "", errorDist, errors.New(err_msg)
-}
-
-func parseLineIntoOpAndDist(line string, vars map[string]Dist) (string, Dist, error) {
-
-	words := strings.Split(strings.TrimSpace(line), " ")
-	op := ""
-	var dist Dist
-
-	switch words[0] {
-	case "*", "/", "+", "-":
-		op = words[0]
-		words = words[1:]
-	default:
-		op = "*" // later, change the below to
-	}
-
-	switch len(words) {
-	case 0:
-		return parseLineErr("Operator must have operand; can't operate on nothing")
-	case 1:
-		var_word, var_word_exists := vars[words[0]]
-		single_float, err1 := strconv.ParseFloat(words[0], 64) // abstract this away to search for K/M/B/T/etc.
-		switch {
-		case var_word_exists:
-			dist = var_word
-		case err1 == nil:
-			dist = Scalar(single_float)
-		case err1 != nil && !var_word_exists:
-			return parseLineErr("Trying to operate on a scalar, but scalar is neither a float nor an assigned variable")
-		}
-	case 2:
-		new_low, err1 := strconv.ParseFloat(words[0], 64)
-		new_high, err2 := strconv.ParseFloat(words[1], 64)
-		if err1 != nil || err2 != nil {
-			return parseLineErr("Trying to operate by a distribution, but distribution is not specified as two floats")
-		}
-		dist = Lognormal{low: new_low, high: new_high}
-	case 3:
-		if words[0] == "beta" || words[0] == "b" {
-			a, err1 := strconv.ParseFloat(words[1], 64)
-			b, err2 := strconv.ParseFloat(words[2], 64)
-			if err1 != nil || err2 != nil {
-				return parseLineErr("Trying to specify a beta distribution? Try beta 1 2")
-			}
-			dist = Beta{a: a, b: b}
-		} else {
-			return parseLineErr("Input not understood or not implemented yet")
-		}
-	default:
-		return parseLineErr("Input not understood or not implemented yet")
-	}
-	return op, dist, nil
-
-}
-
-func multiplyLogDists(l1 Lognormal, l2 Lognormal) Lognormal {
-	logmean1 := (math.Log(l1.high) + math.Log(l1.low)) / 2.0
-	logstd1 := (math.Log(l1.high) - math.Log(l1.low)) / (2.0 * NORMAL90CONFIDENCE)
-
-	logmean2 := (math.Log(l2.high) + math.Log(l2.low)) / 2.0
-	logstd2 := (math.Log(l2.high) - math.Log(l2.low)) / (2.0 * NORMAL90CONFIDENCE)
-
-	logmean_product := logmean1 + logmean2
-	logstd_product := math.Sqrt(logstd1*logstd1 + logstd2*logstd2)
-
-	h := logstd_product * NORMAL90CONFIDENCE
-	loglow := logmean_product - h
-	loghigh := logmean_product + h
-	return Lognormal{low: math.Exp(loglow), high: math.Exp(loghigh)}
-
-}
-
-func multiplyBetaDists(beta1 Beta, beta2 Beta) Beta {
-	return Beta{a: beta1.a + beta2.a, b: beta1.b + beta2.b}
+	fmt.Println(HELP_MSG)
+	return errors.New(err_msg)
 }
 
+/* Operations */
+// Generic operations with samples
 func operateDistsAsSamples(dist1 Dist, dist2 Dist, op string) (Dist, error) {
 
 	xs := dist1.Samples()
 	ys := dist2.Samples()
-	// fmt.Printf("xs: %v\n", xs)
-	// fmt.Printf("ys: %v\n", ys)
 	zs := make([]float64, N_SAMPLES)
 
 	for i := 0; i < N_SAMPLES; i++ {
@@ -239,10 +160,31 @@ func operateDistsAsSamples(dist1 Dist, dist2 Dist, op string) (Dist, error) {
 		}
 	}
 
-	// fmt.Printf("%v\n", zs)
 	return FilledSamples{xs: zs}, nil
 }
 
+// Multiplication
+func multiplyLogDists(l1 Lognormal, l2 Lognormal) Lognormal {
+	logmean1 := (math.Log(l1.high) + math.Log(l1.low)) / 2.0
+	logstd1 := (math.Log(l1.high) - math.Log(l1.low)) / (2.0 * NORMAL90CONFIDENCE)
+
+	logmean2 := (math.Log(l2.high) + math.Log(l2.low)) / 2.0
+	logstd2 := (math.Log(l2.high) - math.Log(l2.low)) / (2.0 * NORMAL90CONFIDENCE)
+
+	logmean_product := logmean1 + logmean2
+	logstd_product := math.Sqrt(logstd1*logstd1 + logstd2*logstd2)
+
+	h := logstd_product * NORMAL90CONFIDENCE
+	loglow := logmean_product - h
+	loghigh := logmean_product + h
+	return Lognormal{low: math.Exp(loglow), high: math.Exp(loghigh)}
+
+}
+
+func multiplyBetaDists(beta1 Beta, beta2 Beta) Beta {
+	return Beta{a: beta1.a + beta2.a, b: beta1.b + beta2.b}
+}
+
 func multiplyDists(old_dist Dist, new_dist Dist) (Dist, error) {
 
 	switch o := old_dist.(type) {
@@ -290,10 +232,16 @@ func divideDists(old_dist Dist, new_dist Dist) (Dist, error) {
 		{
 			switch n := new_dist.(type) {
 			case Lognormal:
-				// to do: check division by zero
+				if n.high == 0 || n.low == 0 {
+					fmt.Println("Error: Can't divide by 0.0")
+					return nil, errors.New("Error: division by zero")
+				}
 				return multiplyLogDists(o, Lognormal{low: 1.0 / n.high, high: 1.0 / n.low}), nil
 			case Scalar:
-				// to do: check division by zero
+				if n == 0.0 {
+					fmt.Println("Error: Can't divide by 0.0")
+					return nil, errors.New("Error: division by zero scalar")
+				}
 				return multiplyLogDists(o, Lognormal{low: 1.0 / float64(n), high: 1.0 / float64(n)}), nil
 			default:
 				return operateDistsAsSamples(old_dist, new_dist, "/")
@@ -305,7 +253,10 @@ func divideDists(old_dist Dist, new_dist Dist) (Dist, error) {
 			case Lognormal:
 				return multiplyLogDists(Lognormal{low: float64(o), high: float64(o)}, Lognormal{low: 1.0 / n.high, high: 1.0 / n.low}), nil
 			case Scalar:
-				// to do: check division by zero
+				if n == 0.0 {
+					fmt.Println("Error: Can't divide by 0.0")
+					return nil, errors.New("Error: division by zero scalar")
+				}
 				return Scalar(float64(o) / float64(n)), nil
 			default:
 				return operateDistsAsSamples(old_dist, new_dist, "/")
@@ -316,95 +267,133 @@ func divideDists(old_dist Dist, new_dist Dist) (Dist, error) {
 	}
 }
 
-/* Combine old dist and new line */
-// We want this as a function to be able to have parenthesis/recusion, possibly functions
-func operateStackWithDist(stack Stack, new_dist Dist, op string) Stack {
-
-	var combined_dist Dist
-	var err error
+// Generic distribution operations
+func operateDists(old_dist Dist, new_dist Dist, op string) (Dist, error) {
 	switch op {
 	case "*":
-		if combined_dist, err = multiplyDists(stack.old_dist, new_dist); err == nil {
-			stack.old_dist = combined_dist
-		}
+		return multiplyDists(old_dist, new_dist)
 	case "/":
-		if combined_dist, err = divideDists(stack.old_dist, new_dist); err == nil {
-			stack.old_dist = combined_dist
-		}
+		return divideDists(old_dist, new_dist)
 	case "+":
-		if combined_dist, err = operateDistsAsSamples(stack.old_dist, new_dist, "+"); err == nil {
-			stack.old_dist = combined_dist
-		}
+		return operateDistsAsSamples(old_dist, new_dist, "+")
 	case "-":
-		if combined_dist, err = operateDistsAsSamples(stack.old_dist, new_dist, "-"); err == nil {
-			stack.old_dist = combined_dist
-		}
+		return operateDistsAsSamples(old_dist, new_dist, "-")
 	default:
-		fmt.Println("Can't combine distributions in this way")
+		return nil, printAndReturnErr("Can't combine distributions in this way")
 	}
-	return stack
-
 }
 
+/* Parser and repl */
+func parseWordsErr(err_msg string) (string, Dist, error) {
+	return "", nil, printAndReturnErr(err_msg)
+}
+func parseWordsIntoOpAndDist(words []string, vars map[string]Dist) (string, Dist, error) {
+
+	op := ""
+	var dist Dist
+
+	switch words[0] {
+	case "*", "/", "+", "-":
+		op = words[0]
+		words = words[1:]
+	default:
+		op = "*" // later, change the below to
+	}
+
+	switch len(words) {
+	case 0:
+		return parseWordsErr("Operator must have operand; can't operate on nothing")
+	case 1:
+		var_word, var_word_exists := vars[words[0]]
+		single_float, err1 := strconv.ParseFloat(words[0], 64) // abstract this away to search for K/M/B/T/etc.
+		switch {
+		case var_word_exists:
+			dist = var_word
+		case err1 == nil:
+			dist = Scalar(single_float)
+		case err1 != nil && !var_word_exists:
+			return parseWordsErr("Trying to operate on a scalar, but scalar is neither a float nor an assigned variable")
+		}
+	case 2:
+		new_low, err1 := strconv.ParseFloat(words[0], 64)
+		new_high, err2 := strconv.ParseFloat(words[1], 64)
+		if err1 != nil || err2 != nil {
+			return parseWordsErr("Trying to operate by a distribution, but distribution is not specified as two floats")
+		}
+		dist = Lognormal{low: new_low, high: new_high}
+	case 3:
+		if words[0] == "beta" || words[0] == "b" {
+			a, err1 := strconv.ParseFloat(words[1], 64)
+			b, err2 := strconv.ParseFloat(words[2], 64)
+			if err1 != nil || err2 != nil {
+				return parseWordsErr("Trying to specify a beta distribution? Try beta 1 2")
+			}
+			dist = Beta{a: a, b: b}
+		} else {
+			return parseWordsErr("Input not understood or not implemented yet")
+		}
+	default:
+		return parseWordsErr("Input not understood or not implemented yet")
+	}
+	return op, dist, nil
+}
+
+/* Combine old dist and new line */
+// We want this as a function (rather than just be in main)
+// to be able to have parenthesis/recusion, possibly functions
+
 func runRepl(stack Stack, reader *bufio.Reader) Stack {
-
-
-	replForLoop:
+replForLoop:
 	for {
 		new_line, _ := reader.ReadString('\n')
 		words := strings.Split(strings.TrimSpace(new_line), " ")
 		switch {
-		/* Empty line case */
-		case strings.TrimSpace(new_line) == "":
-			continue replForLoop
+		case strings.TrimSpace(new_line) == "": /* Empty line case */
 		/* Parenthesis */
 		case len(words) == 2 && (words[0] == "*" || words[0] == "+" || words[0] == "-" || words[0] == "/") && words[1] == "(":
-			new_stack := runRepl(Stack{old_dist: INIT_DIST, vars: stack.vars}, reader) // to do: think if I want to shadow variables or not (right now, variables persist, since I'm copying the map(.
-			stack = operateStackWithDist(stack, new_stack.old_dist, words[0])
-			prettyPrintDist(stack.old_dist)
+			new_stack := runRepl(Stack{old_dist: INIT_DIST, vars: stack.vars}, reader)
+			combined_dist, err := operateDists(stack.old_dist, new_stack.old_dist, words[0])
+			if err == nil {
+				stack.old_dist = combined_dist
+			}
 		case len(words) == 1 && words[0] == ")":
 			return stack
 		/* Special operations */
 		case words[0] == "exit" || words[0] == "e":
 			os.Exit(0)
 		case words[0] == "help" || words[0] == "h":
-			fmt.Println(GENERAL_ERR_MSG)
-			continue replForLoop
+			fmt.Println(HELP_MSG)
 		case words[0] == "debug" || words[0] == "d":
-			fmt.Printf("Old dist: %v\n", stack.old_dist)
-			fmt.Printf("Vars: %v\n", stack.vars)
-			continue replForLoop
+			fmt.Printf("%v", stack)
 		case words[0] == "clear" || words[0] == "c" || words[0] == ".":
 			stack.old_dist = INIT_DIST
 			fmt.Println()
-			continue replForLoop
-		/* Variable assignment */ 
+		/* Variable assignment */
 		case words[0] == "=:" && len(words) == 2:
 			stack.vars[words[1]] = stack.old_dist
 			fmt.Printf("%s ", words[1])
-			prettyPrintDist(stack.old_dist)
-			continue replForLoop
 		case words[0] == "=." && len(words) == 2:
 			stack.vars[words[1]] = stack.old_dist
 			fmt.Printf("%s ", words[1])
 			prettyPrintDist(stack.old_dist)
 			stack.old_dist = INIT_DIST
-			fmt.Println()
-			continue replForLoop
+			// fmt.Println()
+			// continue replForLoop
 		default:
-			op, new_dist, err := parseLineIntoOpAndDist(new_line, stack.vars)
+			op, new_dist, err := parseWordsIntoOpAndDist(words, stack.vars)
 			if err != nil {
 				continue replForLoop
 			}
-			stack = operateStackWithDist(stack, new_dist, op)
-			prettyPrintDist(stack.old_dist)
+			combined_dist, err := operateDists(stack.old_dist, new_dist, op)
+			if err == nil {
+				stack.old_dist = combined_dist
+			}
 		}
+		prettyPrintDist(stack.old_dist)
 	}
 }
 
-/* Main event loop */
 func main() {
-
 	reader := bufio.NewReader(os.Stdin)
 	stack := Stack{old_dist: INIT_DIST, vars: make(map[string]Dist)}
 	runRepl(stack, reader)

makefile

@@ -1,9 +1,9 @@
 run:
-	go run fermi.go
+	go run f.go
 
 build:
-	go build fermi.go
+	go build f.go
 
-install: fermi
-	sudo mv fermi /usr/bin/fermi
-	sudo ln -s /usr/bin/fermi /usr/bin/f
+install: f
+	sudo mv f /usr/bin/f
+	sudo ln -s /usr/bin/f /usr/bin/fermi

pretty/pretty.go

@@ -48,5 +48,4 @@ func PrettyPrint2Floats(low float64, high float64) {
 	PrettyPrintFloat(low)
 	fmt.Printf(" ")
 	PrettyPrintFloat(high)
-	fmt.Printf("\n")
 }

simple/f_minimal.go

@@ -0,0 +1,140 @@
+package main
+
+import (
+	"bufio"
+	"errors"
+	"fmt"
+	"git.nunosempere.com/NunoSempere/fermi/pretty"
+	"math"
+	"os"
+	"strconv"
+	"strings"
+)
+
+type Lognormal struct {
+	low  float64
+	high float64
+}
+
+const HELP_MSG = "  Operation | Variable assignment | Special\n" +
+	"    Operation:                             operator operand\n" +
+	"          operator:                        (empty) | * | / | + | -\n" +
+	"          operand:                         scalar | lognormal\n" +
+	"            lognormal:                     low high\n" +
+	"    Clear stack:                           . | c | clear\n" +
+	"    Other special operations:              help | debug | exit\n" +
+	"  Examples: \n" +
+	"    / 2.5\n" +
+	"    * 1 10 (interpreted as lognormal)\n" +
+	"    / 1 10\n" +
+	"    1 10 (multiplication taken as default operation)\n" +
+	"    .\n" +
+	"    exit\n"
+const NORMAL90CONFIDENCE = 1.6448536269514727
+const N_SAMPLES = 100_000
+
+func prettyPrintLognormal(l Lognormal) {
+	fmt.Printf("=> ")
+	pretty.PrettyPrint2Floats(l.low, l.high)
+	fmt.Println()
+}
+
+func printAndReturnErr(err_msg string) error {
+	fmt.Println(err_msg)
+	fmt.Println(HELP_MSG)
+	return errors.New(err_msg)
+}
+
+func multiplyLogDists(l1 Lognormal, l2 Lognormal) Lognormal {
+	logmean1 := (math.Log(l1.high) + math.Log(l1.low)) / 2.0
+	logstd1 := (math.Log(l1.high) - math.Log(l1.low)) / (2.0 * NORMAL90CONFIDENCE)
+
+	logmean2 := (math.Log(l2.high) + math.Log(l2.low)) / 2.0
+	logstd2 := (math.Log(l2.high) - math.Log(l2.low)) / (2.0 * NORMAL90CONFIDENCE)
+
+	logmean_product := logmean1 + logmean2
+	logstd_product := math.Sqrt(logstd1*logstd1 + logstd2*logstd2)
+
+	h := logstd_product * NORMAL90CONFIDENCE
+	loglow := logmean_product - h
+	loghigh := logmean_product + h
+	return Lognormal{low: math.Exp(loglow), high: math.Exp(loghigh)}
+}
+
+func divideLogDists(l1 Lognormal, l2 Lognormal) (Lognormal, error) {
+	if l2.high == 0 || l2.low == 0 {
+		fmt.Println("Error: Can't divide by 0.0")
+		return Lognormal{}, errors.New("Error: division by zero")
+	}
+	return multiplyLogDists(l1, Lognormal{low: 1.0 / l2.high, high: 1.0 / l2.low}), nil
+}
+
+func parseWordsErr(err_msg string) (string, Lognormal, error) {
+	return "", Lognormal{}, printAndReturnErr(err_msg)
+}
+func parseWordsIntoOpAndDist(words []string) (string, Lognormal, error) {
+	op := ""
+	var dist Lognormal
+
+	switch words[0] {
+	case "*", "/":
+		op = words[0]
+		words = words[1:]
+	default:
+		op = "*"
+	}
+
+	switch len(words) {
+	case 0:
+		return parseWordsErr("Operator must have operand; can't operate on nothing")
+	case 1:
+		single_float, err := strconv.ParseFloat(words[0], 64) // abstract this away to search for K/M/B/T/etc.
+		if err != nil {
+			return parseWordsErr("Trying to operate on a scalar, but scalar is neither a float nor an assigned variable")
+		}
+		dist = Lognormal{low: single_float, high: single_float}
+	case 2:
+		new_low, err1 := strconv.ParseFloat(words[0], 64)
+		new_high, err2 := strconv.ParseFloat(words[1], 64)
+		if err1 != nil || err2 != nil {
+			return parseWordsErr("Trying to operate by a distribution, but distribution is not specified as two floats")
+		}
+		dist = Lognormal{low: new_low, high: new_high}
+	default:
+		return parseWordsErr("Input not understood or not implemented yet")
+	}
+	return op, dist, nil
+}
+
+func main() {
+	reader := bufio.NewReader(os.Stdin)
+	old_dist := Lognormal{low: 1, high: 1}
+
+replForLoop:
+	for {
+
+		new_line, _ := reader.ReadString('\n')
+		words := strings.Split(strings.TrimSpace(new_line), " ")
+
+		if strings.TrimSpace(new_line) == "" {
+			continue replForLoop
+		}
+
+		op, new_dist, err := parseWordsIntoOpAndDist(words)
+		if err != nil {
+			continue replForLoop
+		}
+		switch op {
+		case "*":
+			old_dist = multiplyLogDists(old_dist, new_dist)
+		case "/":
+			result_dist, err := divideLogDists(old_dist, new_dist)
+			if err != nil {
+				continue replForLoop
+			}
+			old_dist = result_dist
+
+		}
+		prettyPrintLognormal(old_dist)
+	}
+}

simple/f_simple.go

@@ -0,0 +1,370 @@
+package main
+
+import (
+	"bufio"
+	"errors"
+	"fmt"
+	"git.nunosempere.com/NunoSempere/fermi/pretty"
+	"git.nunosempere.com/NunoSempere/fermi/sample"
+	"math"
+	"os"
+	"sort"
+	"strconv"
+	"strings"
+)
+
+/* Types and interfaces */
+
+type Dist interface {
+	Samples() []float64
+}
+
+type Scalar float64
+
+type Lognormal struct {
+	low  float64
+	high float64
+}
+
+type Beta struct {
+	a float64
+	b float64
+}
+
+type FilledSamples struct {
+	xs []float64
+}
+
+/* Dist interface functions */
+// https://go.dev/tour/methods/9
+func (p Scalar) Samples() []float64 {
+	xs := make([]float64, N_SAMPLES)
+	for i := 0; i < N_SAMPLES; i++ {
+		xs[i] = float64(p)
+	}
+	return xs
+}
+
+func (ln Lognormal) Samples() []float64 {
+	sampler := func(r sample.Src) float64 { return sample.Sample_to(ln.low, ln.high, r) }
+	return sample.Sample_serially(sampler, N_SAMPLES)
+}
+
+func (beta Beta) Samples() []float64 {
+	sampler := func(r sample.Src) float64 { return sample.Sample_beta(beta.a, beta.b, r) }
+	return sample.Sample_serially(sampler, N_SAMPLES)
+}
+
+func (fs FilledSamples) Samples() []float64 {
+	return fs.xs
+}
+
+/* Constants */
+const HELP_MSG = "  Operation | Variable assignment | Special\n" +
+	"    Operation:                             operator operand\n" +
+	"          operator:                        (empty) | * | / | + | -\n" +
+	"          operand:                         scalar | lognormal | beta\n" +
+	"            lognormal:                     low high\n" +
+	"            beta:                          beta alpha beta\n" +
+	"    Clear stack:                           . | c | clear\n" +
+	"    Other special operations:              help | debug | exit\n" +
+	"  Examples: \n" +
+	"    + 2\n" +
+	"    / 2.5\n" +
+	"    * 1 10 (interpreted as lognormal)\n" +
+	"    + 1 10\n" +
+	"    * beta 1 10\n" +
+	"    1 10 (multiplication taken as default operation)\n" +
+	"    .\n" +
+	"    1 100\n" +
+	"    exit\n"
+const NORMAL90CONFIDENCE = 1.6448536269514727
+const INIT_DIST Scalar = Scalar(1)
+const N_SAMPLES = 100_000
+
+/* Printers */
+func prettyPrintDist(dist Dist) {
+	switch v := dist.(type) {
+	case Lognormal:
+		fmt.Printf("=> ")
+		pretty.PrettyPrint2Floats(v.low, v.high)
+		fmt.Println()
+	case Beta:
+		fmt.Printf("=> beta ")
+		pretty.PrettyPrint2Floats(v.a, v.b)
+		fmt.Println()
+	case Scalar:
+		fmt.Printf("=> scalar ")
+		w := float64(v)
+		pretty.PrettyPrintFloat(w)
+		fmt.Println()
+	case FilledSamples:
+		sorted_xs := make([]float64, N_SAMPLES)
+		copy(sorted_xs, v.xs)
+		sort.Slice(sorted_xs, func(i, j int) bool {
+			return sorted_xs[i] < sorted_xs[j]
+		})
+
+		low := sorted_xs[N_SAMPLES/20]
+		high := sorted_xs[N_SAMPLES*19/20]
+		fmt.Printf("=> ")
+		pretty.PrettyPrint2Floats(low, high)
+
+		fmt.Printf(" (")
+		pretty.PrettyPrintInt(N_SAMPLES)
+		fmt.Printf(" samples)")
+		fmt.Println()
+	default:
+		fmt.Printf("%v\n", v)
+	}
+}
+
+func printAndReturnErr(err_msg string) error {
+	fmt.Println(err_msg)
+	fmt.Println(HELP_MSG)
+	return errors.New(err_msg)
+}
+
+/* Operations */
+// Generic operations with samples
+func operateDistsAsSamples(dist1 Dist, dist2 Dist, op string) (Dist, error) {
+
+	xs := dist1.Samples()
+	ys := dist2.Samples()
+	zs := make([]float64, N_SAMPLES)
+
+	for i := 0; i < N_SAMPLES; i++ {
+		switch op {
+		case "*":
+			zs[i] = xs[i] * ys[i]
+		case "/":
+			if ys[0] != 0 {
+				zs[i] = xs[i] / ys[i]
+			} else {
+				fmt.Println("Error: When dividing as samples, division by zero")
+				return nil, errors.New("Division by zero")
+			}
+		case "+":
+			zs[i] = xs[i] + ys[i]
+		case "-":
+			zs[i] = xs[i] - ys[i]
+		}
+	}
+
+	// fmt.Printf("%v\n", zs)
+	return FilledSamples{xs: zs}, nil
+}
+
+// Multiplication
+func multiplyLogDists(l1 Lognormal, l2 Lognormal) Lognormal {
+	logmean1 := (math.Log(l1.high) + math.Log(l1.low)) / 2.0
+	logstd1 := (math.Log(l1.high) - math.Log(l1.low)) / (2.0 * NORMAL90CONFIDENCE)
+
+	logmean2 := (math.Log(l2.high) + math.Log(l2.low)) / 2.0
+	logstd2 := (math.Log(l2.high) - math.Log(l2.low)) / (2.0 * NORMAL90CONFIDENCE)
+
+	logmean_product := logmean1 + logmean2
+	logstd_product := math.Sqrt(logstd1*logstd1 + logstd2*logstd2)
+
+	h := logstd_product * NORMAL90CONFIDENCE
+	loglow := logmean_product - h
+	loghigh := logmean_product + h
+	return Lognormal{low: math.Exp(loglow), high: math.Exp(loghigh)}
+
+}
+
+func multiplyBetaDists(beta1 Beta, beta2 Beta) Beta {
+	return Beta{a: beta1.a + beta2.a, b: beta1.b + beta2.b}
+}
+
+func multiplyDists(old_dist Dist, new_dist Dist) (Dist, error) {
+
+	switch o := old_dist.(type) {
+	case Lognormal:
+		{
+			switch n := new_dist.(type) {
+			case Lognormal:
+				return multiplyLogDists(o, n), nil
+			case Scalar:
+				return multiplyLogDists(o, Lognormal{low: float64(n), high: float64(n)}), nil
+			default:
+				return operateDistsAsSamples(old_dist, new_dist, "*")
+			}
+		}
+	case Scalar:
+		{
+			if o == 1 {
+				return new_dist, nil
+			}
+			switch n := new_dist.(type) {
+			case Lognormal:
+				return multiplyLogDists(Lognormal{low: float64(o), high: float64(o)}, n), nil
+			case Scalar:
+				return Scalar(float64(o) * float64(n)), nil
+			default:
+				return operateDistsAsSamples(old_dist, new_dist, "*")
+			}
+		}
+	case Beta:
+		switch n := new_dist.(type) {
+		case Beta:
+			return multiplyBetaDists(o, n), nil
+		default:
+			return operateDistsAsSamples(old_dist, new_dist, "*")
+		}
+	default:
+		return operateDistsAsSamples(old_dist, new_dist, "*")
+	}
+}
+
+func divideDists(old_dist Dist, new_dist Dist) (Dist, error) {
+
+	switch o := old_dist.(type) {
+	case Lognormal:
+		{
+			switch n := new_dist.(type) {
+			case Lognormal:
+				// to do: check division by zero
+				if n.high == 0 || n.low == 0 {
+					fmt.Println("Error: Can't divide by 0.0")
+					return nil, errors.New("Error: division by zero")
+				}
+				return multiplyLogDists(o, Lognormal{low: 1.0 / n.high, high: 1.0 / n.low}), nil
+			case Scalar:
+				// to do: check division by zero
+				if n == 0.0 {
+					fmt.Println("Error: Can't divide by 0.0")
+					return nil, errors.New("Error: division by zero scalar")
+				}
+				return multiplyLogDists(o, Lognormal{low: 1.0 / float64(n), high: 1.0 / float64(n)}), nil
+			default:
+				return operateDistsAsSamples(old_dist, new_dist, "/")
+			}
+		}
+	case Scalar:
+		{
+			switch n := new_dist.(type) {
+			case Lognormal:
+				return multiplyLogDists(Lognormal{low: float64(o), high: float64(o)}, Lognormal{low: 1.0 / n.high, high: 1.0 / n.low}), nil
+			case Scalar:
+				// to do: check division by zero
+				if n == 0.0 {
+					fmt.Println("Error: Can't divide by 0.0")
+					return nil, errors.New("Error: division by zero scalar")
+				}
+				return Scalar(float64(o) / float64(n)), nil
+			default:
+				return operateDistsAsSamples(old_dist, new_dist, "/")
+			}
+		}
+	default:
+		return operateDistsAsSamples(old_dist, new_dist, "/")
+	}
+}
+
+// Generic distribution operations
+func operateDists(old_dist Dist, new_dist Dist, op string) (Dist, error) {
+	switch op {
+	case "*":
+		return multiplyDists(old_dist, new_dist)
+	case "/":
+		return divideDists(old_dist, new_dist)
+	case "+":
+		return operateDistsAsSamples(old_dist, new_dist, "+")
+	case "-":
+		return operateDistsAsSamples(old_dist, new_dist, "-")
+	default:
+		return nil, printAndReturnErr("Can't combine distributions in this way")
+	}
+}
+
+/* Parser and repl */
+func parseWordsErr(err_msg string) (string, Dist, error) {
+	return "", nil, printAndReturnErr(err_msg)
+}
+func parseWordsIntoOpAndDist(words []string) (string, Dist, error) {
+
+	op := ""
+	var dist Dist
+
+	switch words[0] {
+	case "*", "/", "+", "-":
+		op = words[0]
+		words = words[1:]
+	default:
+		op = "*" // later, change the below to
+	}
+
+	switch len(words) {
+	case 0:
+		return parseWordsErr("Operator must have operand; can't operate on nothing")
+	case 1:
+		single_float, err := strconv.ParseFloat(words[0], 64) // abstract this away to search for K/M/B/T/etc.
+		if err != nil {
+			return parseWordsErr("Trying to operate on a scalar, but scalar is neither a float nor an assigned variable")
+		}
+		dist = Scalar(single_float)
+	case 2:
+		new_low, err1 := strconv.ParseFloat(words[0], 64)
+		new_high, err2 := strconv.ParseFloat(words[1], 64)
+		if err1 != nil || err2 != nil {
+			return parseWordsErr("Trying to operate by a distribution, but distribution is not specified as two floats")
+		}
+		dist = Lognormal{low: new_low, high: new_high}
+	case 3:
+		if words[0] == "beta" || words[0] == "b" {
+			a, err1 := strconv.ParseFloat(words[1], 64)
+			b, err2 := strconv.ParseFloat(words[2], 64)
+			if err1 != nil || err2 != nil {
+				return parseWordsErr("Trying to specify a beta distribution? Try beta 1 2")
+			}
+			dist = Beta{a: a, b: b}
+		} else {
+			return parseWordsErr("Input not understood or not implemented yet")
+		}
+	default:
+		return parseWordsErr("Input not understood or not implemented yet")
+	}
+	return op, dist, nil
+
+}
+
+func main() {
+	reader := bufio.NewReader(os.Stdin)
+	var old_dist Dist = INIT_DIST
+
+replForLoop:
+	for {
+
+		new_line, _ := reader.ReadString('\n')
+		words := strings.Split(strings.TrimSpace(new_line), " ")
+
+		switch {
+		case strings.TrimSpace(new_line) == "": /* Empty line case */
+			continue replForLoop
+		case words[0] == "exit" || words[0] == "e":
+			os.Exit(0)
+		case words[0] == "help" || words[0] == "h":
+			fmt.Println(HELP_MSG)
+			continue replForLoop
+		case words[0] == "debug" || words[0] == "d":
+			fmt.Printf("%v", old_dist)
+			continue replForLoop
+		case words[0] == "clear" || words[0] == "c" || words[0] == ".":
+			old_dist = INIT_DIST
+			fmt.Println()
+			continue replForLoop
+		}
+
+		op, new_dist, err := parseWordsIntoOpAndDist(words)
+		if err != nil {
+			continue replForLoop
+		}
+		combined_dist, err := operateDists(old_dist, new_dist, op)
+		if err != nil {
+			continue replForLoop
+		}
+		old_dist = combined_dist
+		prettyPrintDist(old_dist)
+	}
+}