Golang plt.PlotOne函数代码示例

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