// PlotTwoNurbs plots two NURBS for comparison
func PlotTwoNurbs(dirout, fn string, b, c *Nurbs, npts int, ids bool, extra func()) {
	plt.Reset()
	if io.FnExt(fn) == ".eps" {
		plt.SetForEps(1.5, 500)
	} else {
		plt.SetForPng(1.5, 500, 150)
	}

	plt.Subplot(3, 1, 1)
	b.DrawCtrl2d(ids, "", "")
	b.DrawElems2d(npts, ids, "", "")
	if extra != nil {
		extra()
	}
	plt.Equal()

	plt.Subplot(3, 1, 2)
	c.DrawCtrl2d(ids, "", "")
	c.DrawElems2d(npts, ids, "", "")
	plt.Equal()

	plt.Subplot(3, 1, 3)
	b.DrawElems2d(npts, ids, ", lw=3", "")
	c.DrawElems2d(npts, ids, ", color='red', marker='+', markevery=10", "color='green', size=7, va='bottom'")
	plt.Equal()

	plt.SaveD(dirout, fn)
}

// 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 Test_tri01(tst *testing.T) {

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

	V := [][]float64{
		{0.0, 0.0},
		{1.0, 0.0},
		{1.0, 1.0},
		{0.0, 1.0},
		{0.5, 0.5},
	}

	C := [][]int{
		{0, 1, 4},
		{1, 2, 4},
		{2, 3, 4},
		{3, 0, 4},
	}

	if chk.Verbose {
		plt.SetForPng(1, 300, 150)
		Draw(V, C, nil)
		plt.Equal()
		plt.AxisRange(-0.1, 1.1, -0.1, 1.1)
		plt.Gll("x", "y", "")
		plt.SaveD("/tmp/gosl/tri", "tri01.png")
	}
}

func do_plot_nurbs_refined(b, c *Nurbs) {
	plt.SetForEps(1.5, 400)
	plt.Subplot(3, 1, 1)
	b.DrawCtrl2D(true)
	b.DrawElems2D(21, true, "", "")
	plt.Equal()

	plt.Subplot(3, 1, 2)
	c.DrawCtrl2D(true)
	c.DrawElems2D(21, true, "", "")
	plt.Equal()

	plt.Subplot(3, 1, 3)
	b.DrawElems2D(21, true, ", lw=3", "")
	c.DrawElems2D(21, true, ", color='red', marker='+', markevery=10", "color='magenta', size=8, va='bottom'")
	plt.Equal()
}

func main() {

	// input data
	fn, fnk := io.ArgToFilename(0, "nurbs01", ".msh", true)
	ctrl := io.ArgToBool(1, true)
	ids := io.ArgToBool(2, true)
	useminmax := io.ArgToBool(3, false)
	axisequal := io.ArgToBool(4, true)
	xmin := io.ArgToFloat(5, 0)
	xmax := io.ArgToFloat(6, 0)
	ymin := io.ArgToFloat(7, 0)
	ymax := io.ArgToFloat(8, 0)
	eps := io.ArgToBool(9, false)
	npts := io.ArgToInt(10, 41)

	// print input table
	io.Pf("\n%s\n", io.ArgsTable("INPUT ARGUMENTS",
		"mesh filename", "fn", fn,
		"show control points", "ctrl", ctrl,
		"show ids", "ids", ids,
		"use xmin,xmax,ymin,ymax", "useminmax", useminmax,
		"enforce axis.equal", "axisequal", axisequal,
		"min(x)", "xmin", xmin,
		"max(x)", "xmax", xmax,
		"min(y)", "ymin", ymin,
		"max(y)", "ymax", ymax,
		"generate eps instead of png", "eps", eps,
		"number of divisions", "npts", npts,
	))

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

	// plot
	if eps {
		plt.SetForEps(0.75, 500)
	} else {
		plt.SetForPng(0.75, 500, 150)
	}
	for _, b := range B {
		if ctrl {
			b.DrawCtrl2d(ids, "", "")
		}
		b.DrawElems2d(npts, ids, "", "")
	}
	if axisequal {
		plt.Equal()
	}
	if useminmax {
		plt.AxisRange(xmin, xmax, ymin, ymax)
	}
	ext := ".png"
	if eps {
		ext = ".eps"
	}
	plt.Save(fnk + ext)
}

// PlotTwoVarsContour plots contour for two variables problem. len(x) == 2
//  Input
//   dirout  -- directory to save files
//   fnkey   -- file name key for eps figure
//   x       -- solution. can be <nil>
//   np      -- number of points for contour
//   extra   -- called just before saving figure
//   axequal -- axis.equal
//   vmin    -- min 0 values
//   vmax    -- max 1 values
//   f       -- function to plot filled contour. can be <nil>
//   gs      -- functions to plot contour @ level 0. can be <nil>
func PlotTwoVarsContour(dirout, fnkey string, x []float64, np int, extra func(), axequal bool,
	vmin, vmax []float64, f TwoVarsFunc_t, gs ...TwoVarsFunc_t) {
	if fnkey == "" {
		return
	}
	chk.IntAssert(len(vmin), 2)
	chk.IntAssert(len(vmax), 2)
	V0, V1 := utl.MeshGrid2D(vmin[0], vmax[0], vmin[1], vmax[1], np, np)
	var Zf [][]float64
	var Zg [][][]float64
	if f != nil {
		Zf = la.MatAlloc(np, np)
	}
	if len(gs) > 0 {
		Zg = utl.Deep3alloc(len(gs), np, np)
	}
	xtmp := make([]float64, 2)
	for i := 0; i < np; i++ {
		for j := 0; j < np; j++ {
			xtmp[0], xtmp[1] = V0[i][j], V1[i][j]
			if f != nil {
				Zf[i][j] = f(xtmp)
			}
			for k, g := range gs {
				Zg[k][i][j] = g(xtmp)
			}
		}
	}
	plt.Reset()
	plt.SetForEps(0.8, 350)
	if f != nil {
		cmapidx := 0
		plt.Contour(V0, V1, Zf, io.Sf("fsz=7, cmapidx=%d", cmapidx))
	}
	for k, _ := range gs {
		plt.ContourSimple(V0, V1, Zg[k], false, 8, "zorder=5, levels=[0], colors=['yellow'], linewidths=[2], clip_on=0")
	}
	if x != nil {
		plt.PlotOne(x[0], x[1], "'r*', label='optimum', zorder=10")
	}
	if extra != nil {
		extra()
	}
	if dirout == "" {
		dirout = "."
	}
	plt.Cross("clr='grey'")
	plt.SetXnticks(11)
	plt.SetYnticks(11)
	if axequal {
		plt.Equal()
	}
	plt.AxisRange(vmin[0], vmax[0], vmin[1], vmax[1])
	args := "leg_out='1', leg_ncol=4, leg_hlen=1.5"
	plt.Gll("$x_0$", "$x_1$", args)
	plt.SaveD(dirout, fnkey+".eps")
}

func main() {

	// GA parameters
	C := goga.ReadConfParams("tsp-simple.json")
	rnd.Init(C.Seed)

	// location / coordinates of stations
	locations := [][]float64{
		{60, 200}, {180, 200}, {80, 180}, {140, 180}, {20, 160}, {100, 160}, {200, 160},
		{140, 140}, {40, 120}, {100, 120}, {180, 100}, {60, 80}, {120, 80}, {180, 60},
		{20, 40}, {100, 40}, {200, 40}, {20, 20}, {60, 20}, {160, 20},
	}
	nstations := len(locations)
	C.SetIntOrd(nstations)
	C.CalcDerived()

	// objective value function
	C.OvaOor = func(ind *goga.Individual, idIsland, time int, report *bytes.Buffer) {
		L := locations
		ids := ind.Ints
		dist := 0.0
		for i := 1; i < nstations; i++ {
			a, b := ids[i-1], ids[i]
			dist += math.Sqrt(math.Pow(L[b][0]-L[a][0], 2.0) + math.Pow(L[b][1]-L[a][1], 2.0))
		}
		a, b := ids[nstations-1], ids[0]
		dist += math.Sqrt(math.Pow(L[b][0]-L[a][0], 2.0) + math.Pow(L[b][1]-L[a][1], 2.0))
		ind.Ovas[0] = dist
		return
	}

	// evolver
	nova, noor := 1, 0
	evo := goga.NewEvolver(nova, noor, C)
	evo.Run()

	// results
	io.Pfgreen("best = %v\n", evo.Best.Ints)
	io.Pfgreen("best OVA = %v  (871.117353844847)\n\n", evo.Best.Ovas[0])

	// plot travelling salesman path
	if C.DoPlot {
		plt.SetForEps(1, 300)
		X, Y := make([]float64, nstations), make([]float64, nstations)
		for k, id := range evo.Best.Ints {
			X[k], Y[k] = locations[id][0], locations[id][1]
			plt.PlotOne(X[k], Y[k], "'r.', ms=5, clip_on=0, zorder=20")
			plt.Text(X[k], Y[k], io.Sf("%d", id), "fontsize=7, clip_on=0, zorder=30")
		}
		plt.Plot(X, Y, "'b-', clip_on=0, zorder=10")
		plt.Plot([]float64{X[0], X[nstations-1]}, []float64{Y[0], Y[nstations-1]}, "'b-', clip_on=0, zorder=10")
		plt.Equal()
		plt.AxisRange(10, 210, 10, 210)
		plt.Gll("$x$", "$y$", "")
		plt.SaveD("/tmp/goga", "test_evo04.eps")
	}
}

// PlotFltFlt plots flt-flt contour
// use iFlt==-1 || jFlt==-1 to plot all combinations
func (o *Optimiser) PlotFltFltContour(sols0 []*Solution, iFlt, jFlt, iOva int, pp *PlotParams) {
	best, _ := GetBestFeasible(o, iOva)
	plotAll := iFlt < 0 || jFlt < 0
	plotCommands := func(i, j int) {
		o.PlotContour(i, j, iOva, pp)
		if sols0 != nil {
			o.PlotAddFltFlt(i, j, sols0, &pp.FmtSols0)
		}
		o.PlotAddFltFlt(i, j, o.Solutions, &pp.FmtSols)
		if best != nil {
			plt.PlotOne(best.Flt[i], best.Flt[j], pp.FmtBest.GetArgs(""))
		}
		if pp.Extra != nil {
			pp.Extra()
		}
		if pp.AxEqual {
			plt.Equal()
		}
	}
	if plotAll {
		idx := 1
		ncol := o.Nflt - 1
		for row := 0; row < o.Nflt; row++ {
			idx += row
			for col := row + 1; col < o.Nflt; col++ {
				plt.Subplot(ncol, ncol, idx)
				plt.SplotGap(0.0, 0.0)
				plotCommands(col, row)
				if col > row+1 {
					plt.SetXnticks(0)
					plt.SetYnticks(0)
				} else {
					plt.Gll(io.Sf("$x_{%d}$", col), io.Sf("$x_{%d}$", row), "leg=0")
				}
				idx++
			}
		}
		idx = ncol*(ncol-1) + 1
		plt.Subplot(ncol, ncol, idx)
		plt.AxisOff()
		// TODO: fix formatting of open marker, add star to legend
		plt.DrawLegend([]plt.Fmt{pp.FmtSols0, pp.FmtSols, pp.FmtBest}, 8, "center", false, "")
	} else {
		plotCommands(iFlt, jFlt)
		if pp.Xlabel == "" {
			plt.Gll(io.Sf("$x_{%d}$", iFlt), io.Sf("$x_{%d}$", jFlt), pp.LegPrms)
		} else {
			plt.Gll(pp.Xlabel, pp.Ylabel, pp.LegPrms)
		}
	}
	plt.SaveD(pp.DirOut, pp.FnKey+pp.FnExt)
}

func Test_invs05(tst *testing.T) {

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

	if SAVEPLOT {
		plt.Reset()
		plt.SetForPng(1, 500, 125)
		PlotRosette(1.1, true, true, true, 7)
	}

	addtoplot := func(σa, σb float64, σ []float64) {
		plt.PlotOne(σa, σb, "'ro', ms=5")
		plt.Text(σa, σb, io.Sf("$\\sigma_{123}=(%g,%g,%g)$", σ[0], σ[1], σ[2]), "size=8")
	}

	dotest := func(σ []float64, σacor, σbcor, σccor, θcor, tolσ float64) {
		w := M_w(σ)
		θ2 := math.Asin(w) * 180.0 / (3.0 * math.Pi)
		θ3 := M_θ(σ)
		σa, σb, σc := L2O(σ[0], σ[1], σ[2])
		σ0, σ1, σ2 := O2L(σa, σb, σc)
		σI, σA := make([]float64, 3), []float64{σa, σb, σc}
		la.MatVecMul(σI, 1, O2Lmat(), σA) // σI := L * σA
		io.Pf("σa σb σc = %v %v %v\n", σa, σb, σc)
		io.Pf("w        = %v\n", w)
		io.Pf("θ2, θ3   = %v, %v\n", θ2, θ3)
		chk.Scalar(tst, "σa", 1e-17, σa, σacor)
		chk.Scalar(tst, "σb", 1e-17, σb, σbcor)
		chk.Scalar(tst, "σc", 1e-17, σc, σccor)
		chk.Scalar(tst, "σ0", tolσ, σ0, σ[0])
		chk.Scalar(tst, "σ1", tolσ, σ1, σ[1])
		chk.Scalar(tst, "σ2", tolσ, σ2, σ[2])
		chk.Scalar(tst, "σI0", tolσ, σI[0], σ[0])
		chk.Scalar(tst, "σI1", tolσ, σI[1], σ[1])
		chk.Scalar(tst, "σI2", tolσ, σI[2], σ[2])
		chk.Scalar(tst, "θ2", 1e-6, θ2, θcor)
		chk.Scalar(tst, "θ3", 1e-17, θ3, θ2)
		addtoplot(σa, σb, σ)
	}

	dotest([]float64{-1, 0, 0, 0}, 0, 2.0/SQ6, 1.0/SQ3, 30, 1e-15)
	dotest([]float64{0, -1, 0, 0}, 1.0/SQ2, -1.0/SQ6, 1.0/SQ3, 30, 1e-15)
	dotest([]float64{0, 0, -1, 0}, -1.0/SQ2, -1.0/SQ6, 1.0/SQ3, 30, 1e-15)

	if SAVEPLOT {
		plt.Gll("$\\sigma_a$", "$\\sigma_b$", "")
		plt.Equal()
		plt.SaveD("/tmp/gosl", "fig_invs05.png")
	}
}

// PlotStar plots star with normalised OVAs
func (o *Optimiser) PlotStar() {
	nf := o.Nf
	dθ := 2.0 * math.Pi / float64(nf)
	θ0 := 0.0
	if nf == 3 {
		θ0 = -math.Pi / 6.0
	}
	for _, ρ := range []float64{0.25, 0.5, 0.75, 1.0} {
		plt.Circle(0, 0, ρ, "ec='gray',lw=0.5,zorder=5")
	}
	arrowM, textM := 1.1, 1.15
	for i := 0; i < nf; i++ {
		θ := θ0 + float64(i)*dθ
		xi, yi := 0.0, 0.0
		xf, yf := arrowM*math.Cos(θ), arrowM*math.Sin(θ)
		plt.Arrow(xi, yi, xf, yf, "sc=10,st='->',lw=0.7,zorder=10,clip_on=0")
		plt.PlotOne(xf, yf, "'k+', ms=0")
		xf, yf = textM*math.Cos(θ), textM*math.Sin(θ)
		plt.Text(xf, yf, io.Sf("%d", i), "size=6,zorder=10,clip_on=0")
	}
	X, Y := make([]float64, nf+1), make([]float64, nf+1)
	clr := false
	neg := false
	step := 1
	count := 0
	colors := []string{"m", "orange", "g", "r", "b", "k"}
	var ρ float64
	for i, sol := range o.Solutions {
		if sol.Feasible() && sol.FrontId == 0 && i%step == 0 {
			for j := 0; j < nf; j++ {
				if neg {
					ρ = 1.0 - sol.Ova[j]/(o.RptFmax[j]-o.RptFmin[j])
				} else {
					ρ = sol.Ova[j] / (o.RptFmax[j] - o.RptFmin[j])
				}
				θ := θ0 + float64(j)*dθ
				X[j], Y[j] = ρ*math.Cos(θ), ρ*math.Sin(θ)
			}
			X[nf], Y[nf] = X[0], Y[0]
			if clr {
				j := count % len(colors)
				plt.Plot(X, Y, io.Sf("'k-',color='%s',markersize=3,clip_on=0", colors[j]))
			} else {
				plt.Plot(X, Y, "'r-',marker='.',markersize=3,clip_on=0")
			}
			count++
		}
	}
	plt.Equal()
	plt.AxisOff()
}

func Test_Mw02(tst *testing.T) {

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

	prms := []string{"φ", "Mfix"}
	vals := []float64{32, 0}
	var o NcteM
	o.Init(prms, vals)

	if SAVE_FIG {
		// rosette
		full, ref := false, true
		r := 1.1 * SQ2 * o.M(1) / 3.0
		PlotRosette(r, full, ref, true, 7)

		// NcteM
		npts := 201
		X := make([]float64, npts)
		Y := make([]float64, npts)
		W := utl.LinSpace(-1, 1, npts)
		for i, w := range W {
			θ := math.Asin(w) / 3.0
			r := SQ2 * o.M(w) / 3.0
			X[i] = -r * math.Sin(math.Pi/6.0-θ)
			Y[i] = r * math.Cos(math.Pi/6.0-θ)
			//plt.Text(X[i], Y[i], io.Sf("$\\\\theta=%.2f$", θ*180.0/math.Pi), "size=8, ha='center', color='red'")
			//plt.Text(X[i], Y[i], io.Sf("$w=%.2f$", w), "size=8, ha='center', color='red'")
		}
		plt.Plot(X, Y, "'b-'")

		// MC
		g := func(θ float64) float64 {
			return SQ2 * o.Sinφ / (SQ3*math.Cos(θ) - o.Sinφ*math.Sin(θ))
		}
		io.Pforan("M( 1) = %v\n", SQ2*o.M(1)/3.0)
		io.Pforan("g(30) = %v\n", g(math.Pi/6.0))
		for i, w := range W {
			θ := math.Asin(w) / 3.0
			r := g(θ)
			X[i] = -r * math.Sin(math.Pi/6.0-θ)
			Y[i] = r * math.Cos(math.Pi/6.0-θ)
		}
		plt.Plot(X, Y, "'k-'")

		// save
		plt.Equal()
		plt.SaveD("/tmp/gosl", "mw02.eps")
	}
}

// PlotNurbsDerivs plots derivatives of basis functions la and lb
func PlotNurbsDerivs(dirout, fn string, b *Nurbs, la, lb int) {
	npts := 41
	plt.Reset()
	if io.FnExt(fn) == ".eps" {
		plt.SetForEps(1.5, 500)
	} else {
		plt.SetForPng(1.5, 600, 150)
	}

	plt.Subplot(4, 2, 1)
	t0 := time.Now()
	b.PlotDeriv(la, 0, "", npts, 0) // 0 => CalcBasisAndDerivs
	io.Pfcyan("time elapsed (calcbasis) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(4, 2, 2)
	t0 = time.Now()
	b.PlotDeriv(la, 0, "", npts, 1) // 1 => NumericalDeriv
	io.Pfcyan("time elapsed (numerical) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(4, 2, 3)
	b.PlotDeriv(la, 1, "", npts, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 4)
	b.PlotDeriv(la, 1, "", npts, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.Subplot(4, 2, 5)
	b.PlotDeriv(lb, 0, "", npts, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 6)
	b.PlotDeriv(lb, 0, "", npts, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.Subplot(4, 2, 7)
	b.PlotDeriv(lb, 1, "", npts, 0) // 0 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(4, 2, 8)
	b.PlotDeriv(lb, 1, "", npts, 1) // 0 => NumericalDeriv
	plt.Equal()

	plt.SaveD(dirout, fn)
}

func Test_bezier02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("bezier02. quadratic Bezier. point-distance")

	bez := BezierQuad{
		Q: [][]float64{
			{-1, 1},
			{0.5, -2},
			{2, 4},
		},
	}

	nx, ny := 5, 5
	xx, yy := utl.MeshGrid2D(-1.5, 2.5, -0.5, 4.5, nx, ny)
	//zz := la.MatAlloc(nx, ny)

	// TODO: finish this test

	doplot := false
	if doplot {
		plt.SetForPng(1, 400, 200)
	}

	C := make([]float64, 2)
	for j := 0; j < ny; j++ {
		for i := 0; i < nx; i++ {
			C[0], C[1] = xx[i][j], yy[i][j]
			d := bez.DistPoint(C, doplot)
			io.Pforan("d = %v\n", d)
		}
	}

	np := 21
	T := utl.LinSpace(0, 1, np)
	X := make([]float64, np)
	Y := make([]float64, np)
	for i, t := range T {
		bez.Point(C, t)
		X[i], Y[i] = C[0], C[1]
	}

	if doplot {
		plt.Plot(X, Y, "'b-', label='Bezier'")
		plt.Gll("x", "y", "")
		plt.Equal()
		plt.SaveD("/tmp", "fig_gm_bezier02.png")
	}
}

// PlotNurbs plots a NURBS
func PlotNurbs(dirout, fn string, b *Nurbs, npts int, ids bool, extra func()) {
	plt.Reset()
	if io.FnExt(fn) == ".eps" {
		plt.SetForEps(1.0, 500)
	} else {
		plt.SetForPng(1.0, 500, 150)
	}
	b.DrawCtrl2d(ids, "", "")
	b.DrawElems2d(npts, ids, "", "")
	if extra != nil {
		extra()
	}
	plt.Equal()
	plt.SaveD(dirout, fn)
}

// Draw2d draws bins' grid
func (o *Bins) Draw2d(withtxt bool) {

	// horizontal lines
	x := []float64{o.Xi[0], o.Xi[0] + o.L[0] + o.S}
	y := make([]float64, 2)
	for j := 0; j < o.N[1]+1; j++ {
		y[0] = o.Xi[1] + float64(j)*o.S
		y[1] = y[0]
		plt.Plot(x, y, "'-', color='#4f3677', clip_on=0")
	}

	// vertical lines
	y[0] = o.Xi[1]
	y[1] = o.Xi[1] + o.L[1] + o.S
	for i := 0; i < o.N[0]+1; i++ {
		x[0] = o.Xi[0] + float64(i)*o.S
		x[1] = x[0]
		plt.Plot(x, y, "'k-', color='#4f3677', clip_on=0")
	}

	// plot items
	for _, bin := range o.All {
		if bin == nil {
			continue
		}
		for _, entry := range bin.Entries {
			plt.PlotOne(entry.X[0], entry.X[1], "'r.', clip_on=0")
		}
	}

	// labels
	if withtxt {
		for j := 0; j < o.N[1]; j++ {
			for i := 0; i < o.N[0]; i++ {
				idx := i + j*o.N[0]
				x := o.Xi[0] + float64(i)*o.S + 0.02*o.S
				y := o.Xi[1] + float64(j)*o.S + 0.02*o.S
				plt.Text(x, y, io.Sf("%d", idx), "size=7")
			}
		}
	}

	// setup
	plt.Equal()
	plt.AxisRange(o.Xi[0]-0.1, o.Xf[0]+o.S+0.1, o.Xi[1]-0.1, o.Xf[1]+o.S+0.1)
}

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_halton01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("halton01. Halton Points")

	dim := 2
	npts := 100
	P := HaltonPoints(dim, npts)
	X := P[0]
	Y := P[1]

	if chk.Verbose {
		plt.SetForEps(1, 400)
		plt.Plot(X, Y, "'r.', clip_on=0")
		plt.Equal()
		plt.SaveD("/tmp/gosl", "halton01.eps")
	}
}

func Test_delaunay01(tst *testing.T) {

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

	// points
	X := []float64{0, 1, 1, 0, 0.5}
	Y := []float64{0, 0, 1, 1, 0.5}

	// generate
	V, C, err := Delaunay(X, Y, chk.Verbose)
	if err != nil {
		tst.Errorf("%v\n", err)
		return
	}

	// check
	xout := make([]float64, len(V))
	yout := make([]float64, len(V))
	for i, v := range V {
		io.Pforan("vert %2d : coords = %v\n", i, v)
		xout[i] = v[0]
		yout[i] = v[1]
	}
	chk.Vector(tst, "X", 1e-15, xout, X)
	chk.Vector(tst, "Y", 1e-15, yout, Y)
	for i, c := range C {
		io.Pforan("cell %2d : verts = %v\n", i, c)
	}
	chk.Ints(tst, "verts of cell 0", C[0], []int{3, 0, 4})
	chk.Ints(tst, "verts of cell 1", C[1], []int{4, 1, 2})
	chk.Ints(tst, "verts of cell 2", C[2], []int{1, 4, 0})
	chk.Ints(tst, "verts of cell 3", C[3], []int{4, 2, 3})

	// plot
	if chk.Verbose {
		plt.SetForPng(1, 300, 150)
		Draw(V, C, nil)
		plt.Equal()
		plt.AxisRange(-0.1, 1.1, -0.1, 1.1)
		plt.Gll("x", "y", "")
		plt.SaveD("/tmp/gosl/tri", "delaunay01.png")
	}
}

// PlotNurbsBasis plots basis functions la and lb
func PlotNurbsBasis(dirout, fn string, b *Nurbs, la, lb int) {
	npts := 41
	plt.Reset()
	if io.FnExt(fn) == ".eps" {
		plt.SetForEps(1.5, 500)
	} else {
		plt.SetForPng(1.5, 600, 150)
	}

	plt.Subplot(3, 2, 1)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	t0 := time.Now()
	b.PlotBasis(la, "", 11, 0) // 0 => CalcBasis
	io.Pfcyan("time elapsed (calcbasis) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(3, 2, 2)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 0) // 0 => CalcBasis
	plt.Equal()

	plt.Subplot(3, 2, 3)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(la, "", 11, 1) // 1 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(3, 2, 4)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 1) // 1 => CalcBasisAndDerivs
	plt.Equal()

	plt.Subplot(3, 2, 5)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	t0 = time.Now()
	b.PlotBasis(la, "", 11, 2) // 2 => RecursiveBasis
	io.Pfcyan("time elapsed (recursive) = %v\n", time.Now().Sub(t0))
	plt.Equal()

	plt.Subplot(3, 2, 6)
	b.DrawCtrl2d(false, "", "")
	b.DrawElems2d(npts, false, "", "")
	b.PlotBasis(lb, "", 11, 2) // 2 => RecursiveBasis
	plt.Equal()

	plt.SaveD(dirout, fn)
}

// PlotHc2d plots 2D hypercube
func PlotHc2d(dirout, fnkey string, x [][]int, xrange [][]float64) {
	m := len(x)
	n := len(x[0])
	dx := make([]float64, m)
	for i := 0; i < m; i++ {
		dx[i] = (xrange[i][1] - xrange[i][0]) / float64(n-1)
	}
	X := utl.DblsAlloc(m, n)
	for i := 0; i < m; i++ {
		for j := 0; j < n; j++ {
			X[i][j] = xrange[i][0] + float64(x[i][j]-1)*dx[i]
		}
	}
	plt.SetForEps(0.8, 300)
	plt.Plot(X[0], X[1], "'r.', clip_on=0, zorder=10")
	plt.Equal()
	plt.Gll("$x$", "$y$", "")
	plt.SaveD(dirout, fnkey+".eps")
}