go-hep.org/x/hep@v0.38.1/hbook/rand_h1d.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
     6  
     7  import (
     8  	"math/rand/v2"
     9  	"sort"
    10  
    11  	"gonum.org/v1/gonum/stat/distuv"
    12  )
    13  
    14  // Rand1D represents a 1D distribution created from a hbook.H1D histogram.
    15  type Rand1D struct {
    16  	y    *distuv.Uniform
    17  	bins []Bin1D
    18  	cdf  []float64
    19  }
    20  
    21  // NewRand1D creates a Rand1D from the provided histogram.
    22  // If src is nil, the global x/exp/rand source will be used.
    23  func NewRand1D(h *H1D, src rand.Source) Rand1D {
    24  	var (
    25  		sum  float64
    26  		n    = len(h.Binning.Bins)
    27  		cdf  = make([]float64, 1, n+1)
    28  		bins = make([]Bin1D, n, n+1)
    29  	)
    30  
    31  	copy(bins, h.Binning.Bins)
    32  	for _, xbin := range bins {
    33  		sum += xbin.SumW()
    34  		cdf = append(cdf, sum)
    35  	}
    36  	cdf[len(cdf)-1] = sum
    37  	norm := 1 / sum
    38  	for i := range cdf {
    39  		cdf[i] *= norm
    40  	}
    41  	// append a bin with similar shape to ease boundary case.
    42  	bins = append(bins, bins[n-1])
    43  	bins[n].Range.Min = bins[n-1].Range.Max
    44  	bins[n].Range.Max = bins[n-1].Range.Max + bins[n-1].Range.Width()
    45  
    46  	return Rand1D{
    47  		y:    &distuv.Uniform{Min: 0, Max: 1, Src: src},
    48  		bins: bins,
    49  		cdf:  cdf,
    50  	}
    51  }
    52  
    53  func (d *Rand1D) search(v float64) int {
    54  	return sort.Search(len(d.cdf)-1, func(i int) bool { return v < d.cdf[i+1] })
    55  }
    56  
    57  // Rand returns a random sample drawn from the distribution.
    58  func (d *Rand1D) Rand() float64 {
    59  	var (
    60  		y = d.y.Rand()
    61  		i = d.search(y)
    62  		x = d.bins[i].XEdges().Min
    63  	)
    64  	if y > d.cdf[i] {
    65  		xbin := d.bins[i+1]
    66  		cdf1 := d.cdf[i+1]
    67  		cdf0 := d.cdf[i]
    68  		x += xbin.XWidth() * (y - cdf0) / (cdf1 - cdf0)
    69  	}
    70  	return x
    71  }
    72  
    73  // CDF computes the value of the cumulative density function at x.
    74  func (d *Rand1D) CDF(x float64) float64 {
    75  	i := Bin1Ds(d.bins).IndexOf(x)
    76  	switch i {
    77  	case UnderflowBin1D:
    78  		return 0
    79  	case OverflowBin1D:
    80  		return 1
    81  	default:
    82  		return d.cdf[i]
    83  	}
    84  }