|
#! /usr/bin/env python
|
|
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
|
|
N = 50
|
|
XMIN, XMAX = 0.0, 2.5
|
|
F = 0.6
|
|
|
|
def f(x):
|
|
"Function to plot"
|
|
#return np.exp(-x) * np.cos(10*x)
|
|
return (2.5 - x) + np.exp(-2*x) * np.cos(20*x**0.7)
|
|
|
|
## Naive version
|
|
xs_lin = np.linspace(XMIN, XMAX, N)
|
|
ys_lin = f(xs_lin)
|
|
|
|
## Densely sample to build an cdf of |f''|
|
|
xs = np.linspace(XMIN, XMAX, 1000)
|
|
ys = f(xs)
|
|
ypps = np.diff(ys, 2)
|
|
#ypps /= (ys[1:-1] + np.mean(ys))
|
|
ayps = np.cumsum(abs(ypps))
|
|
ayps /= ayps[-1]
|
|
|
|
## Mix with linear in F:1-F ratio
|
|
byps = np.cumsum([1. for x in xs[1:-1]])
|
|
byps /= byps[-1]
|
|
ayps = F*ayps + (1-F)*byps
|
|
|
|
## Invert cdf
|
|
# TODO: make more efficient!
|
|
def G(q):
|
|
return np.interp(q, ayps, xs[1:-1])
|
|
qs = np.linspace(0, 1, N)
|
|
xs_opt = G(qs)
|
|
ys_opt = f(xs_opt)
|
|
|
|
## Plot
|
|
for x in xs_opt:
|
|
plt.axvline(x, color="lightgray", linewidth=0.8)
|
|
plt.plot(xs, ys, label="$f(x)$")
|
|
plt.plot(xs[1:-1], ayps, "--", label="CDF, $F(x)$")
|
|
#plt.plot(xs[1:-1], ayps, "*")
|
|
plt.plot(xs_lin, ys_lin, label="Linear")
|
|
plt.plot(xs_opt, ys_opt, label="Optimal")
|
|
plt.legend()
|
|
for f in [".png", ".pdf"]:
|
|
plt.savefig("curvplot2"+f, dpi=150)
|
|
#plt.show()
|