func Test_2dinteg02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("2dinteg02. bidimensional integral")

	// Γ(1/4, 1)
	gamma_1div4_1 := 0.2462555291934987088744974330686081384629028737277219

	x := utl.LinSpace(0, 1, 11)
	y := utl.LinSpace(0, 1, 11)
	m, n := len(x), len(y)
	f := la.MatAlloc(m, n)
	for i := 0; i < m; i++ {
		for j := 0; j < n; j++ {
			f[i][j] = 8.0 * math.Exp(-math.Pow(x[i], 2)-math.Pow(y[j], 4))
		}
	}
	dx, dy := x[1]-x[0], y[1]-y[0]
	Vt := Trapz2D(dx, dy, f)
	Vs := Simps2D(dx, dy, f)
	Vc := math.Sqrt(math.Pi) * math.Erf(1) * (math.Gamma(1.0/4.0) - gamma_1div4_1)
	io.Pforan("Vt = %v\n", Vt)
	io.Pforan("Vs = %v\n", Vs)
	io.Pfgreen("Vc = %v\n", Vc)
	chk.Scalar(tst, "Vt", 0.0114830435645548, Vt, Vc)
	chk.Scalar(tst, "Vs", 1e-4, Vs, Vc)

}

// Plot plots retention model
//  args1 -- arguments for model computed by solving differential equation; e.g. "'b*-'"
//           if args1 == "", plot is skiped
//  args2 -- arguments for model computed by directly calling sl(pc), if available
//           if args2 == "", plot is skiped
func Plot(mdl Model, pc0, sl0, pcf float64, npts int, args1, args2, label string) (err error) {

	// plot using Update
	Pc := utl.LinSpace(pc0, pcf, npts)
	Sl := make([]float64, npts)
	if args1 != "" {
		Sl[0] = sl0
		for i := 1; i < npts; i++ {
			Sl[i], err = Update(mdl, Pc[i-1], Sl[i-1], Pc[i]-Pc[i-1])
			if err != nil {
				return
			}
		}
		plt.Plot(Pc, Sl, io.Sf("%s, label='%s', clip_on=0", args1, label))
	}

	// plot using Sl function
	if args2 != "" {
		if m, ok := mdl.(Nonrate); ok {
			Pc = utl.LinSpace(pc0, pcf, 101)
			Sl = make([]float64, 101)
			for i, pc := range Pc {
				Sl[i] = m.Sl(pc)
			}
			plt.Plot(Pc, Sl, io.Sf("%s, label='%s_direct', clip_on=0", args2, label))
		}
	}
	return
}

func Test_cxdeb01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("cxdeb01. Deb's crossover")

	var ops OpsData
	ops.SetDefault()
	ops.Pc = 1.0
	ops.Xrange = [][]float64{{-3, 3}, {-4, 4}}
	ops.EnfRange = true

	rnd.Init(0)

	A := []float64{-1, 1}
	B := []float64{1, 2}
	a := make([]float64, len(A))
	b := make([]float64, len(A))
	FltCrossoverDeb(a, b, A, B, 0, &ops)
	io.Pforan("A = %v\n", A)
	io.Pforan("B = %v\n", B)
	io.Pfcyan("a = %.6f\n", a)
	io.Pfcyan("b = %.6f\n", b)

	nsamples := 1000
	a0s, a1s := make([]float64, nsamples), make([]float64, nsamples)
	b0s, b1s := make([]float64, nsamples), make([]float64, nsamples)
	for i := 0; i < nsamples; i++ {
		FltCrossoverDeb(a, b, B, A, 0, &ops)
		a0s[i], a1s[i] = a[0], a[1]
		b0s[i], b1s[i] = b[0], b[1]
	}
	ha0 := rnd.Histogram{Stations: []float64{-4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1}}
	hb0 := rnd.Histogram{Stations: []float64{0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 5, 5.5, 6}}
	ha1 := rnd.Histogram{Stations: utl.LinSpace(-4, 4, 11)}
	hb1 := rnd.Histogram{Stations: utl.LinSpace(-4, 4, 11)}
	ha0.Count(a0s, true)
	hb0.Count(b0s, true)
	ha1.Count(a1s, true)
	hb1.Count(b1s, true)

	io.Pforan("\na0s\n")
	io.Pf("%s", rnd.TextHist(ha0.GenLabels("%.1f"), ha0.Counts, 60))
	io.Pforan("b0s\n")
	io.Pf("%s", rnd.TextHist(hb0.GenLabels("%.1f"), hb0.Counts, 60))

	io.Pforan("\na1s\n")
	io.Pf("%s", rnd.TextHist(ha1.GenLabels("%.1f"), ha1.Counts, 60))
	io.Pforan("b1s\n")
	io.Pf("%s", rnd.TextHist(hb1.GenLabels("%.1f"), hb1.Counts, 60))
}

// PlotStress plots stresses along y=0 (horizontal line) and x=0 (vertical line)
func (o *PlateHole) PlotStress(t, L float64, npts int) {

	d := utl.LinSpace(o.r, L, npts)
	Sx := make([]float64, npts)
	Sy := make([]float64, npts)
	Sxy := make([]float64, npts)

	plt.Subplot(2, 1, 1)
	for i := 0; i < npts; i++ {
		Sx[i], Sy[i], _, Sxy[i] = o.Stress(t, []float64{d[i], 0}) // y=0
	}
	plt.Plot(d, Sx, "color='r',label='$\\sigma_x$ @ $y=0$'")
	plt.Plot(d, Sy, "color='g',label='$\\sigma_y$ @ $y=0$'")
	plt.Plot(d, Sxy, "color='b',label='$\\sigma_{xy}$ @ $y=0$'")
	plt.Gll("$x$", "stresses", "")

	plt.Subplot(2, 1, 2)
	for i := 0; i < npts; i++ {
		Sx[i], Sy[i], _, Sxy[i] = o.Stress(t, []float64{0, d[i]}) // x=0
	}
	plt.Plot(Sx, d, "color='r',label='$\\sigma_x$ @ $x=0$'")
	plt.Plot(Sy, d, "color='g',label='$\\sigma_y$ @ $x=0$'")
	plt.Plot(Sxy, d, "color='b',label='$\\sigma_{xy}$ @ $x=0$'")
	plt.Gll("stresses", "$y$", "")
}

func Test_functions03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("functions03")

	eps := 1e-2
	f := func(x float64) float64 { return Sabs(x, eps) }
	ff := func(x float64) float64 { return SabsD1(x, eps) }

	np := 401
	//x  := utl.LinSpace(-5e5, 5e5, np)
	//x  := utl.LinSpace(-5e2, 5e2, np)
	x := utl.LinSpace(-5e1, 5e1, np)
	Y := make([]float64, np)
	y := make([]float64, np)
	g := make([]float64, np)
	h := make([]float64, np)
	tolg, tolh := 1e-6, 1e-5
	with_err := false
	for i := 0; i < np; i++ {
		Y[i] = math.Abs(x[i])
		y[i] = Sabs(x[i], eps)
		g[i] = SabsD1(x[i], eps)
		h[i] = SabsD2(x[i], eps)
		gnum := numderiv(f, x[i])
		hnum := numderiv(ff, x[i])
		errg := math.Abs(g[i] - gnum)
		errh := math.Abs(h[i] - hnum)
		clrg, clrh := "[1;32m", "[1;32m"
		if errg > tolg {
			clrg, with_err = "[1;31m", true
		}
		if errh > tolh {
			clrh, with_err = "[1;31m", true
		}
		io.Pf("errg = %s%23.15e   errh = %s%23.15e[0m\n", clrg, errg, clrh, errh)
	}

	if with_err {
		chk.Panic("errors found")
	}

	if false {
		//if true {
		plt.Subplot(3, 1, 1)
		plt.Plot(x, y, "'k--', label='abs'")
		plt.Plot(x, y, "'b-', label='sabs'")
		plt.Gll("x", "y", "")

		plt.Subplot(3, 1, 2)
		plt.Plot(x, g, "'b-', label='sabs'")
		plt.Gll("x", "dy/dx", "")

		plt.Subplot(3, 1, 3)
		plt.Plot(x, h, "'b-', label='sabs'")
		plt.Gll("x", "d2y/dx2", "")

		plt.Show()
	}
}

// PlotYxe plots the function y(x) implemented by Cb_yxe
func PlotYxe(ffcn Cb_yxe, dirout, fname string, xsol, xa, xb float64, np int, xsolLbl, args string, save, show bool, extra func()) (err error) {
	if !save && !show {
		return
	}
	x := utl.LinSpace(xa, xb, np)
	y := make([]float64, np)
	for i := 0; i < np; i++ {
		y[i], err = ffcn(x[i])
		if err != nil {
			return
		}
	}
	var ysol float64
	ysol, err = ffcn(xsol)
	if err != nil {
		return
	}
	plt.Cross("")
	plt.Plot(x, y, args)
	plt.PlotOne(xsol, ysol, io.Sf("'ro', label='%s'", xsolLbl))
	if extra != nil {
		extra()
	}
	plt.Gll("x", "y(x)", "")
	if save {
		os.MkdirAll(dirout, 0777)
		plt.Save(dirout + "/" + fname)
	}
	if show {
		plt.Show()
	}
	return
}

func main() {

	// filename
	filename, fnkey := io.ArgToFilename(0, "sg111", ".sim", true)

	// results
	out.Start(filename, 0, 0)
	out.Define("tip", out.N{1})
	out.LoadResults(nil)

	// plot FEM results
	out.Plot("t", "uy", "tip", plt.Fmt{C: "r", Ls: "None", M: ".", L: "gofem"}, -1)

	// analytical solution
	tAna := utl.LinSpace(0, 5, 101)
	uyAna := make([]float64, len(tAna))
	for i, t := range tAna {
		uyAna[i] = solution_uy(t, 1.0)
	}

	// save
	plt.SetForPng(0.8, 400, 200)
	out.Draw("/tmp", fnkey+".png", false, func(i, j, n int) {
		plt.Plot(tAna, uyAna, "'g-', clip_on=0, label='analytical'")
	})
}

// PlotDerivs plots derivatives of basis functions in I
// option =  0 : use CalcBasisAndDerivs
//           1 : use NumericalDeriv
func (o *Bspline) PlotDerivs(args string, npts, option int) {
	nmks := 10
	tt := utl.LinSpace(o.tmin, o.tmax, npts)
	I := utl.IntRange(o.NumBasis())
	f := make([]float64, len(tt))
	lbls := []string{"N\\&dN", "numD"}
	var cmd string
	for _, i := range I {
		for j, t := range tt {
			switch option {
			case 0:
				o.CalcBasisAndDerivs(t)
				f[j] = o.GetDeriv(i)
			case 1:
				f[j] = o.NumericalDeriv(t, i)
			}
		}
		if strings.Contains(args, "marker") {
			cmd = io.Sf("label=r'%s:%d', color=GetClr(%d, 2) %s", lbls[option], i, i, args)
		} else {
			cmd = io.Sf("label=r'%s:%d', marker=(None if %d %%2 == 0 else GetMrk(%d/2,1)), markevery=(%d-1)/%d, clip_on=0, color=GetClr(%d, 2) %s", lbls[option], i, i, i, npts, nmks, i, args)
		}
		plt.Plot(tt, f, cmd)
	}
	plt.Gll("$t$", io.Sf(`$\frac{\mathrm{d}N_{i,%d}}{\mathrm{d}t}$`, o.p), io.Sf("leg=1, leg_out=1, leg_ncol=%d, leg_hlen=1.5, leg_fsz=7", o.NumBasis()))
	o.plt_ticks_spans()
}

// CheckDerivT checks derivatives w.r.t to t for fixed coordinates x
func CheckDerivT(tst *testing.T, o Func, t0, tf float64, xcte []float64, np int, tskip []float64, sktol, dtol, dtol2 float64, ver bool) {
	t := utl.LinSpace(t0, tf, np)
	for i := 0; i < np; i++ {
		g := o.G(t[i], xcte)
		h := o.H(t[i], xcte)
		skip := false
		for _, val := range tskip {
			if math.Abs(val-t[i]) < sktol {
				skip = true
				break
			}
		}
		if skip {
			continue
		}
		dnum := num.DerivCen(func(t float64, args ...interface{}) (res float64) {
			return o.F(t, xcte)
		}, t[i])
		chk.AnaNum(tst, io.Sf("G(%10f)", t[i]), dtol, g, dnum, ver)
		dnum2 := num.DerivCen(func(t float64, args ...interface{}) (res float64) {
			return o.G(t, xcte)
		}, t[i])
		chk.AnaNum(tst, io.Sf("H(%10f)", t[i]), dtol2, h, dnum2, ver)
	}
}

func Test_mtdeb01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("mtdeb01. Deb's mutation")

	var ops OpsData
	ops.SetDefault()
	ops.Pm = 1.0
	ops.Xrange = [][]float64{{-3, 3}, {-4, 4}}
	ops.EnfRange = true

	rnd.Init(0)

	A := []float64{-1, 1}
	io.Pforan("before: A = %v\n", A)
	FltMutationDeb(A, 10, &ops)
	io.Pforan("after:  A = %v\n", A)

	ha0 := rnd.Histogram{Stations: utl.LinSpace(-3, 3, 11)}

	nsamples := 1000
	aa := make([]float64, len(A))
	a0s := make([]float64, nsamples)
	for _, t := range []int{0, 50, 100} {
		for i := 0; i < nsamples; i++ {
			copy(aa, A)
			FltMutationDeb(aa, t, &ops)
			a0s[i] = aa[0]
		}
		ha0.Count(a0s, true)
		io.Pf("\ntime = %d\n", t)
		io.Pf("%s", rnd.TextHist(ha0.GenLabels("%.1f"), ha0.Counts, 60))
	}
}

// PlotFltOva plots flt-ova points
func (o *Optimiser) PlotFltOva(sols0 []*Solution, iFlt, iOva int, ovaMult float64, pp *PlotParams) {
	if pp.YfuncX != nil {
		X := utl.LinSpace(o.FltMin[iFlt], o.FltMax[iFlt], pp.NptsYfX)
		Y := make([]float64, pp.NptsYfX)
		for i := 0; i < pp.NptsYfX; i++ {
			Y[i] = pp.YfuncX(X[i])
		}
		plt.Plot(X, Y, pp.FmtYfX.GetArgs(""))
	}
	if sols0 != nil {
		o.PlotAddFltOva(iFlt, iOva, sols0, ovaMult, &pp.FmtSols0)
	}
	o.PlotAddFltOva(iFlt, iOva, o.Solutions, ovaMult, &pp.FmtSols)
	best, _ := GetBestFeasible(o, iOva)
	if best != nil {
		plt.PlotOne(best.Flt[iFlt], best.Ova[iOva]*ovaMult, pp.FmtBest.GetArgs(""))
	}
	if pp.Extra != nil {
		pp.Extra()
	}
	if pp.AxEqual {
		plt.Equal()
	}
	plt.Gll(io.Sf("$x_{%d}$", iFlt), io.Sf("$f_{%d}$", iOva), "leg_out=1, leg_ncol=4, leg_hlen=1.5")
	plt.SaveD(pp.DirOut, pp.FnKey+pp.FnExt)
}

func main() {

	// filename
	filename, fnkey := io.ArgToFilename(0, "sg111", ".sim", true)

	// fem
	if !fem.Start(filename, false, false, false) {
		io.PfRed("Start failed\n")
		return
	}
	dom, sum, ok := fem.AllocSetAndInit(0, true, true)
	if !ok {
		io.PfRed("AllocSetAndInit failed\n")
		return
	}

	// selected node and dof index
	nidx := 1
	didx := 1

	// gofem
	ntout := len(sum.OutTimes)
	uy := make([]float64, ntout)
	for tidx, _ := range sum.OutTimes {

		// read results from file
		if !dom.In(sum, tidx, true) {
			io.PfRed("plot_spo751: cannot read solution\n")
			return
		}

		// collect results for load versus time plot
		nod := dom.Nodes[nidx]
		eq := nod.Dofs[didx].Eq
		uy[tidx] = dom.Sol.Y[eq]

		// check
		if math.Abs(dom.Sol.T-sum.OutTimes[tidx]) > 1e-14 {
			io.PfRed("output times do not match time in solution array\n")
			return
		}
	}

	// plot fem results
	plt.SetForPng(0.8, 400, 200)
	plt.Plot(sum.OutTimes, uy, "'ro-', clip_on=0, label='gofem'")

	// analytical solution
	tAna := utl.LinSpace(0, 5, 101)
	uyAna := make([]float64, len(tAna))
	for i, t := range tAna {
		uyAna[i] = solution_uy(t, 1.0)
	}
	plt.Plot(tAna, uyAna, "'g-', clip_on=0, label='analytical'")

	// save
	plt.Gll("$t$", "$u_y$", "")
	plt.SaveD("/tmp", fnkey+".png")
}

func main() {

	// define function and derivative function
	y_fcn := func(x float64) float64 { return math.Sin(x) }
	dydx_fcn := func(x float64) float64 { return math.Cos(x) }
	d2ydx2_fcn := func(x float64) float64 { return -math.Sin(x) }

	// run test for 11 points
	X := utl.LinSpace(0, 2*math.Pi, 11)
	io.Pf("          %8s %23s %23s %23s\n", "x", "analytical", "numerical", "error")
	for _, x := range X {

		// analytical derivatives
		dydx_ana := dydx_fcn(x)
		d2ydx2_ana := d2ydx2_fcn(x)

		// numerical derivative: dydx
		dydx_num, _ := num.DerivCentral(func(t float64, args ...interface{}) float64 {
			return y_fcn(t)
		}, x, 1e-3)

		// numerical derivative d2ydx2
		d2ydx2_num, _ := num.DerivCentral(func(t float64, args ...interface{}) float64 {
			return dydx_fcn(t)
		}, x, 1e-3)

		// check
		chk.PrintAnaNum(io.Sf("dy/dx   @ %.6f", x), 1e-10, dydx_ana, dydx_num, true)
		chk.PrintAnaNum(io.Sf("d²y/dx² @ %.6f", x), 1e-10, d2ydx2_ana, d2ydx2_num, true)
	}

	// generate 101 points for plotting
	X = utl.LinSpace(0, 2*math.Pi, 101)
	Y := make([]float64, len(X))
	for i, x := range X {
		Y[i] = y_fcn(x)
	}

	// plot
	plt.SetForPng(0.75, 300, 150)
	plt.Plot(X, Y, "'b.-', clip_on=0, markevery=10, label='y(x)=sin(x)'")
	plt.Gll("x", "y", "")
	plt.SaveD("/tmp/gosl", "num_deriv01.png")
}

// CheckDerivs compares analytical with numerical derivatives
func (o *Nurbs) CheckDerivs(tst *testing.T, nn int, tol float64, ver bool) {
	dana := make([]float64, 2)
	dnum := make([]float64, 2)
	for _, u := range utl.LinSpace(o.b[0].tmin, o.b[0].tmax, nn) {
		for _, v := range utl.LinSpace(o.b[1].tmin, o.b[1].tmax, nn) {
			uu := []float64{u, v}
			o.CalcBasisAndDerivs(uu)
			for i := 0; i < o.n[0]; i++ {
				for j := 0; j < o.n[1]; j++ {
					l := i + j*o.n[0]
					o.GetDerivL(dana, l)
					o.NumericalDeriv(dnum, uu, l)
					chk.AnaNum(tst, io.Sf("dR[%d][%d][0](%g,%g)", i, j, uu[0], uu[1]), tol, dana[0], dnum[0], ver)
					chk.AnaNum(tst, io.Sf("dR[%d][%d][1](%g,%g)", i, j, uu[0], uu[1]), tol, dana[1], dnum[1], ver)
				}
			}
		}
	}
}

func main() {

	// define function and derivative function
	y_fcn := func(x float64) float64 { return math.Sin(x) }
	dydx_fcn := func(x float64) float64 { return math.Cos(x) }
	d2ydx2_fcn := func(x float64) float64 { return -math.Sin(x) }

	// run test for 11 points
	X := utl.LinSpace(0, 2*math.Pi, 11)
	for _, x := range X {

		// analytical derivatives
		dydx_ana := dydx_fcn(x)
		d2ydx2_ana := d2ydx2_fcn(x)

		// numerical derivative: dydx
		dydx_num, _ := num.DerivCentral(func(t float64, args ...interface{}) float64 {
			return y_fcn(t)
		}, x, 1e-3)

		// numerical derivative d2ydx2
		d2ydx2_num, _ := num.DerivCentral(func(t float64, args ...interface{}) float64 {
			return dydx_fcn(t)
		}, x, 1e-3)

		// check
		chk.PrintAnaNum(io.Sf("dy/dx   @ %.6f", x), 1e-10, dydx_ana, dydx_num, true)
		chk.PrintAnaNum(io.Sf("d²y/dx² @ %.6f", x), 1e-10, d2ydx2_ana, d2ydx2_num, true)
	}

	// generate 101 points
	X = utl.LinSpace(0, 2*math.Pi, 101)
	Y := make([]float64, len(X))
	for i, x := range X {
		Y[i] = y_fcn(x)
	}

	// plot
	plt.Plot(X, Y, "'b.-'")
	plt.Gll("x", "y", "")
	plt.Show()
}

func (o *Bspline) Draw3d(curveArgs, ctrlArgs string, npts int, first bool) {
	t := utl.LinSpace(o.tmin, o.tmax, npts)
	x := make([]float64, npts)
	y := make([]float64, npts)
	z := make([]float64, npts)
	for i, t := range t {
		C := o.Point(t, 0)
		x[i], y[i], z[i] = C[0], C[1], C[2]
	}
	plt.Plot3dLine(x, y, z, first, "")
}

// PlotPdf plots PDF
func (o VarData) PlotPdf(np int, args string) {
	X := utl.LinSpace(o.Min, o.Max, np)
	Y := make([]float64, np)
	for i := 0; i < np; i++ {
		Y[i] = o.Distr.Pdf(X[i])
	}
	if args == "" {
		args = "'b-'"
	}
	plt.Plot(X, Y, args)
	plt.Gll("$x$", "$f(x)$", "")
}

func (o PressCylin) CalcStresses(Pvals []float64, nr int) (R []float64, Sr, St [][]float64) {
	R = utl.LinSpace(o.a, o.b, nr)
	np := len(Pvals)
	Sr = la.MatAlloc(np, nr)
	St = la.MatAlloc(np, nr)
	for i, P := range Pvals {
		c := o.Calc_c(P)
		for j := 0; j < nr; j++ {
			Sr[i][j], St[i][j] = o.Stresses(c, R[j])
		}
	}
	return
}

func Test_2dinteg04(tst *testing.T) {

	//verbose()
	chk.PrintTitle("2dinteg04. ∫∫y・sin(x)")

	x := utl.LinSpace(0, math.Pi/2.0, 11)
	y := utl.LinSpace(0, 1, 11)
	m, n := len(x), len(y)
	f := la.MatAlloc(m, n)
	for i := 0; i < m; i++ {
		for j := 0; j < n; j++ {
			f[i][j] = math.Sin(x[i]) * y[j]
		}
	}
	dx, dy := x[1]-x[0], y[1]-y[0]
	Vt := Trapz2D(dx, dy, f)
	Vs := Simps2D(dx, dy, f)
	io.Pforan("Vt = %v\n", Vt)
	io.Pforan("Vs = %v\n", Vs)
	chk.Scalar(tst, "Vt", 0.00103, Vt, 0.5)
	chk.Scalar(tst, "Vs", 1e-5, Vs, 0.5)
}

func Test_2dinteg03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("2dinteg03. ∫∫(1+8xy)dxdy")

	x := utl.LinSpace(0, 3, 5)
	y := utl.LinSpace(1, 2, 5)
	m, n := len(x), len(y)
	f := la.MatAlloc(m, n)
	for i := 0; i < m; i++ {
		for j := 0; j < n; j++ {
			f[i][j] = 1.0 + 8.0*x[i]*y[j]
		}
	}
	dx, dy := x[1]-x[0], y[1]-y[0]
	Vt := Trapz2D(dx, dy, f)
	Vs := Simps2D(dx, dy, f)
	io.Pforan("Vt = %v\n", Vt)
	io.Pforan("Vs = %v\n", Vs)
	chk.Scalar(tst, "Vt", 1e-13, Vt, 57.0)
	chk.Scalar(tst, "Vs", 1e-13, Vs, 57.0)
}