// 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()
	}
}

// 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)
}

// Plot plot results
func Plot(dirout, fn string, res *Results, yfcn Cb_ycorr, xa, xb float64, argsAna, argsNum string, extra func()) {

	// data
	if res == nil {
		return
	}
	ndim := len(res.Y)
	if ndim < 1 {
		return
	}

	// closed-form solution
	var xc []float64
	var Yc [][]float64
	if yfcn != nil {
		np := 101
		dx := (xb - xa) / float64(np-1)
		xc = make([]float64, np)
		Yc = utl.DblsAlloc(np, ndim)
		for i := 0; i < np; i++ {
			xc[i] = xa + dx*float64(i)
			yfcn(Yc[i], xc[i])
		}
	}

	// plot
	if argsAna == "" {
		argsAna = "'y-', lw=6, label='analytical', clip_on=0"
	}
	if argsNum == "" {
		argsNum = "'b-', marker='.', lw=1, clip_on=0"
	}
	for j := 0; j < ndim; j++ {
		plt.Subplot(ndim+1, 1, j+1)
		if yfcn != nil {
			plt.Plot(xc, Yc[j], argsAna)
		}
		plt.Plot(res.X, res.Y[j], argsNum+","+io.Sf("label='%s'", res.Method))
		plt.Gll("$x$", "$y$", "")
	}
	plt.Subplot(ndim+1, 1, ndim+1)
	plt.Plot(res.X, res.Dx, io.Sf("'b-', marker='.', lw=1, clip_on=0, label='%s'", res.Method))
	plt.SetYlog()
	plt.Gll("$x$", "$\\log(\\delta x)$", "")

	// write file
	if extra != nil {
		extra()
	}
	plt.SaveD(dirout, fn)
}

// PlotT plots F, G and H for varying t and fixed coordinates x
func PlotT(o Func, dirout, fname string, t0, tf float64, xcte []float64, np int, args string, withG, withH, save, show bool, extra func()) {
	t := utl.LinSpace(t0, tf, np)
	y := make([]float64, np)
	for i := 0; i < np; i++ {
		y[i] = o.F(t[i], xcte)
	}
	var g, h []float64
	nrow := 1
	if withG {
		g = make([]float64, np)
		for i := 0; i < np; i++ {
			g[i] = o.G(t[i], xcte)
		}
		nrow += 1
	}
	if withH {
		h = make([]float64, np)
		for i := 0; i < np; i++ {
			h[i] = o.H(t[i], xcte)
		}
		nrow += 1
	}
	os.MkdirAll(dirout, 0777)
	if withG || withH {
		plt.Subplot(nrow, 1, 1)
	}
	plt.Plot(t, y, args)
	if extra != nil {
		extra()
	}
	plt.Gll("t", "y", "")
	pidx := 2
	if withG {
		plt.Subplot(nrow, 1, pidx)
		plt.Plot(t, g, args)
		plt.Gll("t", "dy/dt", "")
		pidx += 1
	}
	if withH {
		plt.Subplot(nrow, 1, pidx)
		plt.Plot(t, h, args)
		plt.Gll("t", "d2y/dt2", "")
	}
	if save {
		plt.Save(dirout + "/" + fname)
	}
	if show {
		plt.Show()
	}
}

// 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 (o *Plotter) Plot_i_f(x, y []float64, res []*State, sts [][]float64, last bool) {
	if o.m == nil {
		o.set_empty()
		return
	}
	nr := len(res)
	var y2 []float64
	if o.nsurf > 1 {
		y2 = make([]float64, nr)
	}
	for i := 0; i < nr; i++ {
		ys := o.m.YieldFuncs(res[i])
		y[i] = ys[0]
		if o.nsurf > 1 {
			y2[i] = ys[1]
		}
		x[i] = float64(i)
	}
	lbl := "f " + o.Lbl
	plt.Plot(x, y, io.Sf("'r.', ls='-', clip_on=0, color='%s', marker='%s', label=r'%s'", o.Clr, o.Mrk, lbl))
	if o.nsurf > 1 {
		lbl = "F " + o.Lbl
		plt.Plot(x, y2, io.Sf("'b+', ls=':', lw=2, clip_on=0, color='%s', marker='%s', label=r'%s'", o.Clr, o.Mrk, lbl))
	}
	if last {
		plt.Gll("$i$", "$f,\\;F$", "leg_out=1, leg_ncol=4, leg_hlen=2")
		if lims, ok := o.Lims["i,f"]; ok {
			plt.AxisLims(lims)
		}
	}
}

func (o *Plotter) Plot_i_alp(x, y []float64, res []*State, sts [][]float64, last bool) {
	nr := len(res)
	nα := len(res[0].Alp)
	if nα == 0 {
		o.set_empty()
		return
	}
	yy := la.MatAlloc(nα, nr)
	for i := 0; i < nr; i++ {
		x[i] = float64(i)
		for j := 0; j < nα; j++ {
			yy[j][i] = res[i].Alp[j]
		}
	}
	for j := 0; j < nα; j++ {
		lbl := io.Sf("$\\alpha_%d$ "+o.Lbl, j)
		plt.Plot(x, yy[j], io.Sf("'r-', ls='-', clip_on=0, color='%s', marker='%s', label=r'%s'", o.Clr, o.Mrk, lbl))
	}
	if last {
		plt.Gll("$i$", "$\\alpha_k$", "leg_out=1, leg_ncol=4, leg_hlen=2")
		if lims, ok := o.Lims["i,alp"]; ok {
			plt.AxisLims(lims)
		}
	}
}

func (o *Plotter) Plot_ed_q(x, y []float64, res []*State, sts [][]float64, last bool) {
	nr := len(res)
	if len(sts) != nr {
		return
	}
	k := nr - 1
	for i := 0; i < nr; i++ {
		x[i] = o.Ed[i] * 100.0
		if o.QdivP {
			y[i] = o.Q[i] / o.P[i]
		} else {
			y[i] = o.Q[i]
		}
		if o.Multq {
			y[i] *= fun.Sign(o.W[i])
		}
	}
	plt.Plot(x, y, io.Sf("'r.', ls='%s', clip_on=0, color='%s', marker='%s', label=r'%s'", o.Ls, o.Clr, o.Mrk, o.Lbl))
	plt.PlotOne(x[0], y[0], io.Sf("'bo', clip_on=0, color='%s', marker='%s', ms=%d", o.SpClr, o.SpMrk, o.SpMs))
	plt.PlotOne(x[k], y[k], io.Sf("'bs', clip_on=0, color='%s', marker='%s', ms=%d", o.SpClr, o.EpMrk, o.EpMs))
	if last {
		ylbl := "$q$"
		if o.QdivP {
			ylbl = "$q/p$"
		}
		plt.Gll("$\\varepsilon_d\\;[\\%]$", ylbl, "leg_out=1, leg_ncol=4, leg_hlen=1.5")
		if lims, ok := o.Lims["ed,q"]; ok {
			plt.AxisLims(lims)
		}
	}
}

func (o *Plotter) Plot_Dgam_f(x, y []float64, res []*State, sts [][]float64, last bool) {
	if o.m == nil {
		o.set_empty()
		return
	}
	nr := len(res)
	k := nr - 1
	ys := o.m.YieldFuncs(res[0])
	fc0 := ys[0]
	xmi, xma, ymi, yma := res[0].Dgam, res[0].Dgam, fc0, fc0
	for i := 0; i < nr; i++ {
		x[i] = res[i].Dgam
		ys = o.m.YieldFuncs(res[i])
		y[i] = ys[0]
		xmi = min(xmi, x[i])
		xma = max(xma, x[i])
		ymi = min(ymi, y[i])
		yma = max(yma, y[i])
	}
	//o.DrawRamp(xmi, xma, ymi, yma)
	plt.Plot(x, y, io.Sf("'r.', ls='%s', clip_on=0, color='%s', marker='%s', label=r'%s'", o.Ls, o.Clr, o.Mrk, o.Lbl))
	plt.PlotOne(x[0], y[0], io.Sf("'bo', clip_on=0, color='%s', marker='%s', ms=%d", o.SpClr, o.SpMrk, o.SpMs))
	plt.PlotOne(x[k], y[k], io.Sf("'bs', clip_on=0, color='%s', marker='%s', ms=%d", o.SpClr, o.EpMrk, o.EpMs))
	if last {
		plt.Gll("$\\Delta\\gamma$", "$f$", "")
		if lims, ok := o.Lims["Dgam,f"]; ok {
			plt.AxisLims(lims)
		}
	}
}

// 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
}

// 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()
}

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 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")
}

// PlotEnd ends plot and show figure, if show==true
func PlotEnd(show bool) {
	plt.AxisYrange(0, 1)
	plt.Cross("")
	plt.Gll("$p_c$", "$s_{\\ell}$", "")
	if show {
		plt.Show()
	}
}

func main() {
	Nf := []float64{5, 7, 10, 13, 15, 20}
	Eave := []float64{3.5998e-12, 2.9629e-10, 6.0300e-8, 3.3686e-6, 2.5914e-5, 1.1966e-3}
	plt.SetForEps(0.75, 200)
	plt.Plot(Nf, Eave, "'b-', marker='.', clip_on=0")
	plt.SetYlog()
	plt.Gll("$N_f$", "$E_{ave}$", "")
	plt.SaveD("/tmp/goga", "multierror.eps")
}

// 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")
	}
}

// PlotBasis plots basis function (2D only)
// option =  0 : use CalcBasis
//           1 : use CalcBasisAndDerivs
//           2 : use RecursiveBasis
func (o *Nurbs) PlotBasis(l int, args string, npts, option int) {
	lbls := []string{"CalcBasis function", "CalcBasisAndDerivs function", "RecursiveBasis function"}
	switch o.gnd {
	// curve
	case 1:
		U := make([]float64, npts)
		S := make([]float64, npts)
		du := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		uvec := []float64{0}
		for m := 0; m < npts; m++ {
			U[m] = o.b[0].tmin + float64(m)*du
			uvec[0] = U[m]
			switch option {
			case 0:
				o.CalcBasis(uvec)
				S[m] = o.GetBasisL(l)
			case 1:
				o.CalcBasisAndDerivs(uvec)
				S[m] = o.GetBasisL(l)
			case 2:
				S[m] = o.RecursiveBasis(uvec, l)
			}
		}
		plt.Plot(U, S, args)
		plt.Gll("$u$", io.Sf("$S_%d$", l), "")
	// surface
	case 2:
		xx := la.MatAlloc(npts, npts)
		yy := la.MatAlloc(npts, npts)
		zz := la.MatAlloc(npts, npts)
		du0 := (o.b[0].tmax - o.b[0].tmin) / float64(npts-1)
		du1 := (o.b[1].tmax - o.b[1].tmin) / float64(npts-1)
		for m := 0; m < npts; m++ {
			u0 := o.b[0].tmin + float64(m)*du0
			for n := 0; n < npts; n++ {
				u1 := o.b[1].tmin + float64(n)*du1
				u := []float64{u0, u1}
				x := o.Point(u)
				xx[m][n] = x[0]
				yy[m][n] = x[1]
				switch option {
				case 0:
					o.CalcBasis(u)
					zz[m][n] = o.GetBasisL(l)
				case 1:
					o.CalcBasisAndDerivs(u)
					zz[m][n] = o.GetBasisL(l)
				case 2:
					zz[m][n] = o.RecursiveBasis(u, l)
				}
			}
		}
		plt.Contour(xx, yy, zz, "fsz=7")
	}
	plt.Title(io.Sf("%s @ %d", lbls[option], l), "size=7")
}

func Test_igd01(tst *testing.T) {

	//verbose()
	chk.PrintTitle("igd. igd metric with star equal to trial => igd=0")

	// load star values
	prob := "UF1"
	fStar, err := io.ReadMatrix(io.Sf("./examples/mulobj-cec09/cec09/pf_data/%s.dat", prob))
	if err != nil {
		tst.Errorf("cannot read fStar matrix:\n%v", err)
		return
	}
	npts := len(fStar)

	// optimiser
	var opt Optimiser
	opt.Default()
	opt.Nsol = npts
	opt.Ncpu = 1
	opt.FltMin = []float64{0, 0} // used to store fStar
	opt.FltMax = []float64{1, 1} // used to store fStar
	nf, ng, nh := 2, 0, 0

	// generator (store fStar into Flt)
	gen := func(sols []*Solution, prms *Parameters) {
		for i, sol := range sols {
			sol.Flt[0], sol.Flt[1] = fStar[i][0], fStar[i][1]
		}
	}

	// objective function (copy fStar from Flt into Ova)
	obj := func(f, g, h, x []float64, ξ []int, cpu int) {
		f[0], f[1] = x[0], x[1]
	}

	// initialise optimiser
	opt.Init(gen, nil, obj, nf, ng, nh)

	// compute igd
	igd := StatIgd(&opt, fStar)
	io.Pforan("igd = %v\n", igd)
	chk.Scalar(tst, "igd", 1e-15, igd, 0)

	// plot
	if chk.Verbose {
		fmt := &plt.Fmt{C: "red", M: ".", Ms: 1, Ls: "None", L: "solutions"}
		fS0 := utl.DblsGetColumn(0, fStar)
		fS1 := utl.DblsGetColumn(1, fStar)
		io.Pforan("len(fS0) = %v\n", len(fS0))
		plt.SetForEps(0.75, 300)
		opt.PlotAddOvaOva(0, 1, opt.Solutions, true, fmt)
		plt.Plot(fS0, fS1, io.Sf("'b.', ms=2, label='star(%s)', clip_on=0", prob))
		plt.Gll("$f_0$", "$f_1$", "")
		plt.SaveD("/tmp/goga", "igd01.eps")
	}
}