func Test_ind02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("ind02. copy into")

	rnd.Init(0)

	nbases := 1
	A := get_individual(0, nbases)
	B := get_individual(1, nbases)

	fmts := map[string][]string{
		"int": {"%2d", "%4d", "%5d"}, // ints
		"flt": {"%6g", "%6g", "%5g"}, // floats
		"str": {"%4s", "%2s", "%2s"}, // strings
		"key": {"%3x", "%3x", "%3x"}, // keys
		"byt": {"%4s", "%4s", "%4s"}, // bytes
		"fun": {"%3s", "%3s", "%3s"}, // funcs
	}
	io.Pfpink("A = %v\n", A.Output(fmts, false))
	io.Pfcyan("B = %v\n", B.Output(fmts, false))

	var ops OpsData
	ops.SetDefault()
	ops.Pc = 1.0
	ops.Cuts = []int{1, 2}
	ops.Xrange = [][]float64{{0, 1}, {-20, 20}, {-300, 300}}

	a := A.GetCopy()
	b := A.GetCopy()
	IndCrossover(a, b, A, B, 0, &ops)

	io.Pforan("a = %v\n", a.Output(fmts, false))
	io.Pfblue2("b = %v\n", b.Output(fmts, false))

	chk.Ints(tst, "a.Ints   ", a.Ints, []int{1, -20, 300})
	chk.Ints(tst, "b.Ints   ", b.Ints, []int{-1, 20, -300})
	chk.Strings(tst, "a.Strings", a.Strings, []string{"abc", "Y", "c"})
	chk.Strings(tst, "b.Strings", b.Strings, []string{"X", "b", "Z"})
	// TODO: add other tests here
	io.Pf("\n")

	x := get_individual(0, nbases)
	x.Ovas = []float64{0, 0}
	x.Oors = []float64{0, 0, 0}
	io.Pfblue2("x = %v\n", x.Output(fmts, false))
	B.CopyInto(x)

	chk.Scalar(tst, "ova0", 1e-17, x.Ovas[0], 200)
	chk.Scalar(tst, "ova1", 1e-17, x.Ovas[1], 100)
	chk.Scalar(tst, "oor0", 1e-17, x.Oors[0], 15)
	chk.Scalar(tst, "oor1", 1e-17, x.Oors[1], 25)
	chk.Scalar(tst, "oor2", 1e-17, x.Oors[2], 35)

	io.Pforan("x = %v\n", x.Output(fmts, false))
	chk.String(tst, x.Output(fmts, false), B.Output(fmts, false))
}

func Test_pareto01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("pareto01. compare vectors: Pareto-optimal")

	u := []float64{1, 2, 3, 4, 5, 6}
	v := []float64{1, 2, 3, 4, 5, 6}
	io.Pforan("u = %v\n", u)
	io.Pfblue2("v = %v\n", v)
	u_dominates, v_dominates := DblsParetoMin(u, v)
	io.Pfpink("u_dominates = %v\n", u_dominates)
	io.Pfpink("v_dominates = %v\n", v_dominates)
	if u_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if v_dominates {
		tst.Errorf("test failed\n")
		return
	}

	v = []float64{1, 1.8, 3, 4, 5, 6}
	io.Pforan("\nu = %v\n", u)
	io.Pfblue2("v = %v\n", v)
	u_dominates, v_dominates = DblsParetoMin(u, v)
	io.Pfpink("u_dominates = %v\n", u_dominates)
	io.Pfpink("v_dominates = %v\n", v_dominates)
	if u_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if !v_dominates {
		tst.Errorf("test failed\n")
		return
	}

	v = []float64{1, 2.1, 3, 4, 5, 6}
	io.Pforan("\nu = %v\n", u)
	io.Pfblue2("v = %v\n", v)
	u_dominates, v_dominates = DblsParetoMin(u, v)
	io.Pfpink("u_dominates = %v\n", u_dominates)
	io.Pfpink("v_dominates = %v\n", v_dominates)
	if !u_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if v_dominates {
		tst.Errorf("test failed\n")
		return
	}
}

func main() {

	// input data
	matOld := "matOld.mat"
	matNew := "matNew.mat"
	convSymb := true

	// parse flags
	flag.Parse()
	if len(flag.Args()) > 0 {
		matOld = flag.Arg(0)
	}
	if len(flag.Args()) > 1 {
		matNew = flag.Arg(1)
	}
	if len(flag.Args()) > 2 {
		convSymb = io.Atob(flag.Arg(2))
	}

	// print input data
	io.Pf("\nInput data\n")
	io.Pf("==========\n")
	io.Pf("  matOld   = %30s // old material filename\n", matOld)
	io.Pf("  matNew   = %30s // new material filename\n", matNew)
	io.Pf("  convSymb = %30v // do convert symbols\n", convSymb)
	io.Pf("\n")

	// convert old => new
	inp.MatfileOld2New("", matNew, matOld, convSymb)
	io.Pf("conversion successful\n")
	io.Pfblue2("file <matNew.mat> created\n")
}

func Test_GOflt01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("GOflt01. float64")

	Init(1234)

	xmin := 10.0
	xmax := 20.0
	vals := make([]float64, NSAMPLES)

	// using Float64
	t0 := time.Now()
	for i := 0; i < NSAMPLES; i++ {
		vals[i] = Float64(xmin, xmax)
	}
	io.Pforan("time elapsed = %v\n", time.Now().Sub(t0))

	hist := Histogram{Stations: []float64{10, 12.5, 15, 17.5, 20}}
	hist.Count(vals, true)
	io.Pfpink(TextHist(hist.GenLabels("%4g"), hist.Counts, 60))

	// using Float64s
	t0 = time.Now()
	Float64s(vals, xmin, xmax)
	io.Pforan("time elapsed = %v\n", time.Now().Sub(t0))

	hist.Count(vals, true)
	io.Pfblue2(TextHist(hist.GenLabels("%4g"), hist.Counts, 60))
}

func main() {

	// default input data
	fn := "nurbs01.msh"
	ctrl := true
	ids := true
	npts := 41

	// parse flags
	flag.Parse()
	if len(flag.Args()) > 0 {
		fn = flag.Arg(0)
	}
	if len(flag.Args()) > 1 {
		ctrl = io.Atob(flag.Arg(1))
	}
	if len(flag.Args()) > 2 {
		ids = io.Atob(flag.Arg(2))
	}
	if len(flag.Args()) > 3 {
		npts = io.Atoi(flag.Arg(3))
	}

	// print input data
	io.Pforan("Input data\n")
	io.Pforan("==========\n")
	io.Pfblue2("  fn   = %v\n", fn)
	io.Pfblue2("  ctrl = %v\n", ctrl)
	io.Pfblue2("  ids  = %v\n", ids)
	io.Pfblue2("  npts = %v\n", npts)

	// load nurbss
	fnk := io.FnKey(fn)
	B := gm.ReadMsh(fnk)

	// plot
	plt.SetForEps(0.75, 500)
	for _, b := range B {
		if ctrl {
			b.DrawCtrl2d(ids, "", "")
		}
		b.DrawElems2d(npts, ids, "", "")
	}
	plt.Equal()
	plt.Save(fnk + ".eps")
}

func Test_cxint01(tst *testing.T) {

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

	var ops OpsData
	ops.SetDefault()
	ops.Pc = 1
	ops.Ncuts = 1

	A := []int{1, 2}
	B := []int{-1, -2}
	a := make([]int, len(A))
	b := make([]int, len(A))
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2d\n", A)
	io.PfRed("B = %2d\n", B)
	io.Pfcyan("a = %2d\n", a)
	io.Pfblue2("b = %2d\n", b)
	chk.Ints(tst, "a", a, []int{1, -2})
	chk.Ints(tst, "b", b, []int{-1, 2})
	io.Pf("\n")

	A = []int{1, 2, 3, 4, 5, 6, 7, 8}
	B = []int{-1, -2, -3, -4, -5, -6, -7, -8}
	a = make([]int, len(A))
	b = make([]int, len(A))
	ops.Cuts = []int{1, 3}
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2v\n", A)
	io.PfRed("B = %2v\n", B)
	io.Pfcyan("a = %2v\n", a)
	io.Pfblue2("b = %2v\n", b)
	chk.Ints(tst, "a", a, []int{1, -2, -3, 4, 5, 6, 7, 8})
	chk.Ints(tst, "b", b, []int{-1, 2, 3, -4, -5, -6, -7, -8})

	ops.Cuts = []int{5, 7}
	IntCrossover(a, b, A, B, 0, &ops)
	io.Pfred("A = %2v\n", A)
	io.PfRed("B = %2v\n", B)
	io.Pfcyan("a = %2v\n", a)
	io.Pfblue2("b = %2v\n", b)
	chk.Ints(tst, "a", a, []int{1, 2, 3, 4, 5, -6, -7, 8})
	chk.Ints(tst, "b", b, []int{-1, -2, -3, -4, -5, 6, 7, -8})
}

func Test_split01(tst *testing.T) {

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

	r := DblSplit(" 1e4 1 3   8   88   ")
	io.Pfblue2("r = %v\n", r)
	chk.Vector(tst, "r", 1e-16, r, []float64{1e4, 1, 3, 8, 88})
}

func Test_copy01(tst *testing.T) {

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

	v := []float64{1, 2, 3, 4, 4, 5, 5, 6, 6, 6}
	w := DblCopy(v)
	io.Pfblue2("v = %v\n", v)
	chk.Vector(tst, "w==v", 1e-16, w, v)
}

func Test_fileio01(tst *testing.T) {

	chk.PrintTitle("fileio01")

	// start
	if !Start("data/bh16.sim", true, chk.Verbose) {
		tst.Errorf("test failed\n")
	}
	defer End()

	// domain A
	distr := false
	domA := NewDomain(Global.Sim.Regions[0], distr)
	if domA == nil {
		tst.Errorf("test failed\n")
	}
	if !domA.SetStage(0, Global.Sim.Stages[0], distr) {
		tst.Errorf("test failed\n")
	}
	for i, _ := range domA.Sol.Y {
		domA.Sol.Y[i] = float64(i)
	}
	io.Pforan("domA.Sol.Y = %v\n", domA.Sol.Y)

	// write file
	tidx := 123
	if !domA.SaveSol(tidx) {
		tst.Errorf("test failed")
		return
	}
	dir, fnk := Global.Dirout, Global.Fnkey
	io.Pfblue2("file %v written\n", out_nod_path(dir, fnk, tidx, Global.Rank))

	// domain B
	domB := NewDomain(Global.Sim.Regions[0], distr)
	if domB == nil {
		tst.Errorf("test failed\n")
	}
	if !domB.SetStage(0, Global.Sim.Stages[0], distr) {
		tst.Errorf("test failed")
	}
	io.Pfpink("domB.Sol.Y (before) = %v\n", domB.Sol.Y)

	// read file
	if !domB.ReadSol(dir, fnk, tidx) {
		tst.Errorf("test failed")
		return
	}
	io.Pfgreen("domB.Sol.Y (after) = %v\n", domB.Sol.Y)

	// check
	chk.Vector(tst, "Y", 1e-17, domA.Sol.Y, domB.Sol.Y)
	chk.Vector(tst, "dy/dt", 1e-17, domA.Sol.Dydt, domB.Sol.Dydt)
	chk.Vector(tst, "d²y/dt²", 1e-17, domA.Sol.D2ydt2, domB.Sol.D2ydt2)
}

func save_file(filename string, buf *bytes.Buffer, verbose bool) (err error) {
	fil, err := os.Create(filename)
	if err != nil {
		return
	}
	defer func() { err = fil.Close() }()
	_, err = fil.Write(buf.Bytes())
	if verbose {
		io.Pfblue2("file <%s> written\n", filename)
	}
	return
}

func Test_eigenp01(tst *testing.T) {

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

	// constants
	tolP := 1e-14        // eigenprojectors
	tolS := 1e-13        // spectral decomposition
	toldP := 1e-9        // derivatives of eigenprojectors
	ver := chk.Verbose   // check P verbose
	verdP := chk.Verbose // check dPda verbose

	// run test
	nd := test_nd
	for idxA := 0; idxA < len(test_nd); idxA++ {
		//for idxA := 10; idxA < 11; idxA++ {
		//for idxA := 11; idxA < 12; idxA++ {
		//for idxA := 12; idxA < 13; idxA++ {

		// tensor and eigenvalues
		A := test_AA[idxA]
		a := M_Alloc2(nd[idxA])
		Ten2Man(a, A)
		io.PfYel("\n\ntst # %d ###################################################################################\n", idxA)
		io.Pfblue2("a = %v\n", a)
		io.Pfblue2("λ = %v\n", test_λ[idxA])

		// check eigenprojectors
		io.Pforan("\neigenprojectors\n")
		λsorted := CheckEigenprojs(a, tolP, tolS, ver)
		io.Pfyel("λsorted = %v\n", λsorted)
		λchk := utl.DblGetSorted(test_λ[idxA])
		chk.Vector(tst, "λchk", 1e-12, λsorted, λchk)

		// check derivatives of eigenprojectors
		io.Pforan("\nderivatives\n")
		CheckEigenprojsDerivs(a, toldP, verdP, EV_ZERO)

	}
}

func Test_mat02(tst *testing.T) {

	chk.PrintTitle("mat02 (conversion)")

	convertsymbols := true
	MatfileOld2New("/tmp/gofem/inp", "new_layers.mat", "data/old_layers.mat", convertsymbols)

	mdb := ReadMat("/tmp/gofem/inp/", "new_layers.mat")
	if mdb == nil {
		tst.Errorf("test failed\n")
		return
	}
	io.Pfblue2("%v\n", mdb)
}

func Test_hyperelast02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("hyperelast02 (linear)")

	E, ν := 1500.0, 0.25
	K := Calc_K_from_Enu(E, ν)
	G := Calc_G_from_Enu(E, ν)
	io.Pforan("K = %v\n", K)
	io.Pforan("G = %v\n", G)

	var m HyperElast1
	m.Init(2, false, []*fun.Prm{
		&fun.Prm{N: "K0", V: K},
		&fun.Prm{N: "G0", V: G},
		&fun.Prm{N: "le", V: 1},
	})
	io.Pforan("m = %+v\n", m)

	ε := []float64{-0.001, -0.002, -0.003}
	σ := make([]float64, 3)
	m.L_update(σ, ε)
	io.Pfblue2("ε = %v\n", ε)
	io.Pfcyan("σ = %v\n", σ)

	D := la.MatAlloc(3, 3)
	m.L_CalcD(D, ε)
	la.PrintMat("D", D, "%14.6f", false)

	tol := 1e-11
	verb := io.Verbose
	var tmp float64
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			dnum := num.DerivCen(func(x float64, args ...interface{}) (res float64) {
				tmp, ε[j] = ε[j], x
				m.L_update(σ, ε)
				res = σ[i]
				ε[j] = tmp
				return
			}, ε[j])
			chk.AnaNum(tst, io.Sf("D%d%d", i, j), tol, D[i][j], dnum, verb)
		}
	}
}

func Test_hyperelast03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("hyperelast03 (nonlinear)")

	var m HyperElast1
	m.Init(2, false, []*fun.Prm{
		&fun.Prm{N: "kap", V: 0.05},
		&fun.Prm{N: "kapb", V: 20.0},
		&fun.Prm{N: "G0", V: 1500},
		&fun.Prm{N: "pr", V: 2.2},
		&fun.Prm{N: "pt", V: 11.0},
	})
	io.Pforan("m = %+v\n", m)

	ε := []float64{-0.001, -0.002, -0.003}
	σ := make([]float64, 3)
	m.L_update(σ, ε)
	io.Pfblue2("ε = %v\n", ε)
	io.Pfcyan("σ = %v\n", σ)

	D := la.MatAlloc(3, 3)
	m.L_CalcD(D, ε)
	la.PrintMat("D", D, "%14.6f", false)

	tol := 1e-7
	verb := io.Verbose
	var tmp float64
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			dnum := num.DerivCen(func(x float64, args ...interface{}) (res float64) {
				tmp, ε[j] = ε[j], x
				m.L_update(σ, ε)
				res = σ[i]
				ε[j] = tmp
				return
			}, ε[j])
			chk.AnaNum(tst, io.Sf("D%d%d", i, j), tol, D[i][j], dnum, verb)
		}
	}
}

func Test_mat01(tst *testing.T) {

	chk.PrintTitle("mat01")

	mdb1 := ReadMat("data", "bh.mat")
	if mdb1 == nil {
		tst.Errorf("test failed\n")
		return
	}
	io.Pforan("bh.mat just read:\n%v\n", mdb1)

	fn := "test_bh.mat"
	io.WriteFileSD("/tmp/gofem/inp", fn, mdb1.String())

	mdb2 := ReadMat("/tmp/gofem/inp/", fn)
	if mdb2 == nil {
		tst.Errorf("test failed\n")
		return
	}
	io.Pfblue2("\n%v\n", mdb2)
}

func Test_linipm03(tst *testing.T) {

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

	t0 := time.Now()
	defer func() { io.Pfblue2("\ntime elapsed = %v\n", time.Now().Sub(t0)) }()

	// read LP
	A, b, c, l, u := ReadLPfortran("data/afiro.dat")
	//A, b, c, l, u := ReadLPfortran("data/adlittle.dat")
	//A, b, c, l, u := ReadLPfortran("data/share1b.dat")

	// check for unbounded variables
	nx := len(c)
	for i := 0; i < nx; i++ {
		if math.Abs(l[i]) > 1e-15 {
			chk.Panic("cannot handle l != 0 yet")
		}
		if math.Abs(u[i]-1e20) > 1e-15 {
			chk.Panic("cannot handle u != ∞ yet")
		}
	}

	// solve LP
	var ipm LinIpm
	defer ipm.Clean()
	ipm.Init(A, b, c, nil)
	err := ipm.Solve(chk.Verbose)
	if err != nil {
		tst.Errorf("ipm failed:\n%v", err)
		return
	}

	// check
	io.Pf("\n")
	bres := make([]float64, len(b))
	la.MatVecMul(bres, 1, A.ToDense(), ipm.X)
	chk.Vector(tst, "A*x=b", 1e-13, bres, b)
}

func main() {

	// catch errors
	defer func() {
		if err := recover(); err != nil {
			io.PfRed("ERROR: %v\n", err)
		}
	}()

	// input data
	matOld := io.ArgToString(0, "matOld.mat")
	matNew := io.ArgToString(1, "matNew.mat")
	convSymb := io.ArgToBool(2, true)
	io.Pf("\n%s\n", io.ArgsTable(
		"old material filename", "matOld", matOld,
		"new material filenamen", "matNew", matNew,
		"do convert symbols", "convSymb", convSymb,
	))

	// convert old => new
	inp.MatfileOld2New("", matNew, matOld, convSymb)
	io.Pf("conversion successful\n")
	io.Pfblue2("file <matNew.mat> created\n")
}

// Run runs FE simulation
func Run() (runisok bool) {

	// plot functions
	if Global.Sim.PlotF != nil && Global.Root {
		Global.Sim.Functions.PlotAll(Global.Sim.PlotF, Global.Dirout, Global.Fnkey)
	}

	// alloc domains
	var domains []*Domain
	for _, reg := range Global.Sim.Regions {
		dom := NewDomain(reg, Global.Distr)
		if dom == nil {
			break
		}
		domains = append(domains, dom)
	}
	if Stop() {
		return
	}

	// make sure to call linear solver clean up routines upon exit
	defer func() {
		for _, d := range domains {
			if !d.InitLSol {
				d.LinSol.Clean()
			}
		}
	}()

	// current time and output time
	t := 0.0
	tout := 0.0
	tidx := 0

	// summary of outputs; e.g. with output times
	cputime := time.Now()
	var sum Summary
	sum.OutTimes = []float64{t}
	defer func() {
		sum.Save()
		if Global.Verbose && !Global.Debug {
			io.Pf("\nfinal t  = %v\n", t)
			io.Pfblue2("cpu time = %v\n", time.Now().Sub(cputime))
		}
	}()

	// loop over stages
	for stgidx, stg := range Global.Sim.Stages {

		// time incrementers
		Dt := stg.Control.DtFunc
		DtOut := stg.Control.DtoFunc
		tf := stg.Control.Tf
		tout = t + DtOut.F(t, nil)

		// set stage
		for _, d := range domains {
			if LogErrCond(!d.SetStage(stgidx, Global.Sim.Stages[stgidx], Global.Distr), "SetStage failed") {
				break
			}
			d.Sol.T = t
			if !d.Out(tidx) {
				break
			}
		}
		if Stop() {
			return
		}
		tidx += 1

		// log models
		mconduct.LogModels()
		mreten.LogModels()
		mporous.LogModels()
		msolid.LogModels()

		// skip stage?
		if stg.Skip {
			continue
		}

		// time loop using Richardson's extrapolation
		// TODO: works with only one domain for now
		if Global.Sim.Solver.RE {
			var re RichardsonExtrap
			re.Init(domains[0], Dt)
			if !re.Run(domains[0], &sum, DtOut, &t, tf, tout, &tidx) {
				return
			}
			continue
		}

		// time loop
		ndiverg := 0 // number of steps diverging
		md := 1.0    // time step multiplier if divergence control is on
		var Δt, Δtout float64
		var lasttimestep bool
		for t < tf {

//.........这里部分代码省略.........

func Test_ind03(tst *testing.T) {

	//verbose()
	chk.PrintTitle("ind03. comparing")

	nbases := 1
	A := get_individual(0, nbases)
	B := get_individual(1, nbases)
	A_dominates, B_dominates := IndCompareDet(A, B)
	io.Pfblue2("A: ovas = %v\n", A.Ovas)
	io.Pfblue2("A: oors = %v\n", A.Oors)
	io.Pfcyan("B: ovas = %v\n", B.Ovas)
	io.Pfcyan("B: oors = %v\n", B.Oors)
	io.Pforan("A_dominates = %v\n", A_dominates)
	io.Pforan("B_dominates = %v\n", B_dominates)
	if !A_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if B_dominates {
		tst.Errorf("test failed\n")
		return
	}

	A.Oors = []float64{0, 0, 0}
	B.Oors = []float64{0, 0, 0}
	A_dominates, B_dominates = IndCompareDet(A, B)
	io.Pfblue2("\nA: ovas = %v\n", A.Ovas)
	io.Pfblue2("A: oors = %v\n", A.Oors)
	io.Pfcyan("B: ovas = %v\n", B.Ovas)
	io.Pfcyan("B: oors = %v\n", B.Oors)
	io.Pforan("A_dominates = %v\n", A_dominates)
	io.Pforan("B_dominates = %v\n", B_dominates)
	if A_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if B_dominates {
		tst.Errorf("test failed\n")
		return
	}

	A.Ovas = []float64{200, 100}
	A_dominates, B_dominates = IndCompareDet(A, B)
	io.Pfblue2("\nA: ovas = %v\n", A.Ovas)
	io.Pfblue2("A: oors = %v\n", A.Oors)
	io.Pfcyan("B: ovas = %v\n", B.Ovas)
	io.Pfcyan("B: oors = %v\n", B.Oors)
	io.Pforan("A_dominates = %v\n", A_dominates)
	io.Pforan("B_dominates = %v\n", B_dominates)
	if A_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if B_dominates {
		tst.Errorf("test failed\n")
		return
	}

	A.Ovas = []float64{200, 99}
	A_dominates, B_dominates = IndCompareDet(A, B)
	io.Pfblue2("\nA: ovas = %v\n", A.Ovas)
	io.Pfblue2("A: oors = %v\n", A.Oors)
	io.Pfcyan("B: ovas = %v\n", B.Ovas)
	io.Pfcyan("B: oors = %v\n", B.Oors)
	io.Pforan("A_dominates = %v\n", A_dominates)
	io.Pforan("B_dominates = %v\n", B_dominates)
	if !A_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if B_dominates {
		tst.Errorf("test failed\n")
		return
	}

	A.Ovas = []float64{200, 100}
	B.Ovas = []float64{199, 100}
	A_dominates, B_dominates = IndCompareDet(A, B)
	io.Pfblue2("\nA: ovas = %v\n", A.Ovas)
	io.Pfblue2("A: oors = %v\n", A.Oors)
	io.Pfcyan("B: ovas = %v\n", B.Ovas)
	io.Pfcyan("B: oors = %v\n", B.Oors)
	io.Pforan("A_dominates = %v\n", A_dominates)
	io.Pforan("B_dominates = %v\n", B_dominates)
	if A_dominates {
		tst.Errorf("test failed\n")
		return
	}
	if !B_dominates {
		tst.Errorf("test failed\n")
		return
	}
}

func Test_bspline03(tst *testing.T) {

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

	//             0 1 2 3 4 5 6 7 8 9 10
	T := []float64{0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5}
	var s Bspline
	s.Init(T, 2)
	s.SetControl([][]float64{{0, 0}, {0.5, 1}, {1, 0}, {1.5, 0}, {2, 1}, {2.5, 1}, {3, 0.5}, {3.5, 0}})

	// analytical derivatives
	s.CalcBasisAndDerivs(3.99)
	io.Pfpink("ana: dNdt(t=3.99, i=5) = %v\n", s.GetDeriv(5))
	io.Pfpink("ana: dNdt(t=3.99, i=6) = %v\n", s.GetDeriv(6))
	io.Pfpink("ana: dNdt(t=3.99, i=7) = %v\n", s.GetDeriv(7))
	s.CalcBasisAndDerivs(4.0)
	io.Pforan("ana: dNdt(t=4.00, i=5) = %v\n", s.GetDeriv(5))
	io.Pforan("ana: dNdt(t=4.00, i=6) = %v\n", s.GetDeriv(6))
	io.Pforan("ana: dNdt(t=4.00, i=7) = %v\n", s.GetDeriv(7))

	// numerical derivatives
	io.Pfcyan("num: dNdt(t=3.99, i=5) = %v\n", s.NumericalDeriv(3.99, 5))
	io.Pfcyan("num: dNdt(t=3.99, i=6) = %v\n", s.NumericalDeriv(3.99, 6))
	io.Pfcyan("num: dNdt(t=3.99, i=7) = %v\n", s.NumericalDeriv(3.99, 7))
	io.Pfblue2("num: dNdt(t=4.00, i=5) = %v\n", s.NumericalDeriv(4.00, 5))
	io.Pfblue2("num: dNdt(t=4.00, i=6) = %v\n", s.NumericalDeriv(4.00, 6))
	io.Pfblue2("num: dNdt(t=4.00, i=7) = %v\n", s.NumericalDeriv(4.00, 7))

	ver := false
	tol := 1e-5
	tt := utl.LinSpace(0, 5, 11)
	numd := make([]float64, s.NumBasis())
	anad := make([]float64, s.NumBasis())
	for _, t := range tt {
		for i := 0; i < s.NumBasis(); i++ {
			s.CalcBasisAndDerivs(t)
			anad[i] = s.GetDeriv(i)
			numd[i] = s.NumericalDeriv(t, i)
			// numerical fails @ 4 [4,5,6]
			if t == 4 {
				numd[4] = anad[4]
				numd[5] = anad[5]
				numd[6] = anad[6]
			}
			chk.PrintAnaNum(io.Sf("i=%d t=%v", i, t), tol, anad[i], numd[i], ver)
		}
		chk.Vector(tst, io.Sf("derivs @ %v", t), tol, numd, anad)
	}

	if chk.Verbose {

		npts := 201
		plt.SetForPng(1.5, 600, 150)
		plt.SplotGap(0, 0.3)

		str0 := ",lw=2"
		str1 := ",ls='none',marker='+',color='cyan',markevery=10"
		str2 := ",ls='none',marker='x',markevery=10"
		str3 := ",ls='none',marker='+',markevery=10"
		str4 := ",ls='none',marker='4',markevery=10"

		plt.Subplot(3, 1, 1)
		s.Draw2d(str0, "", npts, 0) // 0 => CalcBasis
		s.Draw2d(str1, "", npts, 1) // 1 => RecursiveBasis

		plt.Subplot(3, 1, 2)
		s.PlotBasis("", npts, 0)   // 0 => CalcBasis
		s.PlotBasis(str2, npts, 1) // 1 => CalcBasisAndDerivs
		s.PlotBasis(str3, npts, 2) // 2 => RecursiveBasis

		plt.Subplot(3, 1, 3)
		s.PlotDerivs("", npts, 0)   // 0 => CalcBasisAndDerivs
		s.PlotDerivs(str4, npts, 1) // 1 => NumericalDeriv

		plt.SaveD("/tmp/gosl/gm", "bspline03.png")
	}
}