Skip to content
This repository has been archived by the owner on Jan 2, 2018. It is now read-only.

Commit

Permalink
*: add extra output flags
Browse files Browse the repository at this point in the history 
Signed-off-by: Aleksa Sarai <[email protected]>
  • Loading branch information
@cyphar
cyphar committed Oct 23, 2016
1 parent a0372f8 commit e836cc1
Show file tree
Hide file tree
Showing 3 changed files with 154 additions and 137 deletions.
22 changes: 0 additions & 22 deletions src/searcher.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,33 +22,11 @@
import argparse
import random

import matplotlib
matplotlib.use("TkAgg")
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid.anchored_artists import AnchoredText
import numpy
import numpy.random
import scipy
import scipy.integrate

SPINE_COLOR = "black"

def latexify(ax):
for spine in ["top", "right"]:
ax.spines[spine].set_visible(False)

for spine in ["left", "bottom"]:
ax.spines[spine].set_color(SPINE_COLOR)
ax.spines[spine].set_linewidth(0.5)

ax.xaxis.set_ticks_position("bottom")
ax.yaxis.set_ticks_position("left")

for axis in [ax.xaxis, ax.yaxis]:
axis.set_tick_params(direction="out", color=SPINE_COLOR)

return ax

def csv_column_write(f, cols, fieldnames):
import csv

Expand Down
73 changes: 53 additions & 20 deletions src/shoebox.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,19 @@
import matplotlib.pyplot as plt
import cmocean

matplotlib.rcParams.update({
"backend": "ps",
"text.usetex": True,
"text.latex.preamble": [r"\usepackage[dvips]{graphicx}"],
"axes.labelsize": 13, # fontsize for x and y labels (was 10)
"axes.titlesize": 13,
"font.size": 13, # was 10
"legend.fontsize": 13, # was 10
"xtick.labelsize": 13,
"ytick.labelsize": 13,
"font.family": "serif", # ???
})

def positions(ndarray):
return zip(*numpy.where(numpy.ones_like(ndarray)))

Expand Down Expand Up @@ -131,44 +144,64 @@ def main(config):
theta_space = config.theta_space
config.eta_space = config.eta_space[0]

f, plots = plt.subplots(2, size, sharex='col', sharey='row')
f, plots = plt.subplots(1 + config.best_phi, size, figsize=(5, 5), dpi=80, sharex='col', sharey='row')

if numpy.array([plots]).shape == (1,):
plots = numpy.array([plots])

if len(plots.shape) == 1:
plots = numpy.array([plots]).T

for idx, (ax1, ax2) in enumerate(plots.T):
ax1 = latexify(ax1)
ax2 = latexify(ax2)
for idx, axes in enumerate(plots.T):
ax1 = ax2 = None
if len(axes) >= 1:
ax1 = latexify(axes[0])
if len(axes) >= 2:
ax2 = latexify(axes[1])

config.theta_space = theta_space[idx]
print(config.theta_space)
pcm1 = plot_shoebox(config, ax1, config.file[0], metric=config.metric)
pcm2 = plot_shoebox(config, ax2, config.file[0], metric=config.metric+"_phi")
# pcm1 = plot_shoebox(ax1, config.file[0], metric="linear")
# pcm2 = plot_shoebox(ax2, config.file[0], metric="depth")

f.colorbar(pcm1, ax=ax1)
if config.angle_ticks:
ticks = numpy.linspace(0, 2, num=9)
cbar = f.colorbar(pcm2, ax=ax2, ticks=[x*math.pi for x in ticks])
cbar.ax.set_yticklabels([r"$%s \pi$" % (x,) for x in ticks])
else:
cbar = f.colorbar(pcm2, ax=ax2)

ax1.set_ylabel(r"$\eta$")
ax2.set_xlabel(r"$\theta$")
if ax1 is not None:
pcm1 = plot_shoebox(config, ax1, config.file[0], metric=config.metric)
if ax2 is not None:
pcm2 = plot_shoebox(config, ax2, config.file[0], metric=config.metric+"_phi")

if ax1 is not None:
f.colorbar(pcm1, ax=ax1)
if ax2 is not None:
if config.angle_ticks:
ticks = numpy.linspace(0, 2, num=9)
cbar = f.colorbar(pcm2, ax=ax2, ticks=[x*math.pi for x in ticks])
cbar.ax.set_yticklabels([r"$%s \pi$" % (x,) for x in ticks])
else:
cbar = f.colorbar(pcm2, ax=ax2)

if ax1 is not None:
ax1.set_ylabel(r"$\eta$")
if ax1 is None:
ax1.set_xlabel(r"$\theta$")
if ax2 is not None:
ax2.set_xlabel(r"$\theta$")

# plt.legend()
f.subplots_adjust(hspace=-0.2)
f.tight_layout()
plt.show()

if config.out:
plt.savefig(config.out)
else:
plt.show()

if __name__ == "__main__":
def __wrapped_main__():
parser = argparse.ArgumentParser(description="Plots the shape of the parameter space (eta, theta, phi) with the colour being set by the given metric value.")
# save argument
parser.add_argument("--out", dest="out", type=str, default=None, help="Output to given filename rather than displaying in an interactive window (default: disabled)")
# metric arguments
parser.add_argument("-m", "--metric", dest="metric", type=str, default="depth", help="Metric to plot from {depth, linear} (default: [depth]).")
# plotting arguments
parser.add_argument("-sp", "--show-best-phi", dest="best_phi", action="store_const", const=True, default=False, help="Show best phi plot.")
parser.add_argument("--no-show-best-phi", dest="best_phi", action="store_const", const=False, default=False, help="Do not show best phi plot. (default)")
parser.add_argument("-a", "--angle-ticks", dest="angle_ticks", action="store_const", const=True, default=True, help="Use angle ticks (have multiples of pi when plotting angles) (default).")
parser.add_argument("--no-angle-ticks", dest="angle_ticks", action="store_const", const=False, default=True, help="Do not use angle ticks.")
# shoebox arguments
Expand Down
196 changes: 101 additions & 95 deletions src/single.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,19 @@
import scipy
import scipy.integrate

matplotlib.rcParams.update({
"backend": "ps",
"text.usetex": True,
"text.latex.preamble": [r"\usepackage[dvips]{graphicx}"],
"axes.labelsize": 13, # fontsize for x and y labels (was 10)
"axes.titlesize": 13,
"font.size": 13, # was 10
"legend.fontsize": 13, # was 10
"xtick.labelsize": 13,
"ytick.labelsize": 13,
"font.family": "serif", # ???
})

SPINE_COLOR = "black"

def latexify(ax):
Expand Down Expand Up @@ -169,83 +182,17 @@ def solve(t0, t1, dt, start):
ts, As = numpy.array(As).T
return ts, As

def metric_depth(ts, As, eta, theta, phi):
mid = middle(eta, theta, phi)
Is = numpy.vectorize(abs)(As)

# Cut down time domain so that we don't detect minima before the estimate.
filt = ts >= mid
ts = ts[filt]
Is = Is[filt]

# Get lowest dip value.
i_dip = numpy.argmin(Is)
t_dip = ts[i_dip]
I_dip = Is[i_dip]

# Get highest peak value that happened after mid but before the dip.
filt = ts < t_dip
# ts = ts[filt]
Is = Is[filt]

if not Is.shape[0]:
return None

# TODO: Handle cases where the maximum is before the peak.
I_peak = numpy.amax(Is)

return (I_peak - I_dip) / I_peak

def metric_linear(ts, As, eta, theta, phi):
mid = middle(eta, theta, phi)
Is = numpy.vectorize(abs)(As)

# First-order linear approximation.
n1 = GAMMA * numpy.cos(theta)
n2 = GAMMA * numpy.sin(theta) * numpy.exp(1j * phi)
linAs = eta * numpy.exp(GAMMA * ts) * (n1 * numpy.exp(1j * GAMMA * ts) + n2 * numpy.exp(-1j * GAMMA * ts))
linIs = numpy.vectorize(abs)(linAs)

# Cut down the time domain so we don't take into account anything before the estimated peak.
filt = ts >= mid
ts = ts[filt]
Is = Is[filt]
linIs = linIs[filt]

# Remove entries once the non-linear terms become out-of-hand.
grad = numpy.gradient(Is)
filt = grad < 100
ts = ts[filt]
Is = Is[filt]
linIs = linIs[filt]

# We'll operate on the log values for sanity reasons.
Is = numpy.log(Is)
linIs = numpy.log(linIs)

# Get the difference between linIs and Is and compute a metric from that.
diff = linIs - Is

# We need to ignore parts when Is > linIs.
filt = Is <= linIs
ts = ts[filt]
diff = diff[filt]
return abs(numpy.trapz(diff, ts))
# return numpy.sum(diff)

METRICS = {
"depth": metric_depth,
"linear": metric_linear,
}

def main(config):
eta = config.eta
theta = config.theta
phi = config.phi

fig = plt.figure(figsize=(10, 10), dpi=80)
ax1 = latexify(fig.add_subplot("211"))
ax2 = latexify(fig.add_subplot("212"))
fig = plt.figure(figsize=(5, 5), dpi=80)
ax1 = latexify(fig.add_subplot("%d11" % (1+config.show_phase,)))
if config.show_phase:
ax2 = latexify(fig.add_subplot("212"))
else:
ax2 = None

t0 = 0
mid = middle(eta, theta, phi)
Expand All @@ -255,53 +202,112 @@ def main(config):
ts, As = solve(t0, t1, config.dt, start)

# now plot theoretical
ax1.set_title(r"{$\theta = %s, \eta = %s, \phi = %s$}" % (theta, eta, phi))
# ax1.set_title(r"{$\theta = %s, \eta = %s, \phi = %s$}" % (theta, eta, phi))

# plot integration
ax1.plot(ts, numpy.vectorize(abs)(As), 'k')
ax2.plot(ts, numpy.vectorize(cmath.phase)(As), 'k')

ax1.set_yscale("log")
if ax2 is not None:
ax2.plot(ts, numpy.vectorize(cmath.phase)(As), 'k')

# plot first order approximation.
if config.first_order:
n1 = GAMMA * numpy.cos(theta)
n2 = GAMMA * numpy.sin(theta) * numpy.exp(1j * phi)
linAs = eta * numpy.exp(GAMMA * ts) * (n1 * numpy.exp(1j * GAMMA * ts) + n2 * numpy.exp(-1j * GAMMA * ts))
ax1.plot(ts, numpy.vectorize(abs)(linAs), 'r--')
if ax2 is not None:
ax2.plot(ts, numpy.vectorize(cmath.phase)(linAs), 'r--')

if config.fill_area:
Is = numpy.vectorize(abs)(As)
linIs = numpy.vectorize(abs)(linAs)

# Cut down the time domain so we don't take into account anything before the estimated peak.
filt = ts >= mid
ts = ts[filt]
Is = Is[filt]
linIs = linIs[filt]

# Remove entries once the non-linear terms become out-of-hand.
grad = numpy.gradient(Is)
filt = grad < 100
ts = ts[filt]
Is = Is[filt]
linIs = linIs[filt]

# We need to ignore parts when Is > linIs.
filt = Is <= linIs
ts = ts[filt]
Is = Is[filt]
linIs = linIs[filt]

# Fill.
ax1.fill_between(ts, Is, linIs, color='orange')

ax1.set_yscale(config.scale)
ax1.set_axisbelow(True)
ax1.set_xlabel(r"Time ($\tau$)")
ax1.set_ylabel(r"Amplitude ($|A|$)")
ax1.set_ylabel(r"Amplitude ($|b|$)")
ax1.set_xlim([t0, t1])

ax1.xaxis.set_ticks([mid], minor=True)
ax1.xaxis.grid(True, which="minor", color="k", linestyle=":")
ax2.xaxis.set_ticks([mid], minor=True)
ax2.xaxis.grid(True, which="minor", color="k", linestyle=":")

ax2.set_axisbelow(True)
ax2.set_xlabel(r"Time ($\tau$)")
ax2.set_ylabel(r"Phase")
ax2.set_xlim([t0, t1])
ax2.set_ylim([-math.pi, math.pi])

# Set ticks to represent angle.
ticks = numpy.linspace(-1, 1, num=9)
ax2.set_yticks([x*math.pi for x in ticks])
ax2.yaxis.set_major_formatter(matplotlib.ticker.FixedFormatter([r"$%s \pi$" % (x,) for x in ticks]))
if config.scale == "log":
ax1.set_ylim([None, 1e11])
elif config.scale == "linear":
ax1.set_ylim([0, 3])

# ax1.xaxis.set_ticks([mid], minor=True)
# ax1.xaxis.grid(True, which="minor", color="k", linestyle=":")
# if ax2 is not None:
# ax2.xaxis.set_ticks([mid], minor=True)
# ax2.xaxis.grid(True, which="minor", color="k", linestyle=":")

if ax2 is not None:
ax2.set_axisbelow(True)
ax2.set_xlabel(r"Time ($\tau$)")
ax2.set_ylabel(r"Phase")
ax2.set_xlim([t0, t1])
ax2.set_ylim([-math.pi, math.pi])

# Set ticks to represent angle.
ticks = numpy.linspace(-1, 1, num=9)
ax2.set_yticks([x*math.pi for x in ticks])
ax2.yaxis.set_major_formatter(matplotlib.ticker.FixedFormatter([r"$%s \pi$" % (x,) for x in ticks]))

plt.legend()
fig.tight_layout()

plt.show()
if config.out:
plt.savefig(config.out)
else:
plt.show()


if __name__ == "__main__":
def __wrapped_main__():
parser = argparse.ArgumentParser(description="Searches the parameter space (eta, theta, phi) for the PQS ODE using the set of metrics given.")
# save plot
parser.add_argument("--out", dest="out", type=str, default=None, help="Output to given filename rather than displaying in an interactive window (default: disabled)")
# integrator arguments
parser.add_argument("-dt", dest="dt", type=float, default=0.01, help="Time spacing for integration (default: 0.01).")
# metric arguments
# parser.add_argument("-m", "--metric", dest="metrics", action="append", default=["depth"], help="Set of metrics from {depth, linear} to compute (default: [depth]).")
# plot arguments
parser.add_argument("-fo", "--first-order", dest="first_order", default=False, action="store_const", const=True, help="Plot the first order approximation of the curve.")
parser.add_argument("--no-first-order", dest="first_order", default=False, action="store_const", const=False, help="Do not plot the first order approximation of the curve (default).")
parser.add_argument("-fa", "--fill-area", dest="fill_area", default=False, action="store_const", const=True, help="Fill the area between the first order approximation and the integral.")
parser.add_argument("--no-fill-area", dest="fill_area", default=False, action="store_const", const=False, help="Do not fill the area between the first order approximation and the integral (default).")
parser.add_argument("-pp", "--plot-phase", dest="show_phase", default=False, action="store_const", const=True, help="Show the phase of the curve.")
parser.add_argument("--no-plot-phase", dest="show_phase", default=False, action="store_const", const=False, help="Do not show the phase of the curve. (default)")
parser.add_argument("-sc", "--scale", dest="scale", default="log", help="The scaling to use for the y axis (default: log).")
# parameter arguments
parser.add_argument("-e", "--eta", dest="eta", type=float, required=True, help="Eta value.")
parser.add_argument("-t", "--theta", dest="theta", type=float, required=True, help="Theta value.")
parser.add_argument("-p", "--phi", dest="phi", type=float, required=True, help="Phi value.")

config = parser.parse_args()

if config.fill_area and not config.first_order:
print("--fill-area requires --first-order")
sys.exit(1)

main(config)

__wrapped_main__()

0 comments on commit e836cc1

Please sign in to comment.