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
}

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

// run_rootsol_test runs root solution test
//  Note: xguess is the trial solution for Newton's method (not Brent's)
func run_rootsol_test(tst *testing.T, xa, xb, xguess, tolcmp float64, ffcnA Cb_yxe, ffcnB Cb_f, JfcnB Cb_Jd, fname string, save, show bool) (xbrent float64) {

	// Brent
	io.Pfcyan("\n       - - - - - - - using Brent's method - - -- - - - \n")
	var o Brent
	o.Init(ffcnA)
	var err error
	xbrent, err = o.Solve(xa, xb, false)
	if err != nil {
		chk.Panic("%v", err)
	}
	var ybrent float64
	ybrent, err = ffcnA(xbrent)
	if err != nil {
		chk.Panic("%v", err)
	}
	io.Pforan("x      = %v\n", xbrent)
	io.Pforan("f(x)   = %v\n", ybrent)
	io.Pforan("nfeval = %v\n", o.NFeval)
	io.Pforan("nit    = %v\n", o.It)
	if math.Abs(ybrent) > 1e-10 {
		chk.Panic("Brent failed: f(x) = %g > 1e-10\n", ybrent)
	}

	// Newton
	io.Pfcyan("\n       - - - - - - - using Newton's method - - -- - - - \n")
	var p NlSolver
	p.Init(1, ffcnB, nil, JfcnB, true, false, nil)
	xnewt := []float64{xguess}
	var cnd float64
	cnd, err = p.CheckJ(xnewt, 1e-6, true, !chk.Verbose)
	io.Pforan("cond(J) = %v\n", cnd)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	err = p.Solve(xnewt, false)
	if err != nil {
		chk.Panic("%v", err.Error())
	}
	var ynewt float64
	ynewt, err = ffcnA(xnewt[0])
	if err != nil {
		chk.Panic("%v", err)
	}
	io.Pforan("x      = %v\n", xnewt[0])
	io.Pforan("f(x)   = %v\n", ynewt)
	io.Pforan("nfeval = %v\n", p.NFeval)
	io.Pforan("nJeval = %v\n", p.NJeval)
	io.Pforan("nit    = %v\n", p.It)
	if math.Abs(ynewt) > 1e-9 {
		chk.Panic("Newton failed: f(x) = %g > 1e-10\n", ynewt)
	}

	// compare Brent's and Newton's solutions
	PlotYxe(ffcnA, "results", fname, xbrent, xa, xb, 101, "Brent", "'b-'", save, show, func() {
		plt.PlotOne(xnewt[0], ynewt, "'g+', ms=15, label='Newton'")
	})
	chk.Scalar(tst, "xbrent - xnewt", tolcmp, xbrent, xnewt[0])
	return
}

func (o *Plotter) Plot_ed_ev(x, y []float64, res []*State, sts [][]float64, last bool) {
	nr := len(sts)
	k := nr - 1
	for i := 0; i < nr; i++ {
		x[i], y[i] = o.Ed[i]*100.0, o.Ev[i]*100.0
	}
	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("$\\varepsilon_d\\;[\\%]$", "$\\varepsilon_v\\;[\\%]$", "leg_out=1, leg_ncol=4, leg_hlen=1.5")
		if lims, ok := o.Lims["ed,ev"]; ok {
			plt.AxisLims(lims)
		}
	}
}

// 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 draw_truss(dat *FemData, key string, A *goga.Solution, lef, bot, wid, hei float64) (weight, deflection float64) {
	gap := 0.1
	plt.PyCmds(io.Sf(`
from pylab import axes, setp, sca
ax_current = gca()
ax_new = axes([%g, %g, %g, %g], axisbg='#dcdcdc')
setp(ax_new, xticks=[0,720], yticks=[0,360])
axis('equal')
axis('off')
`, lef, bot, wid, hei))
	_, _, weight, deflection, _, _, _ = dat.RunFEM(A.Int, A.Flt, 1, false)
	plt.PyCmds("sca(ax_current)\n")
	plt.PlotOne(weight, deflection, "'g*', zorder=1000, clip_on=0")
	plt.Text(weight, deflection+gap, key, "")
	return
}

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

	//verbose()
	chk.PrintTitle("cubiceq03. y(x) = x³ + c")

	doplot := false
	np := 41
	var X, Y []float64
	if doplot {
		X = utl.LinSpace(-2, 2, np)
		Y = make([]float64, np)
		plt.SetForPng(0.8, 400, 200)
	}

	a, b := 0.0, 0.0
	colors := []string{"red", "green", "blue"}
	for k, c := range []float64{-1, 0, 1} {
		x1, x2, x3, nx := EqCubicSolveReal(a, b, c)
		io.Pforan("\na=%v b=%v c=%v\n", a, b, c)
		io.Pfcyan("nx=%v\n", nx)
		io.Pfcyan("x1=%v x2=%v x3=%v\n", x1, x2, x3)
		chk.IntAssert(nx, 1)
		chk.Scalar(tst, "x1", 1e-17, x1, -c)
		if doplot {
			for i, x := range X {
				Y[i] = x*x*x + a*x*x + b*x + c
			}
			plt.Plot(X, Y, io.Sf("color='%s', label='c=%g'", colors[k], c))
			plt.PlotOne(x1, 0, io.Sf("'ko', color='%s'", colors[k]))
			plt.Cross("")
			plt.Gll("x", "y", "")
		}
	}
	if doplot {
		plt.SaveD("/tmp", "fig_cubiceq03.png")
	}
}

//.........这里部分代码省略.........
	P := make([]float64, nto)
	Ub := make([]float64, nto)
	R := utl.Deep3alloc(len(Psel), nels, nips)
	Sr := utl.Deep3alloc(len(Psel), nels, nips)
	St := utl.Deep3alloc(len(Psel), nels, nips)
	i := 0
	for tidx, t := range sum.OutTimes {

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

		// collect results for load versus displacement plot
		nod := d.Nodes[nidx]
		eq := nod.Dofs[didx].Eq
		P[tidx] = t * Pcen
		Ub[tidx] = d.Sol.Y[eq]

		// stresses
		if isPsel(Psel, P[tidx], tolPsel) {
			for j, ele := range d.ElemIntvars {
				e := ele.(*ElemU)
				ipsdat := e.OutIpsData()
				for k, dat := range ipsdat {
					res := dat.Calc(d.Sol)
					x, y := dat.X[0], dat.X[1]
					sx := res["sx"] * GPa2MPa
					sy := res["sy"] * GPa2MPa
					sxy := res["sxy"] * GPa2MPa / math.Sqrt2
					R[i][j][k], Sr[i][j][k], St[i][j][k], _ = ana.PolarStresses(x, y, sx, sy, sxy)
				}
			}
			i++
		}
	}

	// auxiliary data for plotting stresses
	colors := []string{"r", "m", "g", "k", "y", "c", "r", "m"}
	markers1 := []string{"o", "s", "x", ".", "^", "*", "o", "s"}
	markers2 := []string{"+", "+", "+", "+", "+", "+", "+", "+"}

	// plot load displacements
	plt.SetForEps(0.8, 300)
	if true {
		//if false {
		plt.Plot(Ub_ana, P_ana, "'b-', ms=2, label='solution', clip_on=0")
		plt.Plot(Ub, P, "'r.--', label='fem: outer', clip_on=0")
		plt.Gll("$u_x\\;\\mathrm{[mm]}$", "$P\\;\\mathrm{[MPa]}$", "")
		plt.SaveD("/tmp", io.Sf("gofem_%s_disp.eps", fnkey))
	}

	// plot radial stresses
	if true {
		//if false {
		plt.Reset()
		for i, Pval := range Psel {
			plt.Plot(R_ana, Sr_ana[i], "'b-'")
			for k := 0; k < nips; k++ {
				for j := 0; j < nels; j++ {
					args := io.Sf("'%s%s'", colors[i], markers1[i])
					if k > 1 {
						args = io.Sf("'k%s', ms=10", markers2[i])
					}
					if k == 0 && j == 0 {
						args += io.Sf(", label='P=%g'", Pval)
					}
					plt.PlotOne(R[i][j][k], Sr[i][j][k], args)
				}
			}
		}
		plt.Gll("$r\\;\\mathrm{[mm]}$", "$\\sigma_r\\;\\mathrm{[MPa]}$", "leg_loc='lower right'")
		plt.AxisXrange(a, b)
		plt.SaveD("/tmp", io.Sf("gofem_%s_sr.eps", fnkey))
	}

	// plot tangential stresses
	if true {
		//if false {
		plt.Reset()
		for i, Pval := range Psel {
			plt.Plot(R_ana, St_ana[i], "'b-'")
			for k := 0; k < nips; k++ {
				for j := 0; j < nels; j++ {
					args := io.Sf("'%s%s'", colors[i], markers1[i])
					if k > 1 {
						args = io.Sf("'k%s', ms=10", markers2[i])
					}
					if k == 0 && j == 0 {
						args += io.Sf(", label='P=%g'", Pval)
					}
					plt.PlotOne(R[i][j][k], St[i][j][k], args)
				}
			}
		}
		plt.Gll("$r\\;\\mathrm{[mm]}$", "$\\sigma_t\\;\\mathrm{[MPa]}$", "leg_loc='upper left'")
		plt.SaveD("/tmp", io.Sf("gofem_%s_st.eps", fnkey))
	}
}

func Test_flt02(tst *testing.T) {

	//verbose()
	chk.PrintTitle("flt02. circle with equality constraint")

	// parameters
	C := NewConfParams()
	C.Eps1 = 1e-3
	C.Pll = false
	C.Nisl = 4
	C.Ninds = 12
	C.Ntrials = 1
	if chk.Verbose {
		C.Ntrials = 40
	}
	C.Verbose = false
	C.Dtmig = 50
	C.CrowdSize = 3
	C.CompProb = false
	C.GAtype = "crowd"
	C.DiffEvol = true
	C.RangeFlt = [][]float64{
		{-1, 3}, // gene # 0: min and max
		{-1, 3}, // gene # 1: min and max
	}
	C.Latin = true
	C.PopFltGen = PopFltGen
	if chk.Verbose {
		C.FnKey = ""
		if C.Ntrials == 1 {
			C.DoPlot = true
		}
	}
	C.Ops.EnfRange = true
	C.NumFmts = map[string][]string{"flt": {"%8.4f", "%8.4f"}}
	C.ShowDem = true
	C.RegTol = 0.01
	C.CalcDerived()
	rnd.Init(C.Seed)

	// geometry
	xe := 1.0                      // centre of circle
	le := -0.4                     // selected level of f(x)
	ys := xe - (1.0+le)/math.Sqrt2 // coordinates of minimum point with level=le
	y0 := 2.0*ys + xe              // vertical axis intersect of straight line defined by c(x)
	xc := []float64{xe, xe}        // centre
	nx := len(xc)

	// functions
	fcn := func(f, g, h []float64, x []float64) {
		res := 0.0
		for i := 0; i < nx; i++ {
			res += (x[i] - xc[i]) * (x[i] - xc[i])
		}
		f[0] = math.Sqrt(res) - 1
		h[0] = x[0] + x[1] + xe - y0
	}

	// simple problem
	sim := NewSimpleFltProb(fcn, 1, 0, 1, C)
	sim.Run(chk.Verbose)

	// stat
	io.Pf("\n")
	sim.Stat(0, 60, -0.4)

	// plot
	sim.PltExtra = func() {
		plt.PlotOne(ys, ys, "'o', markeredgecolor='yellow', markerfacecolor='none', markersize=10")
	}
	sim.Plot("test_flt02")
}

//.........这里部分代码省略.........
	}
	if plot {
		if o.PlotNp < 3 {
			o.PlotNp = 41
		}
		var umin, umax, vmin, vmax float64
		if o.PlotCf < 1 {
			o.PlotCf = 2
		}
		if len(o.PlotUrange) == 0 {
			umin, umax = μ[0]-o.PlotCf*μ[0], μ[0]+o.PlotCf*μ[0]
			vmin, vmax = μ[1]-o.PlotCf*μ[1], μ[1]+o.PlotCf*μ[1]
		} else {
			chk.IntAssert(len(o.PlotUrange), 2)
			chk.IntAssert(len(o.PlotVrange), 2)
			umin, umax = o.PlotUrange[0], o.PlotUrange[1]
			vmin, vmax = o.PlotVrange[0], o.PlotVrange[1]
		}
		o.PlotU, o.PlotV = utl.MeshGrid2D(umin, umax, vmin, vmax, o.PlotNp, o.PlotNp)
		o.PlotZ = la.MatAlloc(o.PlotNp, o.PlotNp)
		plt.SetForEps(0.8, 300)
		for i := 0; i < o.PlotNp; i++ {
			for j := 0; j < o.PlotNp; j++ {
				o.xtmp[0] = o.PlotU[i][j]
				o.xtmp[1] = o.PlotV[i][j]
				o.PlotZ[i][j], err = o.gfcn(o.xtmp, args)
				if err != nil {
					chk.Panic("cannot compute gfcn(%v):\n%v", x, err)
				}
			}
		}
		plt.Contour(o.PlotU, o.PlotV, o.PlotZ, "")
		plt.ContourSimple(o.PlotU, o.PlotV, o.PlotZ, true, 8, "levels=[0], colors=['yellow']")
		plt.PlotOne(x[0], x[1], "'ro', label='initial'")
	}

	// iterations to find β
	var dat VarData
	B := []float64{β}
	itB := 0
	for itB = 0; itB < o.NmaxItB; itB++ {

		// message
		if verbose {
			gx, err := o.gfcn(x, args)
			if err != nil {
				chk.Panic("cannot compute gfcn(%v):\n%v", x, err)
			}
			io.Pf("%s itB=%d β=%g g=%g\n", io.StrThinLine(60), itB, β, gx)
		}

		// plot
		if plot {
			plt.PlotOne(x[0], x[1], "'r.'")
		}

		// compute direction cosines
		itA := 0
		for itA = 0; itA < o.NmaxItA; itA++ {

			// has lognormal random variable (lrv)
			if haslrv {

				// find equivalent normal mean and std deviation for lognormal variables
				for i := 0; i < nx; i++ {
					if o.lrv[i] {

// Plot plots contour
func (o *SimpleFltProb) Plot(fnkey string) {

	// check
	if !o.C.DoPlot {
		return
	}

	// limits and meshgrid
	xmin, xmax := o.C.RangeFlt[0][0], o.C.RangeFlt[0][1]
	ymin, ymax := o.C.RangeFlt[1][0], o.C.RangeFlt[1][1]

	// auxiliary variables
	X, Y := utl.MeshGrid2D(xmin, xmax, ymin, ymax, o.PltNpts, o.PltNpts)
	Zf := utl.DblsAlloc(o.PltNpts, o.PltNpts)
	var Zg [][][]float64
	var Zh [][][]float64
	if o.ng > 0 {
		Zg = utl.Deep3alloc(o.ng, o.PltNpts, o.PltNpts)
	}
	if o.nh > 0 {
		Zh = utl.Deep3alloc(o.nh, o.PltNpts, o.PltNpts)
	}

	// compute values
	x := make([]float64, 2)
	for i := 0; i < o.PltNpts; i++ {
		for j := 0; j < o.PltNpts; j++ {
			x[0], x[1] = X[i][j], Y[i][j]
			o.Fcn(o.ff[0], o.gg[0], o.hh[0], x)
			Zf[i][j] = o.ff[0][o.PltIdxF]
			for k, g := range o.gg[0] {
				Zg[k][i][j] = g
			}
			for k, h := range o.hh[0] {
				Zh[k][i][j] = h
			}
		}
	}

	// prepare plot area
	plt.Reset()
	plt.SetForEps(0.8, 350)

	// plot f
	if o.PltArgs != "" {
		o.PltArgs = "," + o.PltArgs
	}
	if o.PltCsimple {
		plt.ContourSimple(X, Y, Zf, true, 7, "colors=['k'], fsz=7"+o.PltArgs)
	} else {
		plt.Contour(X, Y, Zf, io.Sf("fsz=7, cmapidx=%d"+o.PltArgs, o.PltCmapIdx))
	}

	// plot g
	clr := "yellow"
	if o.PltCsimple {
		clr = "blue"
	}
	for _, g := range Zg {
		plt.ContourSimple(X, Y, g, false, 7, io.Sf("zorder=5, levels=[0], colors=['%s'], linewidths=[%g], clip_on=0", clr, o.PltLwg))
	}

	// plot h
	clr = "yellow"
	if o.PltCsimple {
		clr = "blue"
	}
	for _, h := range Zh {
		plt.ContourSimple(X, Y, h, false, 7, io.Sf("zorder=5, levels=[0], colors=['%s'], linewidths=[%g], clip_on=0", clr, o.PltLwh))
	}

	// initial populations
	l := "initial population"
	for _, pop := range o.PopsIni {
		for _, ind := range pop {
			x := ind.GetFloats()
			plt.PlotOne(x[0], x[1], io.Sf("'k.', zorder=20, clip_on=0, label='%s'", l))
			l = ""
		}
	}

	// final populations
	l = "final population"
	for _, pop := range o.PopsBest {
		for _, ind := range pop {
			x := ind.GetFloats()
			plt.PlotOne(x[0], x[1], io.Sf("'ko', ms=6, zorder=30, clip_on=0, label='%s', markerfacecolor='none'", l))
			l = ""
		}
	}

	// extra
	if o.PltExtra != nil {
		o.PltExtra()
	}

	// best result
	if o.Nfeasible > 0 {
		x, _, _, _ := o.find_best()
//.........这里部分代码省略.........

// Check checks derivatives
func Check(tst *testing.T, mdl Model, pc0, sl0, pcf float64, npts int, tolCc, tolD1a, tolD1b, tolD2a, tolD2b float64, verbose bool, pcSkip []float64, tolSkip float64, doplot bool) {

	// nonrate model
	nr_mdl, is_nonrate := mdl.(Nonrate)
	io.Pforan("is_nonrate = %v\n", is_nonrate)

	// for all pc stations
	Pc := utl.LinSpace(pc0, pcf, npts)
	Sl := make([]float64, npts)
	Sl[0] = sl0
	var err error
	for i := 1; i < npts; i++ {

		// update and plot
		Sl[i], err = Update(mdl, Pc[i-1], Sl[i-1], Pc[i]-Pc[i-1])
		if err != nil {
			tst.Errorf("Update failed: %v\n", err)
			return
		}
		if doplot {
			plt.PlotOne(Pc[i], Sl[i], "'ko', clip_on=0")
		}

		// skip point on checking of derivatives
		if doskip(Pc[i], pcSkip, tolSkip) {
			continue
		}

		// wetting flag
		wet := Pc[i]-Pc[i-1] < 0

		// check Cc = dsl/dpc
		io.Pforan("\npc=%g, sl=%g, wetting=%v\n", Pc[i], Sl[i], wet)
		if is_nonrate {

			// analytical Cc
			Cc_ana, err := mdl.Cc(Pc[i], Sl[i], wet)
			if err != nil {
				tst.Errorf("Cc failed: %v\n", err)
				return
			}

			// numerical Cc
			Cc_num, _ := num.DerivCentral(func(x float64, args ...interface{}) float64 {
				return nr_mdl.Sl(x)
			}, Pc[i], 1e-3)
			chk.AnaNum(tst, "Cc = ∂sl/∂pc    ", tolCc, Cc_ana, Cc_num, verbose)
		}

		// compute all derivatives
		L, Lx, J, Jx, Jy, err := mdl.Derivs(Pc[i], Sl[i], wet)
		if err != nil {
			tst.Errorf("Derivs failed: %v\n", err)
			return
		}
		L_ana_A := L
		L_ana_B, err := mdl.L(Pc[i], Sl[i], wet)
		if err != nil {
			tst.Errorf("L failed: %v\n", err)
			return
		}
		Lx_ana := Lx
		Jx_ana := Jx
		Jy_ana := Jy
		J_ana_A := J
		J_ana_B, err := mdl.J(Pc[i], Sl[i], wet)
		if err != nil {
			tst.Errorf("J failed: %v\n", err)
			return
		}

		// numerical L = ∂Cc/∂pc
		L_num, _ := num.DerivCentral(func(x float64, args ...interface{}) float64 {
			Cctmp, _ := mdl.Cc(x, Sl[i], wet)
			return Cctmp
		}, Pc[i], 1e-3)
		chk.AnaNum(tst, "L  = ∂Cc/∂pc    ", tolD1a, L_ana_A, L_num, verbose)

		// numerical Lx := ∂²Cc/∂pc²
		Lx_num, _ := num.DerivCentral(func(x float64, args ...interface{}) float64 {
			Ltmp, _, _, _, _, _ := mdl.Derivs(x, Sl[i], wet)
			return Ltmp
		}, Pc[i], 1e-3)
		chk.AnaNum(tst, "Lx = ∂²Cc/∂pc²  ", tolD2a, Lx_ana, Lx_num, verbose)

		// numerical J := ∂Cc/∂sl (version A)
		J_num, _ := num.DerivCentral(func(x float64, args ...interface{}) float64 {
			Ccval, _ := mdl.Cc(Pc[i], x, wet)
			return Ccval
		}, Sl[i], 1e-3)
		chk.AnaNum(tst, "J  = ∂Cc/∂sl    ", tolD1b, J_ana_A, J_num, verbose)

		// numerical Jx := ∂²Cc/(∂pc ∂sl)
		Jx_num, _ := num.DerivCentral(func(x float64, args ...interface{}) float64 {
			Ltmp, _, _, _, _, _ := mdl.Derivs(Pc[i], x, wet)
			return Ltmp
		}, Sl[i], 1e-3)
		chk.AnaNum(tst, "Jx = ∂²Cc/∂pc∂sl", tolD2b, Jx_ana, Jx_num, verbose)

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

//.........这里部分代码省略.........
		f0b := math.Log(a0/a1) / (b0 - b1)
		f1b := a0 * math.Exp(-b0*f0b)
		opt.Multi_fcnErr = func(f []float64) float64 {
			if f[0] < f0a {
				return f[1] - math.Exp(-f[0])
			}
			if f[0] < f0b {
				return f[1] - a0*math.Exp(-b0*f[0])
			}
			return f[1] - a1*math.Exp(-b1*f[0])
		}
		extraplot = func() {
			np := 201
			X, Y := utl.MeshGrid2D(0, 1, 0, 1, np, np)
			Z := utl.DblsAlloc(np, np)
			for j := 0; j < np; j++ {
				for i := 0; i < np; i++ {
					Z[i][j] = opt.Multi_fcnErr([]float64{X[i][j], Y[i][j]})
				}
			}
			plt.Contour(X, Y, Z, "levels=[0,0.6],cbar=0,lwd=0.5,fsz=5,cmapidx=6")
			F0 := utl.LinSpace(0, 1, 21)
			F1r := make([]float64, len(F0))
			F1s := make([]float64, len(F0))
			F1t := make([]float64, len(F0))
			for i, f0 := range F0 {
				F1r[i] = math.Exp(-f0)
				F1s[i] = a0 * math.Exp(-b0*f0)
				F1t[i] = a1 * math.Exp(-b1*f0)
			}
			plt.Plot(F0, F1r, "'k--',color='blue'")
			plt.Plot(F0, F1s, "'k--',color='green'")
			plt.Plot(F0, F1t, "'k--',color='gray'")
			plt.PlotOne(f0a, f1a, "'k|', ms=20")
			plt.PlotOne(f0b, f1b, "'k|', ms=20")
		}

	// problem # 2 -- CTP2, Deb 2001, p368/369, fig 226
	case 2:
		f1max = 1.2
		θ, a, b := -0.2*PI, 0.2, 10.0
		c, d, e := 1.0, 6.0, 1.0
		fcn = CTPgenerator(θ, a, b, c, d, e)
		extraplot = CTPplotter(θ, a, b, c, d, e, f1max)
		opt.Multi_fcnErr = CTPerror1(θ, a, b, c, d, e)

	// problem # 3 -- CTP3, Deb 2001, p368/370, fig 227
	case 3:
		f1max = 1.2
		θ, a, b := -0.2*PI, 0.1, 10.0
		c, d, e := 1.0, 0.5, 1.0
		fcn = CTPgenerator(θ, a, b, c, d, e)
		extraplot = CTPplotter(θ, a, b, c, d, e, f1max)
		opt.Multi_fcnErr = CTPerror1(θ, a, b, c, d, e)

	// problem # 4 -- CTP4, Deb 2001, p368/370, fig 228
	case 4:
		f1max = 2.0
		θ, a, b := -0.2*PI, 0.75, 10.0
		c, d, e := 1.0, 0.5, 1.0
		fcn = CTPgenerator(θ, a, b, c, d, e)
		extraplot = CTPplotter(θ, a, b, c, d, e, f1max)
		opt.Multi_fcnErr = CTPerror1(θ, a, b, c, d, e)

	// problem # 5 -- CTP5, Deb 2001, p368/371, fig 229
	case 5: