#!/usr/bin/env python
# coding: utf-8
# $Id: plot_interpolation_test.py
# Author: Daniel R. Reese <daniel.reese@obspm.fr>
# Copyright (C) Daniel R. Reese and contributors
# Copyright license: GNU GPL v3.0
#
# AIMS is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with AIMS. If not, see <http://www.gnu.org/licenses/>.
#
"""
An interactive utility which plots various forms of interpolation error,
stored in a binary file produced by :py:func:`AIMS.test_interpolation`.
It specifically tests the errors from two types of interpolation:
- age interpolation: this is interpolation along a given evolutionary track
- track interpolation: this is interpolation between different evolutionary
tracks
This utility allows various types of plots:
- 3D plots of interpolation errors as a function of grid structural parameters
- 2D slices which show interpolation errors as a function of age for a
given evolutionary track
- interactive plots which allow you to select 2D slices
.. note::
Interpolation errors for models in a given evolutionary track are typically
stored in arrays as follows:
- ``result[model_number,ndim+0]`` = maximum error on the radial modes
- ``result[model_number,ndim+1]`` = RMS error on the radial modes
- ``result[model_number,ndim+2]`` = RMS error on the radial modes near
:math:`\\nu_{\mathrm{max}}`
- ``result[model_number,ndim+3]`` = maximum error on the non radial modes
- ``result[model_number,ndim+4]`` = RMS error on the non radial modes
- ``result[model_number,ndim+5]`` = RMS error on the non radial modes near
:math:`\\nu_{\mathrm{max}}`
- ``result[model_number,ndim+6+[0:nglb]]`` = errors on the global parameters
where:
- ``result`` = the array which containts the interpolation errors
- ``model_numer`` = an index which represents the model (not necessarily the
number of the model along the evolutionary track
- ``ndim`` = the number of dimensions in the grid (including age)
- ``nglb`` = the number of global parameters for stellar models in the grid
.. warning::
This plot utility only works with 3 dimensional grids (incl. the age dimension).
"""
__docformat__ = 'restructuredtext'
import sys
import dill
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
# global variables
ndim = 0
""" number of dimension in grid (including age) """
nglb = 0
""" number of global parameters """
titles = None
""" the grid quantities, which will serve as axis labels """
results_age = None
""" list which contains the arrays with the errors from age interpolation """
results_track = None
""" list which contains the arrays with the errors from track interpolation """
[docs]def all_nan(array):
"""
Test to see if all of the elements of an array are nan's.
:param array: array in which we're checking to see if all elements are nan's.
:type array: np.array
:return: ``True`` if all the elements of ``array`` are nan's, and ``False``
otherwise.
:rtype: boolean
"""
return np.all(np.isnan(array))
[docs]def plot3D(results,error_ndx,tpe="max",title=None,truncate=0):
"""
Create 3D plot showing the error as a function of the two
first grid parameters.
:param results: list of 2D arrays which contain various types of errors as
a function of the model number along a given evolutionary track.
:param error_ndx: value which specifies the type of error to be plotted.
:param tpe: specifies how to combine errors along the
evolutionary track. Options include:
- "max": corresponds to taking the maximum value.
- "avg": takes the root mean-square value.
:param title: the title of the plot
:param truncate: (default = 0): specifies how many models should be omitted
on both ends of the track. This is useful for comparing results from
tests involing different sizes of increments.
:type results: list of np.arrays
:type error_ndx: int
:type tpe: string
:type title: string
:type truncate: int
.. note::
See above introductory description for a more detailled description
of the indices which intervene in the 2D arrays contained in ``results``
and of the relevant values for ``error_ndx``.
"""
n = len(results)
x = []
y = []
z = []
for i in xrange(n):
start = truncate
stop = results[i].shape[0] - truncate
if (tpe == "max"):
value = np.nanmax(results[i][start:stop,ndim+error_ndx])
elif (tpe == "avg"):
num = den = 0.0
for j in xrange(start,stop):
if (not math.isnan(results[i][j,ndim+error_ndx])):
num += results[i][j,ndim+error_ndx]**2
den += 1.0
if (den > 0.0):
value = math.sqrt(num/den)
else:
value = np.nan
else:
print "ERROR: unrecognised type: ",tpe
sys.exit(1)
if (value > 0.0):
z.append(math.log10(value))
x.append(results[i][0,0])
y.append(results[i][0,1])
x = np.array(x,dtype = np.float64)
y = np.array(y,dtype = np.float64)
z = np.array(z,dtype = np.float64)
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(x,y,z,cmap=cm.jet,linewidth=0.2)
ax.set_xlabel(titles[0])
ax.set_ylabel(titles[1])
ax.set_zlabel(r"$\log_{10}$(%s. error)"%(tpe))
if (title is not None): ax.set_title(title)
[docs]def plot_grid(grid):
"""
Make an interactive plot of the grid. Clicking on the blue dots
will produce a 2D slice showing age interpolation errors for
the associated evolutionary track.
:param grid: array containing basic grid parameters (excluding age)
:type grid: np.array
.. warning::
This only works for two-dimensional grids (excluding the age
dimension).
"""
# remember ndim includes the age dimension.
if (ndim != 3):
print "Only able to plot the tessellation in two dimensions."
return
# find bounds:
xmin = np.nanmin(grid[:,0])
xmax = np.nanmax(grid[:,0])
ymin = np.nanmin(grid[:,1])
ymax = np.nanmax(grid[:,1])
dx = xmax-xmin
dy = ymax/ymin
xmin -= dx*0.03
xmax += dx*0.03
ymin /= dy**0.05
ymax *= dy**0.05
fig = plt.figure()
#plt.semilogy(grid[:,0],grid[:,1],'bo', picker=5)
plt.plot(grid[:,0],grid[:,1],'bo', picker=5)
plt.xlim((xmin,xmax))
plt.ylim((ymin,ymax))
plt.xlabel(titles[0])
plt.ylabel(titles[1])
fig.canvas.mpl_connect('pick_event', onpick_age)
[docs]def onpick_age(event):
"""
Event catcher for the grid plot (which shows the positions of
the evolutionary tracks as a function of the grid parameters,
excluding age).
Parameters:
:param event: event caught by the grid plot.
"""
for pos in event.ind:
plot_slice_age(pos)
[docs]def plot_slice_age(pos):
"""
Plot age interpolation error as a function of age for a
given track.
:param pos: index of the relevant track.
:type pos: int
.. note::
This `pos` index applies to results_age, i.e., it is based
on the original track indices.
"""
if (all_nan(results_age[0][pos][:,ndim:ndim+6]) and \
all_nan(results_age[1][pos][:,ndim:ndim+6])):
print "Cowardly refusing to plot nan's"
return
style = ["b:","b-","bv","r:","r-","rv"]
labels = ["Max. (incr=1)", "Avg. (incr=1)", r"$\nu_{\mathrm{max}}$ avg (incr=1)",
"Max. (incr=2)", "Avg. (incr=2)", r"$\nu_{\mathrm{max}}$ avg (incr=2)"]
xmin = np.nanmin(results_age[0][pos][:,ndim-1])
xmax = np.nanmax(results_age[1][pos][:,ndim-1])
xlim = (xmin,xmax+0.5*(xmax-xmin))
plt.figure()
plt.subplot(2,1,1)
if (all_nan(results_age[0][pos][:,ndim:ndim+3]) and \
all_nan(results_age[1][pos][:,ndim:ndim+3])):
plt.plot([0,1,1,0],[0,1,0,1],"k-")
else:
for j in xrange(2):
for i in xrange(3):
plt.plot(results_age[j][pos][:,ndim-1],results_age[j][pos][:,ndim+i],
style[i+3*j],label=labels[i+3*j])
plt.title(r"Error at $(M,Z)=(%f,%f)$"%(results_age[0][pos][0,0],results_age[0][pos][0,1]))
plt.ylabel(r"Error(radial)")
plt.yscale('log')
plt.xlim(xlim)
plt.legend()
plt.subplot(2,1,2)
if (all_nan(results_age[0][pos][:,ndim+3:ndim+6]) and \
all_nan(results_age[1][pos][:,ndim+3:ndim+6])):
plt.plot([0,1,1,0],[0,1,0,1],"k-")
else:
for j in xrange(2):
for i in xrange(3):
plt.plot(results_age[j][pos][:,ndim-1],results_age[j][pos][:,ndim+3+i],
style[i+3*j],label=labels[i+3*j])
plt.ylabel(r"Error(non radial)")
plt.yscale('log')
plt.xlabel(titles[ndim-1])
plt.xlim(xlim)
plt.legend()
plt.show()
[docs]def plot_partition_tessellation(grid, ndx1, ndx2, tessellation):
"""
Make an interactive tessellation plot based on the supplied partition
on the grid. Clicking on the blue dots will produce a 2D slice showing
track interpolation errors for the associated evolutionary track.
:param grid: array containing basic grid parameters (excluding age)
:param ndx1: list with the indices of the first part of the partition.
:param ndx2: list with the indices of the second part of the partition.
:param tessellation: grid tessellation associated with ``ndx2``
:type grid: np.array
:type ndx1: list of int
:type ndx2: list of int
.. warning::
This only works for two-dimensional tessellations.
"""
# remember ndim includes the age dimension which is not included
# in the tessellation:
if (ndim != 3):
print "Only able to plot the tessellation in two dimensions."
return
# find bounds:
xmin = np.nanmin(grid[:,0])
xmax = np.nanmax(grid[:,0])
ymin = np.nanmin(grid[:,1])
ymax = np.nanmax(grid[:,1])
dx = xmax-xmin
dy = ymax/ymin
xmin -= dx*0.03
xmax += dx*0.03
ymin /= dy**0.05
ymax *= dy**0.05
fig = plt.figure()
#plt.semilogy(grid[ndx1,0],grid[ndx1,1],'bo', picker=5)
#plt.semilogy(grid[ndx2,0],grid[ndx2,1],'yo')
plt.plot(grid[ndx1,0],grid[ndx1,1],'bo', picker=5)
plt.plot(grid[ndx2,0],grid[ndx2,1],'yo')
plt.triplot(grid[ndx2,0],grid[ndx2,1],tessellation.simplices.copy())
plt.xlim((xmin,xmax))
plt.ylim((ymin,ymax))
plt.xlabel(titles[0])
plt.ylabel(titles[1])
fig.canvas.mpl_connect('pick_event', onpick_track)
[docs]def onpick_track(event):
"""
Event catcher for the partition tessellation plot (associated with
tests of track interpolation).
:param event: event caught by the partition tessellation plot.
"""
for pos in event.ind:
plot_slice_track(pos)
[docs]def plot_slice_track(pos):
"""
Plot track interpolation error as a function of age for a
given track.
:param pos: index of the relevant track.
:type pos: int
.. note::
This `pos` index applies to results_track, i.e., it is based
on the indices deduced from the grid partition.
"""
if (all_nan(results_track[pos][:,ndim:ndim+6])):
print "Cowardly refusing to plot nan's"
return
style = ["b:","b-","bv"]
labels = ["Max.", "Avg.", r"$\nu_{\mathrm{max}}$ avg"]
xmin = np.nanmin(results_track[pos][:,ndim-1])
xmax = np.nanmax(results_track[pos][:,ndim-1])
xlim = (xmin,xmax+0.3*(xmax-xmin))
plt.figure()
plt.subplot(2,1,1)
if (all_nan(results_track[pos][:,ndim:ndim+3])):
plt.plot([0,1,1,0],[0,1,0,1],"k-")
else:
for i in xrange(3):
plt.plot(results_track[pos][:,ndim-1],results_track[pos][:,ndim+i],
style[i],label=labels[i])
plt.title(r"Error at $(M,Z)=(%f,%f)$"%(results_track[pos][0,0],results_track[pos][0,1]))
plt.ylabel(r"Error(radial)")
plt.yscale('log')
plt.xlim(xlim)
plt.legend()
plt.subplot(2,1,2)
if (all_nan(results_track[pos][:,ndim+3:ndim+6])):
plt.plot([0,1,1,0],[0,1,0,1],"k-")
else:
for i in xrange(3):
plt.plot(results_track[pos][:,ndim-1],results_track[pos][:,ndim+3+i],
style[i],label=labels[i])
plt.ylabel(r"Error(non radial)")
plt.yscale('log')
plt.xlabel(titles[ndim-1])
plt.xlim(xlim)
plt.legend()
plt.show()
if __name__ == "__main__":
"""
Plot interpolation tests.
.. note::
The user may want to edit this.
"""
if (len(sys.argv) < 2):
print "Usage: plot_interpolation_test.py data_file"
sys.exit(1)
filename = sys.argv[1]
input_data = open(filename,"r")
[ndim, nglb, titles, grid, ndx1, ndx2, tessellation, results_age1, \
results_age2, results_track] = dill.load(input_data)
if (ndim != 3):
print "I'm sorry, but I can only plot handle 3D grids."
sys.exit(1)
results_age = [results_age1, results_age2]
input_data.close()
plot3D(results_age1,0,tpe="max",title="Max radial error (nincr = 1)",truncate=1)
plot3D(results_age2,0,tpe="max",title="Max radial error (nincr = 2)")
plot3D(results_track,0,tpe="max",title="Max radial error (struct)")
plot3D(results_age1,1,tpe="avg",title="Avg radial error (nincr = 1)",truncate=1)
plot3D(results_age2,1,tpe="avg",title="Avg radial error (nincr = 2)")
plot3D(results_track,1,tpe="avg",title="Avg radial error (struct)")
plot3D(results_age1,2,tpe="avg",title="numax radial error (nincr = 1)",truncate=1)
plot3D(results_age2,2,tpe="avg",title="numax radial error (nincr = 2)")
plot3D(results_track,2,tpe="avg",title="numax radial error (struct)")
plot_grid(grid)
plot_partition_tessellation(grid, ndx1, ndx2, tessellation)
plt.show()