Add a way to plot the output from generators

For visualization, add a simple script for generating scatter plots and
a binary (via examples) to plot the inputs given various domains.

crates/libm-test/examples/plot_domains.rs

@@ -0,0 +1,105 @@
+//! Program to write all inputs from a generator to a file, then invoke a Julia script to plot
+//! them. Output is in `target/plots`.
+//!
+//! Requires Julia with the `CairoMakie` dependency.
+//!
+//! Note that running in release mode by default generates a _lot_ more datapoints, which
+//! causes plotting to be extremely slow (some simplification to be done in the script).
+
+use std::fmt::Write as _;
+use std::io::{BufWriter, Write};
+use std::path::Path;
+use std::process::Command;
+use std::{env, fs};
+
+use libm_test::domain::HasDomain;
+use libm_test::gen::{domain_logspace, edge_cases};
+use libm_test::{MathOp, op};
+
+const JL_PLOT: &str = "examples/plot_file.jl";
+
+fn main() {
+    let manifest_env = env::var("CARGO_MANIFEST_DIR").unwrap();
+    let manifest_dir = Path::new(&manifest_env);
+    let out_dir = manifest_dir.join("../../target/plots");
+    if !out_dir.exists() {
+        fs::create_dir(&out_dir).unwrap();
+    }
+
+    let jl_script = manifest_dir.join(JL_PLOT);
+    let mut config = format!(r#"out_dir = "{}""#, out_dir.display());
+    config.write_str("\n\n").unwrap();
+
+    // Plot a few domains with some functions that use them.
+    plot_one_operator::<op::sqrtf::Routine>(&out_dir, &mut config);
+    plot_one_operator::<op::cosf::Routine>(&out_dir, &mut config);
+    plot_one_operator::<op::cbrtf::Routine>(&out_dir, &mut config);
+
+    let config_path = out_dir.join("config.toml");
+    fs::write(&config_path, config).unwrap();
+
+    // The script expects a path to `config.toml` to be passed as its only argument
+    let mut cmd = Command::new("julia");
+    if cfg!(optimizations_enabled) {
+        cmd.arg("-O3");
+    }
+    cmd.arg(jl_script).arg(config_path);
+
+    println!("launching script... {cmd:?}");
+    cmd.status().unwrap();
+}
+
+/// Run multiple generators for a single operator.
+fn plot_one_operator<Op>(out_dir: &Path, config: &mut String)
+where
+    Op: MathOp<FTy = f32> + HasDomain<f32>,
+{
+    plot_one_generator(
+        out_dir,
+        Op::BASE_NAME.as_str(),
+        "logspace",
+        config,
+        domain_logspace::get_test_cases::<Op>(),
+    );
+    plot_one_generator(
+        out_dir,
+        Op::BASE_NAME.as_str(),
+        "edge_cases",
+        config,
+        edge_cases::get_test_cases::<Op, _>(),
+    );
+}
+
+/// Plot the output of a single generator.
+fn plot_one_generator(
+    out_dir: &Path,
+    fn_name: &str,
+    gen_name: &str,
+    config: &mut String,
+    gen: impl Iterator<Item = (f32,)>,
+) {
+    let text_file = out_dir.join(format!("input-{fn_name}-{gen_name}.txt"));
+
+    let f = fs::File::create(&text_file).unwrap();
+    let mut w = BufWriter::new(f);
+    let mut count = 0u64;
+
+    for input in gen {
+        writeln!(w, "{:e}", input.0).unwrap();
+        count += 1;
+    }
+
+    w.flush().unwrap();
+    println!("generated {count} inputs for {fn_name}-{gen_name}");
+
+    writeln!(
+        config,
+        r#"[[input]]
+function = "{fn_name}"
+generator = "{gen_name}"
+input_file = "{}"
+"#,
+        text_file.to_str().unwrap()
+    )
+    .unwrap()
+}

crates/libm-test/examples/plot_file.jl

@@ -0,0 +1,157 @@
+"A quick script for plotting a list of floats.
+
+Takes a path to a TOML file (Julia has builtin TOML support but not JSON) which
+specifies a list of source files to plot. Plots are done with both a linear and
+a log scale.
+
+Requires [Makie] (specifically CairoMakie) for plotting.
+
+[Makie]: https://docs.makie.org/stable/
+"
+
+using CairoMakie
+using TOML
+
+function main()::Nothing
+    CairoMakie.activate!(px_per_unit=10)
+    config_path = ARGS[1]
+
+    cfg = Dict()
+    open(config_path, "r") do f
+        cfg = TOML.parse(f)
+    end
+
+    out_dir = cfg["out_dir"]
+    for input in cfg["input"]
+        fn_name = input["function"]
+        gen_name = input["generator"]
+        input_file = input["input_file"]
+
+        plot_one(input_file, out_dir, fn_name, gen_name)
+    end
+end
+
+"Read inputs from a file, create both linear and log plots for one function"
+function plot_one(
+    input_file::String,
+    out_dir::String,
+    fn_name::String,
+    gen_name::String,
+)::Nothing
+    fig = Figure()
+
+    lin_out_file = joinpath(out_dir, "plot-$fn_name-$gen_name.png")
+    log_out_file = joinpath(out_dir, "plot-$fn_name-$gen_name-log.png")
+
+    # Map string function names to callable functions
+    if fn_name == "cos"
+        orig_func = cos
+        xlims = (-6.0, 6.0)
+        xlims_log = (-pi * 10, pi * 10)
+    elseif fn_name == "cbrt"
+        orig_func = cbrt
+        xlims = (-2.0, 2.0)
+        xlims_log = (-1000.0, 1000.0)
+    elseif fn_name == "sqrt"
+        orig_func = sqrt
+        xlims = (-1.1, 6.0)
+        xlims_log = (-1.1, 5000.0)
+    else
+        println("unrecognized function name `$fn_name`; update plot_file.jl")
+        exit(1)
+    end
+
+    # Edge cases don't do much beyond +/-1, except for infinity.
+    if gen_name == "edge_cases"
+        xlims = (-1.1, 1.1)
+        xlims_log = (-1.1, 1.1)
+    end
+
+    # Turn domain errors into NaN
+    func(x) = map_or(x, orig_func, NaN)
+
+    # Parse a series of X values produced by the generator
+    inputs = readlines(input_file)
+    gen_x = map((v) -> parse(Float32, v), inputs)
+
+    do_plot(
+        fig, gen_x, func, xlims[1], xlims[2],
+        "$fn_name $gen_name (linear scale)",
+        lin_out_file, false,
+    )
+
+    do_plot(
+        fig, gen_x, func, xlims_log[1], xlims_log[2],
+        "$fn_name $gen_name (log scale)",
+        log_out_file, true,
+    )
+end
+
+"Create a single plot"
+function do_plot(
+    fig::Figure,
+    gen_x::Vector{F},
+    func::Function,
+    xmin::AbstractFloat,
+    xmax::AbstractFloat,
+    title::String,
+    out_file::String,
+    logscale::Bool,
+)::Nothing where F<:AbstractFloat
+    println("plotting $title")
+
+    # `gen_x` is the values the generator produces. `actual_x` is for plotting a
+    # continuous function.
+    input_min = xmin - 1.0
+    input_max = xmax + 1.0
+    gen_x = filter((v) -> v >= input_min && v <= input_max, gen_x)
+    markersize = length(gen_x) < 10_000 ? 6.0 : 4.0
+
+    steps = 10_000
+    if logscale
+        r = LinRange(symlog10(input_min), symlog10(input_max), steps)
+        actual_x = sympow10.(r)
+        xscale = Makie.pseudolog10
+    else
+        actual_x = LinRange(input_min, input_max, steps)
+        xscale = identity
+    end
+
+    gen_y = @. func(gen_x)
+    actual_y = @. func(actual_x)
+
+    ax = Axis(fig[1, 1], xscale=xscale, title=title)
+
+    lines!(
+        ax, actual_x, actual_y, color=(:lightblue, 0.6),
+        linewidth=6.0, label="true function",
+    )
+    scatter!(
+        ax, gen_x, gen_y, color=(:darkblue, 0.9),
+        markersize=markersize, label="checked inputs",
+    )
+    axislegend(ax, position=:rb, framevisible=false)
+
+    save(out_file, fig)
+    delete!(ax)
+end
+
+"Apply a function, returning the default if there is a domain error"
+function map_or(
+    input::AbstractFloat,
+    f::Function,
+    default::Any
+)::Union{AbstractFloat,Any}
+    try
+        return f(input)
+    catch
+        return default
+    end
+end
+
+# Operations for logarithms that are symmetric about 0
+C = 10
+symlog10(x::Number) = sign(x) * (log10(1 + abs(x)/(10^C)))
+sympow10(x::Number) = (10^C) * (10^x - 1)
+
+main()