go-hep.org/x/hep@v0.38.1/hbook/rand_h1d_example_test.go (about)

     1  // Copyright ©2020 The go-hep Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package hbook_test
     6  
     7  import (
     8  	"fmt"
     9  	"image/color"
    10  	"log"
    11  	"math/rand/v2"
    12  	"os"
    13  	"path"
    14  	"testing"
    15  
    16  	"go-hep.org/x/hep/hbook"
    17  	"go-hep.org/x/hep/hplot"
    18  	"gonum.org/v1/gonum/stat/distuv"
    19  	"gonum.org/v1/plot/cmpimg"
    20  	"gonum.org/v1/plot/vg"
    21  )
    22  
    23  func ExampleRand1D() {
    24  	const N = 10000
    25  	var (
    26  		h1  = hbook.NewH1D(100, -10, 10)
    27  		src = rand.New(rand.NewPCG(1234, 1234))
    28  		rnd = distuv.Normal{
    29  			Mu:    0,
    30  			Sigma: 2,
    31  			Src:   src,
    32  		}
    33  	)
    34  
    35  	for range N {
    36  		h1.Fill(rnd.Rand(), 1)
    37  	}
    38  
    39  	var (
    40  		h2 = hbook.NewH1D(100, -10, 10)
    41  		hr = hbook.NewRand1D(h1, rand.NewPCG(5678, 5678))
    42  	)
    43  
    44  	for range N {
    45  		h2.Fill(hr.Rand(), 1)
    46  	}
    47  
    48  	fmt.Printf(
    49  		"h1: mean=%+8f std-dev=%+8f +/- %8f\n",
    50  		h1.XMean(), h1.XStdDev(), h1.XStdErr(),
    51  	)
    52  	fmt.Printf(
    53  		"cdf(0)= %+1.1f\ncdf(1)= %+1.1f\n",
    54  		rnd.CDF(0), rnd.CDF(1),
    55  	)
    56  	fmt.Printf(
    57  		"h2: mean=%+8f std-dev=%+8f +/- %8f\n",
    58  		h2.XMean(), h2.XStdDev(), h2.XStdErr(),
    59  	)
    60  	fmt.Printf(
    61  		"cdf(0)= %+1.1f\ncdf(1)= %+1.1f\n",
    62  		hr.CDF(0), hr.CDF(1),
    63  	)
    64  
    65  	h1.Scale(1. / h1.Integral(h1.XMin(), h1.XMax()))
    66  	h2.Scale(1. / h2.Integral(h2.XMin(), h2.XMax()))
    67  
    68  	{
    69  		rp := hplot.NewRatioPlot()
    70  		rp.Ratio = 0.3
    71  
    72  		rp.Top.Title.Text = "Distributions"
    73  		rp.Top.Y.Label.Text = "Y"
    74  
    75  		hh1 := hplot.NewH1D(h1)
    76  		hh1.FillColor = color.NRGBA{R: 255, A: 100}
    77  
    78  		hh2 := hplot.NewH1D(h2)
    79  		hh2.FillColor = color.NRGBA{B: 255, A: 100}
    80  
    81  		rp.Top.Add(
    82  			hh1, hh2,
    83  			hplot.NewGrid(),
    84  		)
    85  
    86  		rp.Top.Legend.Add("template", hh1)
    87  		rp.Top.Legend.Add("monte-carlo", hh2)
    88  		rp.Top.Legend.Top = true
    89  
    90  		rp.Bottom.X.Label.Text = "X"
    91  		rp.Bottom.Y.Label.Text = "Diff"
    92  		rp.Bottom.Add(
    93  			hplot.NewH1D(hbook.SubH1D(h1, h2)),
    94  			hplot.NewGrid(),
    95  		)
    96  
    97  		err := hplot.Save(rp, 15*vg.Centimeter, -1, "testdata/rand_h1d.png")
    98  		if err != nil {
    99  			log.Fatal(err)
   100  		}
   101  	}
   102  
   103  	// Output:
   104  	// h1: mean=+0.011317 std-dev=+2.027279 +/- 0.020273
   105  	// cdf(0)= +0.5
   106  	// cdf(1)= +0.7
   107  	// h2: mean=-0.026749 std-dev=+2.022967 +/- 0.020230
   108  	// cdf(0)= +0.5
   109  	// cdf(1)= +0.7
   110  }
   111  
   112  func TestRand1DExample(t *testing.T) {
   113  	checkPlot(cmpimg.CheckPlot)(ExampleRand1D, t, "rand_h1d.png")
   114  }
   115  
   116  type chkplotFunc func(ExampleFunc func(), t *testing.T, filenames ...string)
   117  
   118  func checkPlot(f chkplotFunc) chkplotFunc {
   119  	return func(ex func(), t *testing.T, filenames ...string) {
   120  		t.Helper()
   121  		f(ex, t, filenames...)
   122  		if t.Failed() {
   123  			return
   124  		}
   125  		for _, fname := range filenames {
   126  			_ = os.Remove(path.Join("testdata", fname))
   127  		}
   128  	}
   129  }
   130  
   131  var (
   132  	_ distuv.Rander = (*hbook.Rand1D)(nil)
   133  )