Golang math.Exp函数代码示例

func TestLegendre(t *testing.T) {
	for i, test := range []struct {
		f        func(float64) float64
		min, max float64
		n        []int
		tol      []float64
		ans      float64
	}{
		// Tolerances determined from intuition and a bit of post-hoc tweaking.
		{
			f:   func(x float64) float64 { return math.Exp(x) },
			min: -3,
			max: 5,
			n:   []int{3, 4, 6, 7, 15, 16, 300, 301},
			tol: []float64{5e-2, 5e-3, 5e-6, 1e-7, 1e-14, 1e-14, 1e-14, 1e-14},
			ans: math.Exp(5) - math.Exp(-3),
		},
	} {
		for j, n := range test.n {
			ans := Fixed(test.f, test.min, test.max, n, Legendre{}, 0)
			if !floats.EqualWithinAbsOrRel(ans, test.ans, test.tol[j], test.tol[j]) {
				t.Errorf("Mismatch. Case = %d, n = %d. Want %v, got %v", i, n, test.ans, ans)
			}
			ans2 := Fixed(test.f, test.min, test.max, n, Legendre{}, 3)
			if !floats.EqualWithinAbsOrRel(ans2, test.ans, test.tol[j], test.tol[j]) {
				t.Errorf("Mismatch concurrent. Case = %d, n = %d. Want %v, got %v", i, n, test.ans, ans)
			}
		}
	}
}

func NewActivationFunction(name ActivationName) ActivationFunction {
	switch name {
	case ActivationName_LINEAR:
		return func(x mat64.Matrix, y *mat64.Dense) { y.Clone(x) }
	case ActivationName_LOGISTIC:
		return func(x mat64.Matrix, y *mat64.Dense) {
			y.Apply(func(r, c int, v float64) float64 {
				return 1 / (1 + math.Exp(-v))
			}, x)
		}
	case ActivationName_RELU:
		return func(x mat64.Matrix, y *mat64.Dense) {
			y.Apply(func(r, c int, v float64) float64 { return math.Max(0, v) }, x)
		}
	case ActivationName_TANH:
		return func(x mat64.Matrix, y *mat64.Dense) {
			y.Apply(func(r, c int, v float64) float64 { return math.Tanh(v) }, x)
		}
	case ActivationName_SOFTMAX:
		return func(x mat64.Matrix, y *mat64.Dense) {
			r, c := x.Dims()
			for i := 0; i < r; i++ {
				exp_sum := 0.0
				for j := 0; j < c; j++ {
					exp_sum = exp_sum + math.Exp(x.At(i, j))
				}
				for j := 0; j < c; j++ {
					y.Set(i, j, math.Exp(x.At(i, j))/exp_sum)
				}
			}
		}
	}
	return nil
}

func gradient(x float64, y float64) (xd float64, yd float64) {
	xd = (2 * x * math.Cos(math.Pow(x, 2)+math.Pow(y, 2)) * math.Cos(y+math.Exp(x))) -
		(math.Exp(x) * math.Sin(math.Pow(x, 2)+math.Pow(y, 2)) * math.Sin(y+math.Exp(x)))
	yd = (2 * y * math.Cos(math.Pow(x, 2)+math.Pow(y, 2)) * math.Cos(y+math.Exp(x))) -
		(math.Sin(math.Pow(x, 2)+math.Pow(y, 2)) * math.Sin(y+math.Exp(x)))
	return
}

/* Test integrating e^(-x) over infinite domains*/
func TestNegativeExponential(t *testing.T) {
	const (
		h       = 1e-8
		correct = 1
	)

	f := func(x float64) float64 { return math.Exp(-x) }

	// Check (-Inf, 0]; should be +Inf
	if msg, ok := test_integral(f, math.Inf(-1), 0, h, math.Inf(1)); !ok {
		t.Error(msg)
	}

	// Check [0, +Inf); should be 1
	if msg, ok := test_integral(f, 0, math.Inf(1), h, 1); !ok {
		t.Error(msg)
	}

	// Now check that these results hold for -f
	f = func(x float64) float64 { return -math.Exp(-x) }

	// Check (-Inf, 0]; should be -Inf
	if msg, ok := test_integral(f, math.Inf(-1), 0, h, math.Inf(-1)); !ok {
		t.Error(msg)
	}

	// Check [0, +Inf); should be -1
	if msg, ok := test_integral(f, 0, math.Inf(1), h, -1); !ok {
		t.Error(msg)
	}
}

func Exp_mult_e(x float64, y float64, result *Result) error {
	ay := math.Abs(y)

	if y == 0.0 {
		result.val = 0.0
		result.err = 0.0
		return err.SUCCESS
	} else if (x < 0.5*gsl.LOG_DBL_MAX && x > 0.5*gsl.LOG_DBL_MIN) && (ay < 0.8*gsl.SQRT_DBL_MAX && ay > 1.2*gsl.SQRT_DBL_MIN) {
		ex := math.Exp(x)
		result.val = y * ex
		result.err = (2.0 + math.Abs(x)) * gsl.DBL_EPSILON * math.Abs(result.val)
		return err.SUCCESS
	}

	ly := math.Log(ay)
	lnr := x + ly

	if lnr > gsl.LOG_DBL_MAX-0.01 {
		return OverflowError(result)
	} else if lnr < gsl.LOG_DBL_MIN+0.01 {
		return UnderflowError(result)
	}

	sy := gsl.Sign(y)
	M := math.Floor(x)
	N := math.Floor(ly)
	a := x - M
	b := ly - N
	berr := 2.0 * gsl.DBL_EPSILON * (math.Abs(ly) + math.Abs(N))
	result.val = float64(sy) * math.Exp(M+N) * math.Exp(a+b)
	result.err = berr * math.Abs(result.val)
	result.err += 2.0 * gsl.DBL_EPSILON * (M + N + 1.0) * math.Abs(result.val)
	return err.SUCCESS
}

// NormalArea calculates the area below the normal curve.
//
// den = 2^n*n!*(2*n+1)
//
// x - x(3,5,7...) / den
//
func NormalArea_(x float64) float64 {

	a := x
	// xx = x*x
	n := 1.0
	f := 1.0

	sign := true

	for {

		f = float64(Factorial(int64(n))) * math.Pow(2, n) * (2*n + 1)
		n += 2
		println(n, f, a)

		if sign {
			a -= math.Pow(x, n) / f
			sign = false
		} else {
			a += math.Pow(x, n) / f
			sign = true
		}

		if n > 10 {
			break
		}
	}

	println("a", a)
	println("/", math.Exp(-x*x/2.0)/k)

	return a * math.Exp(-x*x/2.0) / k
}

func rbmExactExpectation(r *RBM, layer []bool, hidden bool) linalg.Vector {
	var normalizer kahan.Summer64
	var outcomeSum []kahan.Summer64
	if hidden {
		outcomeSum = make([]kahan.Summer64, rbmTestHiddenSize)
	} else {
		outcomeSum = make([]kahan.Summer64, rbmTestVisibleSize)
	}

	for i := 0; i < (1 << uint(len(outcomeSum))); i++ {
		variableVec := boolVecFromInt(i, len(outcomeSum))
		var prob float64
		if hidden {
			prob = math.Exp(-rbmEnergy(r, layer, variableVec))
		} else {
			prob = math.Exp(-rbmEnergy(r, variableVec, layer))
		}
		normalizer.Add(prob)
		for j, b := range variableVec {
			if b {
				outcomeSum[j].Add(prob)
			}
		}
	}

	expectation := make(linalg.Vector, len(outcomeSum))
	norm := 1.0 / normalizer.Sum()
	for i, s := range outcomeSum {
		expectation[i] = norm * s.Sum()
	}
	return expectation
}

func (sh *ssHeap) Less(i, j int) bool {
	ss := sh.ss
	a, b := &ss.buckets[sh.h[i]], &ss.buckets[sh.h[j]]
	rateA, rateB := a.rate, b.rate
	lastA, lastB := a.lastTs, b.lastTs

	// Formula the same as recount(), inline is faster
	if lastA >= lastB {
		// optimization. if rateB is already smaller than rateA, there
		// is no need to compute real rates. It ain't gonna grow, and
		// we can avoid running expensive math.Exp().
		if rateB >= rateA {
			rateB *= math.Exp(float64(lastA-lastB) * ss.weightHelper)
		}
	} else {
		if rateA >= rateB {
			rateA *= math.Exp(float64(lastB-lastA) * ss.weightHelper)
		}
	}

	if rateA != rateB {
		return rateA < rateB
	} else {
		// This makes difference for unitialized buckets. Rate is
		// zero, but lastTs is modified. In such case make sure to use
		// the unintialized bucket first.
		return lastA < lastB
	}
}

//Boost performs categorical adaptive boosting using the specified partition and
//returns the weight that tree that generated the partition should be given.
func (t *AdaBoostTarget) Boost(leaves *[][]int) (weight float64) {
	weight = 0.0
	for _, cases := range *leaves {
		weight += t.Impurity(&cases, nil)
	}
	if weight >= .5 {
		return 0.0
	}
	weight = .5 * math.Log((1-weight)/weight)

	for _, cases := range *leaves {
		m := t.Modei(&cases)
		for _, c := range cases {
			if t.IsMissing(c) == false {
				cat := t.Geti(c)
				if cat != m {
					t.Weights[c] = t.Weights[c] * math.Exp(weight)
				} else {
					t.Weights[c] = t.Weights[c] * math.Exp(-weight)
				}
			}

		}
	}
	normfactor := 0.0
	for _, v := range t.Weights {
		normfactor += v
	}
	for i, v := range t.Weights {
		t.Weights[i] = v / normfactor
	}
	return
}

// Init initialises the model
func (o *RefIncRL1) Init(prms Prms) (err error) {

	// parameters
	for _, p := range prms {
		switch p.N {
		case "lam0":
			o.λ0 = p.V
		case "lam1":
			o.λ1 = p.V
		case "alp":
			o.α = p.V
		case "bet":
			o.β = p.V
		default:
			return chk.Err("ref-inc-rl1: parameter named %q is invalid", p.N)
		}
	}

	// set b
	o.b = -1.0 // not flipped
	if o.λ1 < o.λ0 {
		o.b = 1.0 // flipped
	}

	// constants
	o.c1 = o.β * o.b * (o.λ1 - o.λ0)
	o.c2 = math.Exp(o.β * o.b * o.α)
	o.c3 = math.Exp(o.β*o.b*(1.0-o.λ0)) - o.c2*math.Exp(o.c1)
	return
}

// Upper incomplete gamma.
func ugamma(x, s float64, regularized bool) float64 {
	if x <= 1.1 || x <= s {
		if regularized {
			return 1 - lgamma(x, s, regularized)
		}
		return math.Gamma(s) - lgamma(x, s, regularized)
	}

	f := 1.0 + x - s
	C := f
	D := 0.0
	var a, b, chg float64

	for i := 1; i < 10000; i++ {
		a = float64(i) * (s - float64(i))
		b = float64(i<<1) + 1.0 + x - s
		D = b + a*D
		C = b + a/C
		D = 1.0 / D
		chg = C * D
		f *= chg
		if math.Abs(chg-1) < eps {
			break
		}
	}
	if regularized {
		logg, _ := math.Lgamma(s)
		return math.Exp(s*math.Log(x) - x - logg - math.Log(f))
	}
	return math.Exp(s*math.Log(x) - x - math.Log(f))
}

// Evolve object within likelihood constraint
// logLstar: Likelihood constraint L > Lstar
func (Obj *Object) Explore(logLstar float64) {
	step := 0.1    // Initial guess suitable step-size in (0,1)
	m := 20        // MCMC counter (pre-judged # steps)
	accept := 0    // # MCMC acceptances
	reject := 0    // # MCMC rejections
	var Try Object // Trial object
	for ; m > 0; m-- {
		// Trial object
		Try.u = Obj.u + step*(2.*rand.Float64()-1.) // |move| < step
		Try.v = Obj.v + step*(2.*rand.Float64()-1.) // |move| < step
		Try.u -= math.Floor(Try.u)                  // wraparound to stay within (0,1)
		Try.v -= math.Floor(Try.v)                  // wraparound to stay within (0,1)
		Try.x = 4.0*Try.u - 2.0                     // map to x
		Try.y = 2.0 * Try.v                         // map to y
		Try.logL = logLhood(Try.x, Try.y)           // trial likelihood value
		// Accept if and only if within hard likelihood constraint
		if Try.logL > logLstar {
			*Obj = Try
			accept++
		} else {
			reject++
		}
		// Refine step-size to let acceptance ratio converge around 50%
		if accept > reject {
			step *= math.Exp(1.0 / float64(accept))
		}
		if accept < reject {
			step /= math.Exp(1.0 / float64(reject))
		}
	}
}

// Init initialises the model
func (o *RefDecSp1) Init(prms Prms) (err error) {

	// parameters
	e := prms.Connect(&o.β, "bet", "ref-dec-sp1 function")
	e += prms.Connect(&o.λ1, "lam1", "ref-dec-sp1 function")
	e += prms.Connect(&o.ya, "ya", "ref-dec-sp1 function")
	e += prms.Connect(&o.yb, "yb", "ref-dec-sp1 function")
	if e != "" {
		err = chk.Err("%v\n", e)
		return
	}

	// check
	if o.yb >= o.ya {
		return chk.Err("yb(%g) must be smaller than ya(%g)", o.yb, o.ya)
	}

	// constants
	o.c1 = o.β * o.λ1
	o.c2 = math.Exp(-o.β * o.ya)
	o.c3 = math.Exp(-o.β*o.yb) - o.c2
	o.c1timestmax = 400

	// check
	if math.IsInf(o.c2, 0) || math.IsInf(o.c3, 0) {
		return chk.Err("β*ya or β*yb is too large:\n β=%v, ya=%v, yb=%v\n c1=%v, c2=%v, c3=%v", o.β, o.ya, o.yb, o.c1, o.c2, o.c3)
	}
	return
}

func Anneal(state Annealable, maxTemp, minTemp float64, steps int) Annealable {
	factor := -math.Log(maxTemp / minTemp)
	state = state.Copy()
	bestState := state.Copy()
	bestEnergy := state.Energy()
	previousEnergy := bestEnergy
	for step := 0; step < steps; step++ {
		if step%100000 == 0 {
			showProgress(step, steps)
		}
		pct := float64(step) / float64(steps-1)
		temp := maxTemp * math.Exp(factor*pct)
		undo := state.DoMove()
		energy := state.Energy()
		change := energy - previousEnergy
		if change > 0 && math.Exp(-change/temp) < rand.Float64() {
			state.UndoMove(undo)
		} else {
			previousEnergy = energy
			if energy < bestEnergy {
				// pct := float64(step*100) / float64(steps)
				// fmt.Printf("step: %d of %d (%.1f%%), temp: %.3f, energy: %.3f\n",
				//     step, steps, pct, temp, energy)
				bestEnergy = energy
				bestState = state.Copy()
			}
		}
	}
	showProgress(steps, steps)
	return bestState
}

// Init initialises the function
func (o *RefDecGen) Init(prms Prms) (err error) {

	// parameters
	for _, p := range prms {
		switch p.N {
		case "bet":
			o.β = p.V
		case "a":
			o.a = p.V
		case "b":
			o.b = p.V
		case "c":
			o.c = p.V
		case "A":
			o.A = p.V
		case "B":
			o.B = p.V
		case "xini":
			o.xini = p.V
		case "yini":
			o.yini = p.V
		default:
			return chk.Err("ref-dec-gen: parameter named %q is invalid", p.N)
		}
	}

	// constants
	o.c1 = o.β * (o.b*o.A - o.a)
	o.c2 = ((o.A - o.B) / (o.A - o.a/o.b)) * math.Exp(-o.β*o.c)
	o.c3 = math.Exp(o.β*o.b*(o.yini+o.A*o.xini)) - o.c2*math.Exp(o.c1*o.xini)
	return
}