Commits

Sebastien Binet  committed 52f4f9c Merge

applied go1rename

  • Participants
  • Parent commits 0eb91ca, 0e8b94e
  • Tags go.weekly.2011-11-09

Comments (0)

Files changed (8)

File pkg/hep/yoda/Makefile

 
 TARG=bitbucket.org/binet/go-hep/pkg/hep/yoda
 GOFILES=\
+	axis.go\
+	bin.go\
+	histo1d.go\
+	mathutils.go\
+	object.go\
+	point2d.go\
 	yoda.go\
 
 include ${GOROOT}/src/Make.pkg

File pkg/hep/yoda/axis.go

+package yoda
+
+import "sort"
+
+// Axis1D is a container of ordered bins
+type Axis1D struct {
+	// the bins contained in this histogram
+	bins []hbin1d
+	// a distribution counter for underflow fills
+	underflow dbn1d
+	// a distribution counter for overflow fills
+	overflow dbn1d
+
+	// a distribution counter for the whole histogram
+	dbn dbn1d
+
+	// bin edges: lower edges, except last entry, which is the high edge of the last entry
+	edges []float64
+
+	// hash table for fast bin lookup
+	hash map[float64]int
+}
+
+// Create a new Axis1D from a list of bin edges
+func NewAxis1DFromEdges(edges []float64) *Axis1D {
+	nbins := len(edges) - 1
+	a := &Axis1D{bins: make([]hbin1d, 0, nbins),
+		underflow: dbn1d{},
+		overflow:  dbn1d{},
+		dbn:       dbn1d{},
+		edges:     make([]float64, 0, len(edges)),
+		hash:      make(map[float64]int),
+	}
+	copy(a.edges, edges)
+	for i := 0; i < nbins; i++ {
+		a.bins = append(a.bins, hbin1d{*NewBin1D(edges[i], edges[i+1])})
+	}
+	sort.Sort((sorted_hbin1ds)(a.bins))
+	for i := range a.bins {
+		// insert upper bound mapped to bin id
+		a.hash[edges[i+1]] = i
+	}
+	return a
+}
+
+// Create a new Axis1D from a number of bins and a bin distribution
+func NewAxis1D(nbins int, lower, upper float64) *Axis1D {
+	return NewAxis1DFromEdges(linspace(lower, upper, nbins))
+}
+
+// Returns the number of bins (not counting under|over-flows)
+func (a *Axis1D) NumBins() uint64 {
+	return uint64(len(a.bins))
+}
+
+// Returns the bins' axis
+func (a *Axis1D) Bins() []hbin1d {
+	return a.bins
+}
+
+// Returns the edges of the bin number 'id'
+func (a *Axis1D) BinEdges(id int) (float64, float64) {
+	return a.edges[id], a.edges[id+1]
+}
+
+// Returns the low edge of the axis
+func (a *Axis1D) LowEdge() float64 {
+	return a.bins[0].LowEdge()
+}
+
+// Returns the high edge of the axis
+func (a *Axis1D) HighEdge() float64 {
+	return a.bins[len(a.bins)-1].LowEdge()
+}
+
+// Returns the bin at bin number 'id'
+func (a *Axis1D) Bin(id int) *hbin1d {
+	return &a.bins[id]
+}
+
+// Returns the bin at coordinate 'x'
+func (a *Axis1D) BinByCoord(x float64) *hbin1d {
+	id, ok := a.hash[x]
+	if ok {
+		return &a.bins[id]
+	}
+	//FIXME: O(N) search
+	for i := range a.edges[1:] {
+		if a.edges[i-1] <= x && a.edges[i] < x {
+			id = i - 1
+			a.hash[x] = id
+			return &a.bins[id]
+		}
+	}
+	return nil
+}
+
+// Reset the axis content
+func (a *Axis1D) Reset() {
+	a.dbn.Reset()
+	a.underflow.Reset()
+	a.overflow.Reset()
+	for i := range a.bins {
+		a.bins[i].Reset()
+	}
+}
+
+// Scale the axis coordinates
+func (a *Axis1D) ScaleW(scale float64) {
+	a.dbn.scaleW(scale)
+	a.underflow.scaleW(scale)
+	a.overflow.scaleW(scale)
+	for i := range a.bins {
+		a.bins[i].scaleW(scale)
+	}
+}
+
+// In-place add of 2 axes
+func Axis1D_IAdd(a, b *Axis1D) error {
+	//FIXME error handling
+	if len(a.bins) != len(b.bins) {
+		panic("axes lengthes differ")
+	}
+	for i := range a.bins {
+		err := hbin1d_iadd(&a.bins[i], &b.bins[i])
+		if err != nil {
+			return err
+		}
+	}
+	return dbn1d_iadd(&a.dbn, &b.dbn)
+}

File pkg/hep/yoda/bin.go

+package yoda
+
+import (
+	"errors"
+	"math"
+
+	//"sort"
+)
+
+// Base class for bins in 1D normal and profile histograms.
+// The lower bin edge is inclusive
+type Bin1D struct {
+	// the bin limits
+	edges [2]float64
+
+	// distribution of weighted x-values
+	xdbn dbn1d
+}
+
+// Create a new Bin1D given low and high edges
+func NewBin1D(low, hi float64) *Bin1D {
+	return &Bin1D{edges: [2]float64{low, hi}, xdbn: dbn1d{}}
+}
+
+// XMin returns the lower limit of the bin (inclusive)
+func (b *Bin1D) XMin() float64 {
+	return b.edges[0]
+}
+
+// Returns the lower limit of the bin (inclusive)
+func (b *Bin1D) LowEdge() float64 {
+	return b.edges[0]
+}
+
+// XMax returns the upper limit of the bin (exclusive)
+func (b *Bin1D) XMax() float64 {
+	return b.edges[1]
+}
+
+// Returns the upper limit of the bin (exclusive)
+func (b *Bin1D) HighEdge() float64 {
+	return b.edges[1]
+}
+
+// Edges returns the lower and upper limits of the bin
+func (b *Bin1D) Edges() (float64, float64) {
+	return b.edges[0], b.edges[1]
+}
+
+// Reset this object
+func (b *Bin1D) Reset() {
+	b.xdbn.Reset()
+}
+
+// Returns the separation of low and high edges, i,e, high - low
+func (b *Bin1D) Width() float64 {
+	return b.edges[1] - b.edges[0]
+}
+
+// Returns the mean position in the bin, or the midpoint if that is not available
+func (b *Bin1D) Focus() float64 {
+	if b.xdbn.sumw != 0.0 {
+		return b.XMean()
+	}
+	return b.MidPoint()
+}
+
+// Returns the geometric centre of the bin, i.e. high+low/2.0
+func (b *Bin1D) MidPoint() float64 {
+	return (b.edges[1] + b.edges[0]) / 2.0
+}
+
+// Returns the mean value of x-values in the bin.
+func (b *Bin1D) XMean() float64 {
+	return b.xdbn.mean()
+}
+
+// Returns the variance of x-values in the bin
+func (b *Bin1D) XVariance() float64 {
+	v, err := b.xdbn.variance()
+	if err != nil {
+		panic(err)
+	}
+	return v
+}
+
+// Returns the standard deviation (spread) of x-values in the bin
+func (b *Bin1D) XStdDev() float64 {
+	v, err := b.xdbn.stdDev()
+	if err != nil {
+		panic(err)
+	}
+	return v
+}
+
+// Returns the standard error on the bin focus
+func (b *Bin1D) XStdError() float64 {
+	v, err := b.xdbn.stdErr()
+	if err != nil {
+		panic(err)
+	}
+	return v
+}
+
+// Returns the number of entries in the bin
+func (b *Bin1D) NumEntries() uint64 {
+	return b.xdbn.nfills
+}
+
+// Returns the sum of weights
+func (b *Bin1D) SumW() float64 {
+	return b.xdbn.sumw
+}
+
+// Returns the sum of weights squared
+func (b *Bin1D) SumW2() float64 {
+	return b.xdbn.sumw2
+}
+
+// Returns the sum of x*weight
+func (b *Bin1D) SumWX() float64 {
+	return b.xdbn.sumwx
+}
+
+// Returns the sum of x**2 * weight
+func (b *Bin1D) SumWX2() float64 {
+	return b.xdbn.sumwx2
+}
+
+func (b *Bin1D) scaleW(scale float64) {
+	b.xdbn.scaleW(scale)
+}
+
+// Returns a new bin, the sum of the input bins
+// if the edges do not match, nil is returned
+func Bin1D_Add(a, b *Bin1D) *Bin1D {
+	if a.edges[0] != b.edges[0] ||
+		a.edges[1] != b.edges[1] {
+		return nil
+	}
+	return &Bin1D{edges: a.edges, xdbn: *dbn1d_add(&a.xdbn, &b.xdbn)}
+}
+
+// in-place sum of the input bins: a += b
+func Bin1D_IAdd(a, b *Bin1D) error {
+	if a.edges[0] != b.edges[0] ||
+		a.edges[1] != b.edges[1] {
+		return errors.New("iadd: axes' edges to not match")
+	}
+	return dbn1d_iadd(&a.xdbn, &b.xdbn)
+}
+
+// Returns a new bin, the subtraction of the input bins
+func Bin1D_Sub(a, b *Bin1D) *Bin1D {
+	if a.edges[0] != b.edges[0] ||
+		a.edges[1] != b.edges[1] {
+		return nil
+	}
+	return &Bin1D{edges: a.edges, xdbn: *dbn1d_sub(&a.xdbn, &b.xdbn)}
+}
+
+// Compares 2 Bin1Ds, by lower edge position
+// FIXME: check for overlap somewhere...
+func Bin1D_Less(a, b *Bin1D) bool {
+	return a.edges[0] < b.edges[0]
+}
+
+// dbn1d is a 1D distribution.
+// This class is used to centralize the calculation of statistics of unbounded,
+// unbinned sampled distributions.
+// Each distribution fill contributes a weight 'w' and a value 'x'.
+// By storing the total number of fills (ignoring weights), Sum(w), Sum(w^2),
+// Sum(w.x) and Sum(w.x^2), the dbn1d can calculate the mean and spread 
+// \sigma^2, \sigma and \hat{\sigma} of the sampled distribution.
+// It is used to provide this information in bins for the "hidden" 'y' 
+// distribution in profile histogram bins.
+type dbn1d struct {
+	nfills uint64
+	sumw   float64
+	sumw2  float64
+	sumwx  float64
+	sumwx2 float64
+}
+
+func (d *dbn1d) scaleW(factor float64) {
+	sf2 := factor * factor
+	d.sumw *= factor
+	d.sumw2 *= sf2
+	d.sumwx *= factor
+	d.sumwx2 *= sf2
+}
+
+func (d *dbn1d) fill(val, weight float64) {
+	d.nfills += 1
+	d.sumw += weight
+	w2 := weight * weight
+	if weight < 0. {
+		w2 *= -1.0
+	}
+	d.sumw2 += w2
+	d.sumwx += weight * val
+	d.sumwx2 += weight * val * val
+}
+
+func (d *dbn1d) Reset() {
+	d.nfills = 0
+	d.sumw = 0.0
+	d.sumw2 = 0.0
+	d.sumwx = 0.0
+	d.sumwx2 = 0.0
+}
+
+func (d *dbn1d) effNumEntries() float64 {
+	return d.sumw * d.sumw / d.sumw2
+}
+
+func (d *dbn1d) mean() float64 {
+	return d.sumwx / d.sumw
+}
+
+// The weighted variance is defined as:
+//  sig2 = (sum(wx**2) * sum(w) - sum(wx)**2) / (sum(w)**2 - sum(w**2))
+//  http://en.wikipedia.org/wiki/Weighted_mean
+func (d *dbn1d) variance() (float64, error) {
+	effn := d.effNumEntries()
+	if effn == 0.0 {
+		return 0.0, errors.New("requested width of a distribution with no net fill weights")
+	}
+	if effn <= 1.0 {
+		return 0.0, errors.New("requested width of a distribution with only one effective entry")
+	}
+	num := d.sumwx2*d.sumw - d.sumwx*d.sumwx
+	den := d.sumw*d.sumw - d.sumw2
+	if den == 0 {
+		panic("undefined weighted variance")
+	}
+	if math.Abs(num) < 1e-10 && math.Abs(den) < 1e-10 {
+		panic("numerically unstable weights in width calculation")
+	}
+	return num / den, nil
+}
+
+func (d *dbn1d) stdDev() (float64, error) {
+	v, err := d.variance()
+	if err != nil {
+		return 0.0, err
+	}
+	return math.Sqrt(v), nil
+}
+
+func (d *dbn1d) stdErr() (float64, error) {
+	effnum := d.effNumEntries()
+	if effnum == 0.0 {
+		return 0.0, errors.New("requested std error of a distribution with no net fill weights")
+	}
+	v, err := d.variance()
+	if err != nil {
+		return 0.0, err
+	}
+	return math.Sqrt(v / effnum), nil
+}
+
+func dbn1d_add(a, b *dbn1d) *dbn1d {
+	return &dbn1d{
+		nfills: a.nfills + b.nfills,
+		sumw:   a.sumw + b.sumw,
+		sumw2:  a.sumw2 + b.sumw2,
+		sumwx:  a.sumwx + b.sumwx,
+		sumwx2: a.sumwx2 + b.sumwx2,
+	}
+}
+
+func dbn1d_iadd(a, b *dbn1d) error {
+	a.nfills += b.nfills
+	a.sumw += b.sumw
+	a.sumw2 += b.sumw2
+	a.sumwx += b.sumwx
+	a.sumwx2 += b.sumwx2
+	return nil
+}
+
+func dbn1d_sub(a, b *dbn1d) *dbn1d {
+	return &dbn1d{
+		nfills: a.nfills + b.nfills,
+		sumw:   a.sumw - b.sumw,
+		sumw2:  a.sumw2 - b.sumw2,
+		sumwx:  a.sumwx - b.sumwx,
+		sumwx2: a.sumwx2 - b.sumwx2,
+	}
+}
+
+func dbn1d_isub(a, b *dbn1d) error {
+	a.nfills += b.nfills
+	a.sumw -= b.sumw
+	a.sumw2 -= b.sumw2
+	a.sumwx -= b.sumwx
+	a.sumwx2 -= b.sumwx2
+	return nil
+}
+
+// a list of sorted Bin1Ds
+type SortedBin1Ds []Bin1D
+
+func (s SortedBin1Ds) Len() int {
+	return len(s)
+}
+
+func (s SortedBin1Ds) Less(i, j int) bool {
+	return Bin1D_Less(&s[i], &s[j])
+}
+
+func (s SortedBin1Ds) Swap(i, j int) {
+	s[i], s[j] = s[j], s[i]
+}

File pkg/hep/yoda/histo1d.go

+package yoda
+
+import (
+	"errors"
+	"math"
+)
+
+type hbin1d struct {
+	Bin1D
+}
+
+// Fill this bin with weight 'weight' at position 'coord'
+func (h *hbin1d) fill(coord, weight float64) {
+	h.Bin1D.xdbn.fill(coord, weight)
+}
+
+// Fill this bin with weight 'weight'
+func (h *hbin1d) fillBin(weight float64) {
+	h.Bin1D.xdbn.fill(h.Bin1D.MidPoint(), weight)
+}
+
+// The area of a bin is the sum of weights in the bin, i.e. the
+// width of the bin has no influence on this figure
+func (h *hbin1d) area() float64 {
+	return h.Bin1D.SumW()
+}
+
+// The height is defined as area/width
+func (h *hbin1d) height() float64 {
+	return h.area() / h.Bin1D.Width()
+}
+
+// Error computed using binomial statistics on the sum of bin weights
+// i.e. err_area = sqrt{sum{weights}}
+func (h *hbin1d) areaError() float64 {
+	return math.Sqrt(h.Bin1D.SumW())
+}
+
+// Error includes scaling factor of the bin width
+// i.e. err_height = sqrt{sum{weight}} / width
+func (h *hbin1d) heightError() float64 {
+	return h.areaError() / h.Bin1D.Width()
+}
+
+// Returns a new bin, the sum of the input bins
+// if the edges do not match, nil is returned
+func hbin1d_add(a, b *hbin1d) *hbin1d {
+	if a.Bin1D.edges[0] != b.Bin1D.edges[0] ||
+		a.Bin1D.edges[1] != b.Bin1D.edges[1] {
+		return nil
+	}
+	return &hbin1d{Bin1D{edges: a.Bin1D.edges, xdbn: *dbn1d_add(&a.Bin1D.xdbn, &b.Bin1D.xdbn)}}
+}
+
+// in-place sum of the input bins: a += b
+func hbin1d_iadd(a, b *hbin1d) error {
+	if a.Bin1D.edges[0] != b.Bin1D.edges[0] ||
+		a.Bin1D.edges[1] != b.Bin1D.edges[1] {
+		return errors.New("iadd: axes' edges to not match")
+	}
+	return dbn1d_iadd(&a.Bin1D.xdbn, &b.Bin1D.xdbn)
+}
+
+// Returns a new bin, the subtraction of the input bins
+func hbin1d_sub(a, b *hbin1d) *hbin1d {
+	if a.Bin1D.edges[0] != b.Bin1D.edges[0] ||
+		a.Bin1D.edges[1] != b.Bin1D.edges[1] {
+		return nil
+	}
+	return &hbin1d{Bin1D{edges: a.Bin1D.edges, xdbn: *dbn1d_sub(&a.Bin1D.xdbn, &b.Bin1D.xdbn)}}
+}
+
+// Compares 2 hbin1ds, by lower edge position
+// FIXME: check for overlap somewhere...
+func hbin1d_less(a, b *hbin1d) bool {
+	return a.Bin1D.edges[0] < b.Bin1D.edges[0]
+}
+
+// a list of sorted hbin1d
+type sorted_hbin1ds []hbin1d
+
+func (s sorted_hbin1ds) Len() int {
+	return len(s)
+}
+
+func (s sorted_hbin1ds) Less(i, j int) bool {
+	return hbin1d_less(&s[i], &s[j])
+}
+
+func (s sorted_hbin1ds) Swap(i, j int) {
+	s[i], s[j] = s[j], s[i]
+}
+// A one-dimensional histogram
+type Histo1D struct {
+	axis Axis1D
+}

File pkg/hep/yoda/mathutils.go

+package yoda
+
+import (
+	"math"
+)
+
+// Compare 2 floating point numbers for equality with a degree of fuzziness
+// The tolreance parameter is fractional.
+func fFuzzyEq(a, b float64) bool {
+	return fFuzzyEqWithTolerance(a, b, 1e-5)
+}
+
+// Compare 2 floating point numbers for equality with a degree of fuzziness
+// The tolreance parameter is fractional.
+func fFuzzyEqWithTolerance(a, b, tolerance float64) bool {
+	absavg := math.Abs(a+b) / 2.0
+	absdiff := math.Abs(a - b)
+	return (absavg == 0.0 && absdiff == 0.0) || (absdiff < tolerance*absavg)
+}
+
+// Returns a list of nbins+1 values equally spaced between 'start' and 'end' inclusive
+func linspace(start, end float64, nbins int) []float64 {
+	if end < start {
+		//FIXME(binet) proper error handling
+		panic("invalid range")
+	}
+	if nbins <= 0 {
+		//FIXME(binet) proper error handling
+		panic("invalid range")
+	}
+	o := make([]float64, 0, nbins+1)
+	interval := (end - start) / float64(nbins)
+	for edge := start; edge <= end; edge += interval {
+		o = append(o, edge)
+	}
+	if len(o) != nbins+1 {
+		//FIXME(binet) proper error handling
+		panic("internal error")
+	}
+	return o
+}
+
+/// Calculates the mean of a sample
+func mean(sample []int) float64 {
+	m := 0.0
+	for _, i := range sample {
+		m += float64(i)
+	}
+	return m / float64(len(sample))
+}
+
+/// Calculates the covariance (variance) between two samples
+func covariance(s1, s2 []int) float64 {
+	m1 := mean(s1)
+	m2 := mean(s2)
+	n := len(s1)
+	if n > len(s2) {
+		n = len(s2)
+	}
+	cov := 0.0
+	for i := 0; i < n; i++ {
+		c := (float64(s1[i]) - m1) * (float64(s2[i]) - m2)
+		cov += c
+	}
+	if n > 1 {
+		return cov / float64(n-1)
+	}
+	return 0.0
+}
+
+/// Calculates the correlation strength between two samples
+func correlation(s1, s2 []int) float64 {
+	//FIXME: make this cache-friendly by iterating only ONCE!
+	cov := covariance(s1, s2)
+	var1 := covariance(s1, s1)
+	var2 := covariance(s2, s2)
+	corr := cov / math.Sqrt(var1*var2)
+	return corr * math.Sqrt(var2/var1)
+}

File pkg/hep/yoda/object.go

+package yoda
+
+type Annotations map[string]interface{}
+
+type Object interface {
+	// Reset this analysis object
+	Reset()
+
+	// Retrieve the annotations attached to this analysis object
+	Annotations() Annotations
+
+	// Check if an annotation 'name' is defined
+	HasAnnotation(name string) bool
+}
+
+type Bin interface {
+	// Reset this bin
+	Reset()
+
+	// NumEntries returns the number of entries for this bin
+	NumEntries() uint64
+
+	// SumW returns the sum of weights
+	SumW() float64
+
+	// SumW2 returns the sum of weights squared
+	SumW2() float64
+}
+
+type obj_impl struct {
+	ann Annotations
+}
+

File pkg/hep/yoda/point2d.go

+package yoda
+
+// A 2D data point to be contained in a Scatter2D
+type Point2D struct {
+	coord [2]float64
+	err   [2][2]float64
+}
+
+func NewPoint2D(x, y float64) *Point2D {
+	return NewPoint2DErr(x, y, 0., 0.)
+}
+
+func NewPoint2DErr(x, y, ex, ey float64) *Point2D {
+	return &Point2D{
+		coord: [2]float64{x, y},
+		err:   [2][2]float64{{ex, ex}, {ey, ey}},
+	}
+}
+
+func NewPoint2DAsymErr(x, y, exmin, exmax, eymin, eymax float64) *Point2D {
+	return &Point2D{
+		coord: [2]float64{x, y},
+		err:   [2][2]float64{{exmin, exmax}, {eymin, eymax}},
+	}
+}
+
+// Get x-coordinate
+func (p *Point2D) X() float64 {
+	return p.coord[0]
+}
+
+// Get x-value minus negative x-error
+func (p *Point2D) XMin() float64 {
+	return p.coord[0] - p.err[0][0]
+}
+
+// Get x-value plus positive x-error
+func (p *Point2D) XMax() float64 {
+	return p.coord[0] + p.err[0][1]
+}
+
+// Get x-error values
+func (p *Point2D) XErrs() (float64, float64) {
+	return p.err[0][0], p.err[0][1]
+}
+
+// Get negative x-error value
+func (p *Point2D) XErrMinus() float64 {
+	return p.err[0][0]
+}
+
+// Get positive x-error value
+func (p *Point2D) XErrPlus() float64 {
+	return p.err[0][1]
+}
+
+// Get average x-error value
+func (p *Point2D) XErrAvg() float64 {
+	return (p.err[0][0] + p.err[0][1]) / float64(2.)
+}
+
+// Set x-coordinate
+func (p *Point2D) SetX(x float64) {
+	p.coord[0] = x
+}
+
+// Set symmetric x-error
+func (p *Point2D) SetXErr(ex float64) {
+	p.err[0][0] = ex
+	p.err[0][1] = ex
+}
+
+// Set asymmetric x-error
+func (p *Point2D) SetXErrs(exmin, exmax float64) {
+	p.err[0][0] = exmin
+	p.err[0][1] = exmax
+}
+
+// Get y-coordinate
+func (p *Point2D) Y() float64 {
+	return p.coord[1]
+}
+
+// Get y-value minus negative y-error
+func (p *Point2D) YMin() float64 {
+	return p.coord[1] - p.err[1][0]
+}
+
+// Get y-value plus positive y-error
+func (p *Point2D) YMax() float64 {
+	return p.coord[1] + p.err[1][1]
+}
+
+// Get y-error values
+func (p *Point2D) YErrs() (float64, float64) {
+	return p.err[1][0], p.err[1][1]
+}
+
+// Get negative y-error value
+func (p *Point2D) YErrMinus() float64 {
+	return p.err[1][0]
+}
+
+// Get positive y-error value
+func (p *Point2D) YErrPlus() float64 {
+	return p.err[1][1]
+}
+
+// Get average y-error value
+func (p *Point2D) YErrAvg() float64 {
+	return (p.err[1][0] + p.err[1][1]) / float64(2.)
+}
+
+// Set y-coordinate
+func (p *Point2D) SetY(y float64) {
+	p.coord[1] = y
+}
+
+// Set symmetric y-error
+func (p *Point2D) SetYErr(ey float64) {
+	p.err[1][0] = ey
+	p.err[1][1] = ey
+}
+
+// Set asymmetric y-error
+func (p *Point2D) SetYErrs(eymin, eymax float64) {
+	p.err[1][0] = eymin
+	p.err[1][1] = eymax
+}
+
+/// Equality test (of x-characteristics only)
+func Point2D_Eq(p1, p2 *Point2D) bool {
+	return fFuzzyEq(p1.X(), p2.X()) &&
+		fFuzzyEq(p1.XErrMinus(), p2.XErrMinus()) &&
+		fFuzzyEq(p1.XErrPlus(), p2.XErrPlus())
+}

File pkg/hep/yoda/yoda.go

 package yoda
 
-type Point2D struct {
-	coord [2]float64
-	errs  [2][2]float64
-}
-
-func NewPoint2D(x, y float64) *Point2D {
-	return &Point2D{coord: [2]float64{x, y}, errs: [2][2]float64{}}
-}
-
-func NewPoint2DErr(x, y, ex, ey float64) *Point2D {
-	return &Point2D{
-		coord: [2]float64{x, y},
-		errs:  [2][2]float64{{-ex, +ex}, {-ey, +ey}},
-	}
-}