diff --git a/bin/old_viewing_software/qtView.py b/bin/old_viewing_software/qtView.py index 90e8b5a..d61387c 100644 --- a/bin/old_viewing_software/qtView.py +++ b/bin/old_viewing_software/qtView.py @@ -5,18 +5,20 @@ import argparse import random -import cPickle +import pickle import numpy as np -from PyQt4.QtCore import * -from PyQt4.QtGui import * +from PyQt5.QtCore import * +from PyQt5.QtGui import * +from PyQt5.QtWidgets import * import matplotlib from matplotlib.figure import Figure +from matplotlib.backends.backend_qt5agg import ( + FigureCanvasQTAgg as FigCanvas, + NavigationToolbar2QT as NavigationToolbar +) -from matplotlib.backends.backend_qt4agg import \ - FigureCanvasQTAgg as FigCanvas, \ - NavigationToolbar2QTAgg as NavigationToolbar class ViewException(Exception): pass @@ -66,7 +68,7 @@ def makeDataObj(frac=None, extnList=['.model', '.results']): 'Array Index'] if len(textlabels) != len(axislabels): - print "Label list lengths not identical." + print("Label list lengths not identical.") sys.exit() dataObjList = [] @@ -78,10 +80,10 @@ def makeDataObj(frac=None, extnList=['.model', '.results']): try: f = open(filename, 'rb') except IOError: - print "Could not open file {0}.".format(filename) + print("Could not open file {0}.".format(filename)) sys.exit() - pop = cPickle.load(f) + pop = pickle.load(f) f.close() # create numpy array of right size diff --git a/bin/old_viewing_software/tkView.py b/bin/old_viewing_software/tkView.py index 924292d..21fe2b2 100644 --- a/bin/old_viewing_software/tkView.py +++ b/bin/old_viewing_software/tkView.py @@ -4,17 +4,21 @@ import os import argparse import random +import warnings -import cPickle +import pickle import numpy as np +from scipy.stats import norm +from scipy.optimize import curve_fit +from scipy.signal import argrelmax -import Tkinter as tk +import tkinter as tk import matplotlib matplotlib.use('TKAgg') from matplotlib.figure import Figure from matplotlib.backends.backend_tkagg import \ - FigureCanvasTkAgg as FigCanvas, \ - NavigationToolbar2TkAgg as NavigationToolbar + FigureCanvasTkAgg, \ + NavigationToolbar2Tk as NavigationToolbar class ViewException(Exception): pass @@ -38,8 +42,9 @@ def makeDataObj(frac=None, extnList=['.model', '.results']): 'Gal Y', 'Gal Z', 'W', - 'alpha', - 'rho', + 'Pdot', + 'Age', + 'B', 'SI', 'S1400', 'gl', @@ -52,34 +57,36 @@ def makeDataObj(frac=None, extnList=['.model', '.results']): 'X (kpc)', 'Y (kpc)', 'Z (kpc)', - 'Width (degrees)', - 'alpha (deg)', - 'rho (deg)', + 'Width (°)', + 'Pdot', + 'Age [P/2Pdot] (yr)', + 'B [(PPdot)^(1/2)] (G)', 'Spectral Index', 'S1400 (mJy)', - 'Galactic Longitude (degrees)', - 'Galactic Latitude (degrees)', + 'Galactic Longitude (°)', + 'Galactic Latitude (°)', 'Distance (kpc)', 'GalacticRadius (kpc)', 'Array Index'] if len(textlabels) != len(axislabels): - print "Label list lengths not identical." + print("Label list lengths not identical.") sys.exit() dataObjList = [] - for filename in os.listdir(os.getcwd()): + for filename in os.listdir("/PsrPopPy/"): for extn in extnList: if filename.endswith(extn): # read in population file to population self.pop + # filename = "/PsrPopPy/populate.model" try: f = open(filename, 'rb') except IOError: - print "Could not open file {0}.".format(filename) + print("Could not open file {0}.".format(filename)) sys.exit() - pop = cPickle.load(f) + pop = pickle.load(f) f.close() # for each of those labels, get the data from population into lists @@ -104,19 +111,19 @@ def makeDataObj(frac=None, extnList=['.model', '.results']): dataArray[3][npsr] = psr.galCoords[1] dataArray[4][npsr] = psr.galCoords[2] dataArray[5][npsr] = psr.width_degree - dataArray[6][npsr] = psr.alpha - dataArray[7][npsr] = psr.rho - dataArray[8][npsr] = psr.spindex - dataArray[9][npsr] = psr.s_1400() - dataArray[10][npsr] = psr.gl - dataArray[11][npsr] = psr.gb - dataArray[12][npsr] = psr.dtrue - dataArray[13][npsr] = psr.r0 - dataArray[14][npsr] = npsr + dataArray[6][npsr] = psr.pdot + dataArray[7][npsr] = (psr.period / (2 * psr.pdot) / 31556952 ) if psr.pdot else 0 + dataArray[8][npsr] = (np.sqrt(psr.period * psr.pdot) * 3.2e16) if psr.pdot else 0 + dataArray[9][npsr] = psr.spindex + dataArray[10][npsr] = psr.s_1400() + dataArray[11][npsr] = psr.gl + dataArray[12][npsr] = psr.gb + dataArray[13][npsr] = psr.dtrue + dataArray[14][npsr] = psr.r0 + dataArray[15][npsr] = npsr npsr+=1 - dataObjList.append( - DataObj(filename, textlabels, axislabels, dataArray, pop.size())) + dataObjList.append(DataObj(filename, textlabels, axislabels, dataArray, pop.size())) if len(dataObjList) == 0: raise ViewException('No files found matching the extensions') @@ -131,179 +138,182 @@ class VisualizeFrame(tk.Frame): def __init__(self, dataObjList): self.dataObjList = dataObjList - tk.Frame.__init__(self, None, -1) - - self.colour_list = ['b.', 'r.', 'g.', 'c.', 'm.', 'y.', 'k.'] + tk.Frame.__init__(self, master=root) + + self.colour_list = ['r.', 'g.', 'y.', 'm.', 'c.', 'b.', 'k.'] + self.xIndex = tk.IntVar() + self.yIndex = tk.IntVar(value=1) + self.logx_var = tk.BooleanVar() + self.logy_var = tk.BooleanVar() + self.grid_var = tk.BooleanVar() + self.hist_var = tk.BooleanVar() + self.gaus_var = tk.IntVar() + self.plot_var = tk.BooleanVar(value=True) self.create_main_panel() - def create_main_panel(self): - self.panel = tk.Panel(self) + self.panel = tk.Frame(self) + self.panel.pack(fill=tk.BOTH, expand=True) # Pack the panel to fill the entire space of the master window + self.dpi = 100 - self.fig = Figure((5., 5.), dpi = self.dpi) - self.canvas = FigCanvas(self.panel, -1, self.fig) + self.fig = Figure((5., 5.), dpi=self.dpi) + self.canvas = FigureCanvasTkAgg(self.fig, master=self.panel) + self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=True) # Pack the canvas to the top of the panel, filling both x and y directions + + self.toolbar = NavigationToolbar(self.canvas, self.panel) + # self.toolbar.update() + # self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=True) # Pack the canvas again (it's the same canvas but now with toolbar), to the top of the panel, filling both x and y directions self.axes = self.fig.add_subplot(111) - self.drawbutton = tk.Button(self.panel, -1, "Plot") - self.Bind(tk.EVT_BUTTON, self.on_draw_button, self.drawbutton) + self.buttons1 = tk.Frame(self.panel) + self.buttons1.pack(side=tk.TOP) - self.logx = tk.CheckBox(self.panel, -1, - "log X", - style=tk.ALIGN_RIGHT) - self.Bind(tk.EVT_CHECKBOX, self.on_logx, self.logx) + self.buttons2 = tk.Frame(self.panel) + self.buttons2.pack(side=tk.TOP) - self.logy = tk.CheckBox(self.panel, -1, - "log Y", - style=tk.ALIGN_RIGHT) - self.Bind(tk.EVT_CHECKBOX, self.on_logy, self.logy) + tk.Button(self.buttons1, text="Plot", command=self.on_draw_button).pack(side=tk.LEFT, padx=5, pady=5) + tk.Checkbutton(self.buttons1, text="logX", variable=self.logx_var).pack(side=tk.LEFT, padx=5, pady=5) + tk.Checkbutton(self.buttons1, text="logY", variable=self.logy_var).pack(side=tk.LEFT, padx=5, pady=5) + tk.Checkbutton(self.buttons1, text="grid", variable=self.grid_var).pack(side=tk.LEFT, padx=5, pady=5) - self.grid = tk.CheckBox(self.panel, -1, - "grid", - style=tk.ALIGN_RIGHT) - self.Bind(tk.EVT_CHECKBOX, self.on_grid, self.grid) - - modelList = [d.name for d in self.dataObjList] + tk.Checkbutton(self.buttons2, text="hist", variable=self.hist_var, command=self._hist).pack(side=tk.LEFT, padx=5, pady=5) + self.gauss1 = tk.Radiobutton(self.buttons2, text="gauss1", variable=self.gaus_var, value=0, state=tk.DISABLED) + self.gauss1.pack(side=tk.LEFT, padx=5, pady=5) + self.gauss2 = tk.Radiobutton(self.buttons2, text="gauss2", variable=self.gaus_var, value=1, state=tk.DISABLED) + self.gauss2.pack(side=tk.LEFT, padx=5, pady=5) # create list of population models. Set all "on" by default - self.modelCheckList = tk.CheckListBox(self.panel, -1, - choices=modelList, - style=tk.ALIGN_RIGHT) - self.modelCheckList.SetChecked(range(len(modelList))) - - # Create the navigation toolbar, tied to the canvas - # - self.toolbar = NavigationToolbar(self.canvas) - - # - # Layout with box sizers - # - - self.vbox = tk.BoxSizer(wx.VERTICAL) - self.v_buttonbox_1 = tk.BoxSizer(wx.VERTICAL) - self.v_buttonbox_2 = tk.BoxSizer(wx.VERTICAL) - self.hpltbox = tk.BoxSizer(wx.HORIZONTAL) # for plot and plot selection - self.htoolbox = tk.BoxSizer(wx.HORIZONTAL) # for drawing and log toggles - - # fill the top box first, radio buttons in a vbox, - # next to the canvas - self.radioBoxX = tk.RadioBox(self.panel, 1, 'X Axis', - choices = self.dataObjList[0].textlabels, - majorDimension=1, - style = tk.RA_SPECIFY_COLS) - - self.radioBoxY = tk.RadioBox(self.panel, 2, 'Y Axis', - choices = self.dataObjList[0].textlabels, - majorDimension=1, - style = tk.RA_SPECIFY_COLS) - - self.xIndex=0 - self.yIndex=0 - - # event for radio box with ID 1 - tk.EVT_RADIOBOX(self.panel, 1, self.onXRadioClick) - # event for radiobox with ID 2 - tk.EVT_RADIOBOX(self.panel, 2, self.onYRadioClick) - - - # top horizontal panel - canvas and radio boxes - self.hpltbox.Add(self.radioBoxX) - self.hpltbox.Add(self.radioBoxY) - self.hpltbox.Add(self.canvas, 1, tk.LEFT | wx.TOP | wx.GROW) - self.hpltbox.Add(self.modelCheckList, 0, border=3) - - # add the matplotlib toolbar - self.vbox.Add(self.hpltbox) - self.vbox.Add(self.toolbar, 0, tk.EXPAND) - self.vbox.AddSpacer(10) - - - # bottom horizontal panel - check buttons and plot button - flags = tk.ALIGN_LEFT | wx.ALL | wx.ALIGN_CENTER_VERTICAL - self.htoolbox.Add(self.drawbutton, 0, border=3, flag=flags) - self.htoolbox.Add(self.logx, 0, border=3, flag=flags) - self.htoolbox.Add(self.logy, 0, border=3, flag=flags) - self.htoolbox.AddSpacer(20) - self.htoolbox.Add(self.grid, 0, border=3, flag=flags) + modelList = [d.name for d in self.dataObjList] + self.modelCheckList = tk.Listbox(self.panel, selectmode=tk.MULTIPLE) + for model in modelList: + self.modelCheckList.insert(tk.END, model) + self.modelCheckList.select_set(0, tk.END) # Select all items by default + self.modelCheckList.pack(side=tk.LEFT, padx=5, pady=5, fill=tk.Y) - self.vbox.Add(self.htoolbox, 0, flag = tk.ALIGN_LEFT | wx.TOP) + # create LabelFrames for holding Radiobuttons + self.radioBoxX = tk.LabelFrame(self.panel, text='X Axis') + self.radioBoxX.pack(side=tk.LEFT, padx=5, pady=5, fill=tk.X) - self.panel.SetSizer(self.vbox) - self.vbox.Fit(self) + self.radioBoxY = tk.LabelFrame(self.panel, text='Y Axis') + self.radioBoxY.pack(side=tk.LEFT, padx=5, pady=5, fill=tk.Y) - def onRadioClick(self, event): - radioBox = event.GetEventObject() - print event.GetId() - print radioBox + # Create and pack Radiobuttons for pulsar props inside the LabelFrames + labels = self.dataObjList[0].textlabels + for i in range(len(labels)): + tk.Radiobutton(self.radioBoxX, text=labels[i], variable=self.xIndex, value=i).pack(anchor=tk.W) + tk.Radiobutton(self.radioBoxY, text=labels[i], variable=self.yIndex, value=i).pack(anchor=tk.W) - def onXRadioClick(self, event): - # get x index of selected button - radioBox = event.GetEventObject() - self.xIndex=radioBox.GetSelection() + # top horizontal panel - canvas and radio boxes + self.radioBoxX.pack(side=tk.LEFT, padx=5, pady=5, fill=tk.Y) + self.radioBoxY.pack(side=tk.LEFT, padx=5, pady=5, fill=tk.Y) + self.canvas.get_tk_widget().pack(side=tk.LEFT, padx=5, pady=5, fill=tk.BOTH, expand=True) + # self.modelCheckList.pack(side=tk.LEFT, padx=5, pady=5, fill=tk.Y) - def onYRadioClick(self, event): - # gets yIndex of selected button - radioBox = event.GetEventObject() - self.yIndex=radioBox.GetSelection() + # add the matplotlib toolbar + self.toolbar.pack(side=tk.TOP, fill=tk.X) def draw_figure(self): self.axes.clear() - for dataObjIndex in self.modelCheckList.GetChecked(): - try: - self.axes.plot(self.dataObjList[dataObjIndex].dataArray[self.xIndex], - self.dataObjList[dataObjIndex].dataArray[self.yIndex], - self.colour_list[dataObjIndex], - label=self.dataObjList[dataObjIndex].name) - except IndexError: - self.axes.plot(self.dataObjList[dataObjIndex].dataArray[self.xIndex], - self.dataObjList[dataObjIndex].dataArray[self.yIndex], - label=self.dataObjList[dataObjIndex].name) - #self.axes.legend(loc='upper center', bbox_to_anchor=(0.5,-0.05), - # ncol=len(self.modelCheckList.GetChecked()), - # ) - - if len(self.modelCheckList.GetChecked())>0: - self.axes.set_xlabel(self.dataObjList[0].axislabels[self.xIndex], + if not self.plot_var.get(): + self.axes.text(0.5, 0.5, "Not plottable") + self.plot_var.set(True) + return + + for dataObjIndex in self.modelCheckList.curselection(): + dataX = self.dataObjList[dataObjIndex].dataArray[self.xIndex.get()] + dataY = self.dataObjList[dataObjIndex].dataArray[self.yIndex.get()] + color = self.colour_list[dataObjIndex] + + if not self.hist_var.get(): + self.axes.plot(dataX, dataY, color, + label=self.dataObjList[dataObjIndex].name) + else: + # try: + hist, bins_edges, _ = \ + self.axes.hist(dataX, color=color[0], + label=self.dataObjList[dataObjIndex].name, + bins=100, density=True, alpha=0.5, edgecolor='k') + # except: + # self.plot_var.set(False) + # self.draw_figure() + # return + + if not self.gaus_var.get(): + # try: + μ, σ = norm.fit(dataX) + xmin, xmax = self.axes.get_xlim() + # warnings.simplefilter('error') + # warnings.warn("Not allowed") + x = np.linspace(xmin, xmax, 1000) + p = norm.pdf(x, μ, σ) + self.axes.plot(x, p, color[0]) + # except: + # self.plot_var.set(False) + # self.draw_figure() + # return + else: + # try: + bins_centers = (bins_edges[:-1] + bins_edges[1:]) / 2 + peak_indices = argrelmax(hist)[0] + peak_values = hist[peak_indices] + top_two_peaks = np.argsort(peak_values)[::-1][:2] + μ1, μ2 = bins_centers[peak_indices[top_two_peaks]] + + def gaussian_mixture(x, A1, μ1, σ1, A2, μ2, σ2): + return (A1 * norm.pdf(x, μ1, σ1) + + A2 * norm.pdf(x, μ2, σ2)) + + initial_guess = [1, μ1, 1, 1, μ2, 1] + params, cov = curve_fit(gaussian_mixture, bins_centers, hist, p0=initial_guess) + # warnings.simplefilter('error') + # warnings.warn("Not allowed") + x = np.linspace(min(dataX), max(dataX), 1000) + self.axes.plot(x, gaussian_mixture(x, *params), color[0]) + # except: + # self.plot_var.set(False) + # self.draw_figure() + # return + + # self.axes.set_title(f"μ={μ:.2f}, σ={σ:.2f}") + + self.axes.legend(loc='upper center', bbox_to_anchor=(0.5, -0.075), + ncol=len(self.modelCheckList.curselection()), + ) + + if len(self.modelCheckList.curselection()) > 0: + self.axes.set_xlabel(self.dataObjList[0].axislabels[self.xIndex.get()], fontsize=10) - self.axes.set_ylabel(self.dataObjList[0].axislabels[self.yIndex], + self.axes.set_ylabel(self.dataObjList[0].axislabels[self.yIndex.get()], fontsize=10) - for label in self.axes.get_xticklabels(): - label.set_fontsize(7) - for label in self.axes.get_yticklabels(): - label.set_fontsize(7) - - if self.logx.IsChecked(): + if self.logx_var.get(): self.axes.set_xscale('log') - for label in self.axes.get_xticklabels(): - label.set_fontsize(7) - if self.logy.IsChecked(): + if self.logy_var.get(): self.axes.set_yscale('log') - for label in self.axes.get_yticklabels(): - label.set_fontsize(7) - self.axes.grid(self.grid.IsChecked()) - self.canvas.draw() + # for label in self.axes.get_xticklabels(): + # label.set_fontsize(7) + # for label in self.axes.get_yticklabels(): + # label.set_fontsize(7) - def on_draw_button(self, event): - # redraw the canvas - self.draw_figure() - - def on_logx(self, event): - # placeholder functions in case I want to do something else - # decided I don't want plot to redraw until the button is - # clicke, rather than on selection of log axis - pass + self.axes.grid(self.grid_var.get()) + self.canvas.draw() - def on_logy(self, event): - pass + def _hist(self): + if self.hist_var.get(): + self.gauss1.config(state=tk.NORMAL) + self.gauss2.config(state=tk.NORMAL) + else: + self.gauss1.config(state=tk.DISABLED) + self.gauss2.config(state=tk.DISABLED) - def on_grid(self, event): + def on_draw_button(self): + # redraw the canvas self.draw_figure() - if __name__ == '__main__': """ 'Main' function for calling from command line""" parser = argparse.ArgumentParser(description='Visualize a population object') @@ -318,10 +328,10 @@ def on_grid(self, event): help='extension(s) to look for when finding population models') args = parser.parse_args() - dataObj = makeDataObj(frac=args.frac, extnList=args.extn) + dataObj = makeDataObj(frac=None, extnList=['.results', '.model']) - app = tk.PySimpleApp() - app.frame = VisualizeFrame(dataObj) - app.frame.Show() - app.MainLoop() + root = tk.Tk() + frame = VisualizeFrame(dataObj) + frame.pack(fill=tk.BOTH, expand=True) + root.mainloop() diff --git a/bin/old_viewing_software/tkView_old.py b/bin/old_viewing_software/tkView_old.py new file mode 100644 index 0000000..455a473 --- /dev/null +++ b/bin/old_viewing_software/tkView_old.py @@ -0,0 +1,328 @@ +#!/usr/bin/env python + +import sys +import os +import argparse +import random + +import pickle +import numpy as np + +import tkinter as tk +import matplotlib +matplotlib.use('TKAgg') +from matplotlib.figure import Figure +from matplotlib.backends.backend_tkagg import \ + FigureCanvasTkAgg, \ + NavigationToolbar2Tk as NavigationToolbar + +class ViewException(Exception): + pass + +class DataObj: + def __init__(self, name, textlabels, axislabels, dataarray, size): + self.name = name + self.textlabels = textlabels + self.axislabels = axislabels + self.dataArray = dataarray + self.count = size + +def makeDataObj(frac=None, extnList=['.model', '.results']): + """ + Function to get different parameters of a population file + + """ + textlabels = ['Period', + 'DM', + 'Gal X', + 'Gal Y', + 'Gal Z', + 'W', + 'alpha', + 'rho', + 'SI', + 'S1400', + 'gl', + 'gb', + 'D', + 'r0', + 'n'] + axislabels = ['Period (ms)', + 'DM (cm^-3 pc)', + 'X (kpc)', + 'Y (kpc)', + 'Z (kpc)', + 'Width (degrees)', + 'alpha (deg)', + 'rho (deg)', + 'Spectral Index', + 'S1400 (mJy)', + 'Galactic Longitude (degrees)', + 'Galactic Latitude (degrees)', + 'Distance (kpc)', + 'GalacticRadius (kpc)', + 'Array Index'] + + if len(textlabels) != len(axislabels): + print("Label list lengths not identical.") + sys.exit() + + dataObjList = [] + + for filename in os.listdir(os.getcwd()): + for extn in extnList: + if filename.endswith(extn): + # read in population file to population self.pop + try: + f = open(filename, 'rb') + except IOError: + print("Could not open file {0}.".format(filename)) + sys.exit() + + pop = pickle.load(f) + f.close() + + # for each of those labels, get the data from population into lists + # I suppose most sense is to read it into arrays + + # create numpy array of right size + dataArray = np.zeros((len(textlabels), pop.size()), float) + # loop over pulsars and fill array + npsr = 0 + + # going to throw in a factor to reduce the number of plotted pulsars + # to improve speed + # make this an option + if frac is None or frac > 1.0: + frac = 1.0 + + for psr in pop.population: + if random.random() < frac: + dataArray[0][npsr] = psr.period + dataArray[1][npsr] = psr.dm + dataArray[2][npsr] = psr.galCoords[0] + dataArray[3][npsr] = psr.galCoords[1] + dataArray[4][npsr] = psr.galCoords[2] + dataArray[5][npsr] = psr.width_degree + dataArray[6][npsr] = psr.alpha + dataArray[7][npsr] = psr.rho + dataArray[8][npsr] = psr.spindex + dataArray[9][npsr] = psr.s_1400() + dataArray[10][npsr] = psr.gl + dataArray[11][npsr] = psr.gb + dataArray[12][npsr] = psr.dtrue + dataArray[13][npsr] = psr.r0 + dataArray[14][npsr] = npsr + npsr+=1 + + dataObjList.append( + DataObj(filename, textlabels, axislabels, dataArray, pop.size())) + + if len(dataObjList) == 0: + raise ViewException('No files found matching the extensions') + + # sort list from largest -> smallest so plotting is always done that way + dataObjList.sort(key=lambda x: x.count, reverse=True) + + return dataObjList + + +class VisualizeFrame(tk.Frame): + def __init__(self, dataObjList): + + self.dataObjList = dataObjList + tk.Frame.__init__(self, master=tk.Tk()) + + self.colour_list = ['b.', 'r.', 'g.', 'c.', 'm.', 'y.', 'k.'] + self.xIndex = tk.IntVar() + self.yIndex = tk.IntVar() + self.create_main_panel() + + + + def create_main_panel(self): + self.panel = tk.Frame(self) + self.dpi = 100 + self.fig = Figure((5., 5.), dpi = self.dpi) + self.canvas = FigureCanvasTkAgg(self.fig, master=self.panel) + + self.axes = self.fig.add_subplot(111) + + self.drawbutton = tk.Button(self.panel, text="Plot") + self.drawbutton.bind("", self.on_draw_button) + + self.logx = tk.Checkbutton(self.panel, text="log X", variable=tk.IntVar()) + self.logx.bind("", self.on_logx) + + self.logy = tk.Checkbutton(self.panel, text="log Y", variable=tk.IntVar()) + self.logy.bind("", self.on_logy) + + self.grid = tk.Checkbutton(self.panel, text="grid", variable=tk.IntVar()) + self.grid.bind("", self.on_grid) + + modelList = [d.name for d in self.dataObjList] + + # create list of population models. Set all "on" by default + self.modelCheckList = tk.Listbox(self.panel, selectmode=tk.MULTIPLE) + for model in modelList: + self.modelCheckList.insert(tk.END, model) + self.modelCheckList.select_set(0, tk.END) # Select all items by default + + # Create the navigation toolbar, tied to the canvas + self.toolbar = NavigationToolbar(self.canvas) + + # + # Layout with box sizers + # + + self.vbox = tk.Frame(self) + self.v_buttonbox_1 = tk.Frame(self.vbox) + self.v_buttonbox_2 = tk.Frame(self.vbox) + self.hpltbox = tk.Frame(self.vbox) # for plot and plot selection + self.htoolbox = tk.Frame(self.vbox) # for drawing and log toggles + + self.v_buttonbox_1.pack() + self.v_buttonbox_2.pack() + self.hpltbox.pack() + self.htoolbox.pack() + + # fill the top box first, radio buttons in a vbox, + # next to the canvas + self.radioBoxX = tk.LabelFrame(self.panel, text='X Axis') + self.radioBoxY = tk.LabelFrame(self.panel, text='Y Axis') + + # Create and pack Radiobuttons inside the LabelFrames + for label in self.dataObjList[0].textlabels: + tk.Radiobutton(self.radioBoxX, text=label, variable=self.xIndex, value=label).pack(anchor=tk.W) + + for label in self.dataObjList[0].textlabels: + tk.Radiobutton(self.radioBoxY, text=label, variable=self.yIndex, value=label).pack(anchor=tk.W) + + self.xIndex=0 + self.yIndex=0 + + # event for radio box with ID 1 + self.radioBoxX.bind("", self.onXRadioClick) + # event for radiobox with ID 2 + self.radioBoxY.bind("", self.onYRadioClick) + + + # top horizontal panel - canvas and radio boxes + self.hpltbox.Add(self.radioBoxX) + self.hpltbox.Add(self.radioBoxY) + self.hpltbox.Add(self.canvas, 1, tk.LEFT | tk.TOP | tk.GROW) + self.hpltbox.Add(self.modelCheckList, 0, border=3) + + # add the matplotlib toolbar + self.vbox.Add(self.hpltbox) + self.vbox.Add(self.toolbar, 0, tk.EXPAND) + self.vbox.AddSpacer(10) + + + # bottom horizontal panel - check buttons and plot button + flags = tk.ALIGN_LEFT | tk.ALL | tk.ALIGN_CENTER_VERTICAL + self.htoolbox.Add(self.drawbutton, 0, border=3, flag=flags) + self.htoolbox.Add(self.logx, 0, border=3, flag=flags) + self.htoolbox.Add(self.logy, 0, border=3, flag=flags) + self.htoolbox.AddSpacer(20) + self.htoolbox.Add(self.grid, 0, border=3, flag=flags) + + self.vbox.Add(self.htoolbox, 0, flag = tk.ALIGN_LEFT | tk.TOP) + + self.panel.SetSizer(self.vbox) + self.vbox.Fit(self) + + def onRadioClick(self, event): + radioBox = event.GetEventObject() + print(event.GetId()) + print(radioBox) + + def onXRadioClick(self, event): + # get x index of selected button + radioBox = event.GetEventObject() + self.xIndex=radioBox.GetSelection() + + def onYRadioClick(self, event): + # gets yIndex of selected button + radioBox = event.GetEventObject() + self.yIndex=radioBox.GetSelection() + + def draw_figure(self): + self.axes.clear() + for dataObjIndex in self.modelCheckList.GetChecked(): + try: + self.axes.plot(self.dataObjList[dataObjIndex].dataArray[self.xIndex], + self.dataObjList[dataObjIndex].dataArray[self.yIndex], + self.colour_list[dataObjIndex], + label=self.dataObjList[dataObjIndex].name) + except IndexError: + self.axes.plot(self.dataObjList[dataObjIndex].dataArray[self.xIndex], + self.dataObjList[dataObjIndex].dataArray[self.yIndex], + label=self.dataObjList[dataObjIndex].name) + #self.axes.legend(loc='upper center', bbox_to_anchor=(0.5,-0.05), + # ncol=len(self.modelCheckList.GetChecked()), + # ) + + if len(self.modelCheckList.GetChecked())>0: + self.axes.set_xlabel(self.dataObjList[0].axislabels[self.xIndex], + fontsize=10) + self.axes.set_ylabel(self.dataObjList[0].axislabels[self.yIndex], + fontsize=10) + + for label in self.axes.get_xticklabels(): + label.set_fontsize(7) + for label in self.axes.get_yticklabels(): + label.set_fontsize(7) + + if self.logx.IsChecked(): + self.axes.set_xscale('log') + for label in self.axes.get_xticklabels(): + label.set_fontsize(7) + + if self.logy.IsChecked(): + self.axes.set_yscale('log') + for label in self.axes.get_yticklabels(): + label.set_fontsize(7) + + self.axes.grid(self.grid.IsChecked()) + self.canvas.draw() + + def on_draw_button(self, event): + # redraw the canvas + self.draw_figure() + + def on_logx(self, event): + # placeholder functions in case I want to do something else + # decided I don't want plot to redraw until the button is + # clicke, rather than on selection of log axis + pass + + def on_logy(self, event): + pass + + def on_grid(self, event): + self.draw_figure() + + +if __name__ == '__main__': + """ 'Main' function for calling from command line""" + parser = argparse.ArgumentParser(description='Visualize a population object') + parser.add_argument('-f', metavar='fname', default='populate.model', + help='file containing population model (def="populate.model")') + + parser.add_argument('-frac', nargs=1, type=float, default=None, + help='plot only this fraction of pulsars') + + parser.add_argument('-extn', nargs='+', type=str, + default=['.results', '.model'], + help='extension(s) to look for when finding population models') + args = parser.parse_args() + + dataObj = makeDataObj(frac=args.frac, extnList=args.extn) + + app = tk.Tk() + app.frame = VisualizeFrame(dataObj) + app.frame.Show() + app.pack() + app.MainLoop() + diff --git a/bin/old_viewing_software/visualize.py b/bin/old_viewing_software/visualize.py index 5edc22d..0a32cba 100644 --- a/bin/old_viewing_software/visualize.py +++ b/bin/old_viewing_software/visualize.py @@ -5,12 +5,17 @@ import math import random -import cPickle +import pickle +import matplotlib +matplotlib.use('TkAgg') import matplotlib.pyplot as plt from matplotlib.widgets import Button, RadioButtons, CheckButtons import numpy as np +from scipy.stats import norm +from scipy.optimize import curve_fit +from scipy.signal import argrelmax from psrpoppy.population import Population from psrpoppy.pulsar import Pulsar @@ -22,21 +27,21 @@ class Visualize: using the matplotlib module. """ - def __init__(self, popfile='populate.model', frac=None): + def __init__(self, popfile, frac): """Initialise the Visualize object.""" # read in population file to population self.pop try: f = open(popfile, 'rb') except IOError: - print "Could not open file {0}.".format(popfile) + print(f"Could not open file {popfile}.") sys.exit() - self.pop = cPickle.load(f) + self.pop = pickle.load(f) f.close() # print out the population parameters - print self.pop + print(self.pop) self.textlabels = ['Period', 'DM', @@ -44,33 +49,33 @@ def __init__(self, popfile='populate.model', frac=None): 'Gal Y', 'Gal Z', 'W', - 'alpha', - 'rho', + 'Age', + 'Pdot', 'SI', 'S1400', 'gl', 'gb', 'D', 'r0', - 'n'] + 'B'] self.axislabels = ['Period (ms)', 'DM (cm^-3 pc)', 'X (kpc)', 'Y (kpc)', 'Z (kpc)', 'Width (degrees)', - 'alpha (deg)', - 'rho (deg)', + 'Age [P/2Pdot] (years)', + 'Pdot', 'Spectral Index', 'S1400 (mJy)', 'Galactic Longitude (degrees)', 'Galactic Latitude (degrees)', 'Distance (kpc)', 'GalacticRadius (kpc)', - 'Array Index'] + 'B [(PPdot)^(1/2)] (G)'] if len(self.textlabels) != len(self.axislabels): - print "Label list lengths not identical." + print("Label list lengths not identical.") sys.exit() # for each of those labels, get the data from population into lists @@ -87,6 +92,8 @@ def __init__(self, popfile='populate.model', frac=None): if frac is None or frac > 1.0: frac = 1.0 + # print(self.pop.population[0].__dir__()) + for psr in self.pop.population: if random.random() < frac: dataArray[0][npsr] = psr.period @@ -95,15 +102,15 @@ def __init__(self, popfile='populate.model', frac=None): dataArray[3][npsr] = psr.galCoords[1] dataArray[4][npsr] = psr.galCoords[2] dataArray[5][npsr] = psr.width_degree - dataArray[6][npsr] = psr.alpha - dataArray[7][npsr] = psr.rho + dataArray[6][npsr] = psr.period / (2 * psr.pdot) / 31556952 + dataArray[7][npsr] = psr.pdot dataArray[8][npsr] = psr.spindex dataArray[9][npsr] = psr.s_1400() dataArray[10][npsr] = psr.gl dataArray[11][npsr] = psr.gb dataArray[12][npsr] = psr.dtrue dataArray[13][npsr] = psr.r0 - dataArray[14][npsr] = npsr + dataArray[14][npsr] = np.sqrt(psr.period * psr.pdot) * 3.2e16 npsr+=1 self.dataArray = dataArray @@ -113,70 +120,55 @@ def __init__(self, popfile='populate.model', frac=None): def display(self): """Method to create the plotting window and fill it with buttons.""" + # Create matplotlib window dressing - self.fig = plt.figure(figsize=(11,9), facecolor='lightgoldenrodyellow') - #self.fig = plt.figure(facecolor='lightgoldenrodyellow') - self.fig.canvas.set_window_title('PyPop: Visualization') + self.fig = plt.figure(figsize=(10, 10), facecolor='lightgoldenrodyellow') + # self.fig.canvas.set_window_title('PyPop: Visualization') self.fig.clear() - # put some buttons on the figure for the plot axes - - buttonW = 0.055 - buttonH = 0.045 # default button sizes + # give titles + self.fig.text(0.01, 0.97, 'PLOT SELECTION') + self.fig.text(0.01, 0.80, 'X Axis') + self.fig.text(0.01, 0.50, 'Y Axis') - # arguments are button left, bottom, width, height and label + # arguments are [left, bottom, width, height] and label + + # radio buttons for x and y axis values + radioXAxis = RadioButtons(plt.axes([0.070, 0.650, 0.150, 0.300]), self.textlabels, + active=0, activecolor='blue') + self.xindex = 0 - # check box to do log scale - buttonXLog = CheckButtons(plt.axes([0.001, .3, 1.5*buttonW, 1.5*buttonH]), - ['log x'], + radioYAxis = RadioButtons(plt.axes([0.070, 0.350, 0.150, 0.300]), self.textlabels, + active=1, activecolor='blue') + self.yindex = 1 + + # checkboxes to do log scale + buttonXLog = CheckButtons(plt.axes([0.025, 0.225, 0.075, 0.050]), ['log x'], actives=[False]) self.xlog = False - buttonYLog = CheckButtons(plt.axes([0.001+2.5*buttonW, .3, - 1.5*buttonW, 1.5*buttonH]), - ['log y'], + buttonYLog = CheckButtons(plt.axes([0.125, 0.225, 0.075, 0.050]), ['log y'], actives=[False]) self.ylog = False - # button to do plot - buttonPlot = Button(plt.axes([0.001+2.0*buttonW, 0.2, 1.5*buttonW, buttonH]), - 'Plot') - - self.fig.text(0.003, 0.98, 'PLOT SELECTION') - self.fig.text(0.003, 0.95, 'X Axis') - self.fig.text(0.003, 0.65, 'Y Axis') - - # radio buttons for x and y axis values - radioXAxis = RadioButtons(plt.axes([0.001, - 0.72, - len(self.textlabels)*0.02, - len(self.textlabels)*0.02]), - self.textlabels, - active=0, - activecolor='blue') - self.xindex = 0 + buttonPlot = Button (plt.axes([0.075, 0.150, 0.075, 0.050]), 'Plot') - radioYAxis = RadioButtons(plt.axes([0.001, - 0.42, - len(self.textlabels)*0.02, - len(self.textlabels)*0.02]), - self.textlabels, - active=1, - activecolor='blue') - self.yindex = 1 + # button to do histogram plot + buttonHist = Button (plt.axes([0.075, 0.050, 0.075, 0.050]), 'Hist') # if radio button is clicked, change active X/Y index radioXAxis.on_clicked(self._radioXClick) radioYAxis.on_clicked(self._radioYClick) - # add callback to the plot button(s?) - buttonPlot.on_clicked(self._scatterplot) - # callbacks for the log plot switches buttonXLog.on_clicked(self._logClick) buttonYLog.on_clicked(self._logClick) + # add callback to the plot button(s?) + buttonPlot.on_clicked(self._scatterplot) + buttonHist.on_clicked(self._histplot) + plt.show() def _logClick(self, label): @@ -186,7 +178,7 @@ def _logClick(self, label): elif label == 'log y': self.ylog = not self.ylog else: - print "Weird log label!" + print("Weird log label!") sys.exit() def _radioXClick(self, label): @@ -196,52 +188,113 @@ def _radioYClick(self, label): self.yindex = self.textlabels.index(label) def _scatterplot(self, event): - # if there's already a scatter plot, delete it! - try: - self.fig.delaxes(self.scatterplt) - except: - pass + # if there's already a scatter/histogram plot, delete it! + try: self.fig.delaxes(self.histplot) + except: pass + + try: self.fig.delaxes(self.scatterplt) + except: pass # define axis position and dimensions - self.scatterplt = self.fig.add_axes([0.3, 0.15, 0.65, 0.8]) + self.scatterplt = self.fig.add_axes([0.32, 0.15, 0.65, 0.8]) # axis labels self.scatterplt.set_xlabel(self.axislabels[self.xindex]) self.scatterplt.set_ylabel(self.axislabels[self.yindex]) - # plot stuff. + # plot stuff self.scatterplt.plot(self.dataArray[self.xindex], self.dataArray[self.yindex], 'r.') - + # if log switches are on, do a log plot if self.xlog: try: self.scatterplt.set_xscale('log') except ValueError: - print "matplotlib refuses to do a log plot of the X data!" + print("matplotlib refuses to do a log plot of the X data!") pass if self.ylog: try: self.scatterplt.set_yscale('log') except ValueError: - print "matplotlib refuses to do a log plot of the Y data!" + print("matplotlib refuses to do a log plot of the Y data!") pass # finally, display the plot plt.show() + def _histplot(self, event): + # if there's already a scatter/histogram plot, delete it! + try: self.fig.delaxes(self.scatterplt) + except: pass + + try: self.fig.delaxes(self.histplot) + except: pass + + # define axis position and dimensions + self.histplot = self.fig.add_axes([0.32, 0.15, 0.65, 0.8]) + + # axis labels + self.histplot.set_xlabel(self.axislabels[self.xindex]) + + # define data + data = self.dataArray[self.xindex] + + if self.xindex in [2, 3, 10, 12]: + # plotting two gaussians for better fit + hist, bins_edges, _ = self.histplot.hist(data, bins=100, density=True, alpha=0.5, edgecolor="black") + + # Calculate bin centers + bins_centers = (bins_edges[:-1] + bins_edges[1:]) / 2 + + # Find peaks in histogram + peak_indices = argrelmax(hist)[0] # Find local maximums + peak_values = hist[peak_indices] + top_two_peaks = np.argsort(peak_values)[::-1][:2] # Select the two highest peaks + μ1, μ2 = bins_centers[peak_indices[top_two_peaks]] + + # Define the model function as a mixture of two Gaussians + # with corresponding parameters (amplitude A, mean μ, standard deviation σ) + def gaussian_mixture(x, A1, μ1, σ1, A2, μ2, σ2): + return (A1 * norm.pdf(x, μ1, σ1) + + A2 * norm.pdf(x, μ2, σ2)) + + # Initial guess for the parameters + initial_guess = [1, μ1, 1, 1, μ2, 1] + + # Fit the model to the data [params (A1, μ1, σ1, A2, μ2, σ2), covariance matrix] + params, cov = curve_fit(gaussian_mixture, bins_centers, hist, p0=initial_guess) + + # Plot the fitted curve [*params --> unpacking params] + x = np.linspace(min(data), max(data), 1000) + plt.plot(x, gaussian_mixture(x, *params), "r") + plt.title("A1={0:.2f}, μ1={1:.2f}, σ1={2:.2f}, A2={3:.2f}, μ2={4:.2f}, σ2={5:.2f}".format(*params)) + + else: + # plot histogram along with gaussian + self.histplot.hist(data, bins=100, density=True, alpha=0.5, edgecolor="black") + + μ, σ = norm.fit(data) + xmin, xmax = plt.xlim() + x = np.linspace(xmin, xmax, 1000) + p = norm.pdf(x, μ, σ) + plt.plot(x, p, "k") + plt.title(f"μ={μ:.2f}, σ={σ:.2f}") + + plt.show() + if __name__ == '__main__': """ 'Main' function for calling from command line""" parser = argparse.ArgumentParser(description='Visualize a population object') - parser.add_argument('-f', metavar='fname', default='populate.model', - help='file containing population model (def="populate.model")') + parser.add_argument('-f', metavar='filename', required=True, + help='file containing population model') - parser.add_argument('-frac', nargs=1, type=float, default=None, - help='plot only this fraction of pulsars') + parser.add_argument('-frac', metavar='fraction', type=float, default=None, + help='plot only this fraction of pulsars (default=None)') args = parser.parse_args() - v = Visualize(popfile = args.f, frac=args.frac) + v = Visualize(popfile=args.f, frac=args.frac) v.display() diff --git a/bin/popheader b/bin/popheader index 4f3176c..5bb1bf0 100755 --- a/bin/popheader +++ b/bin/popheader @@ -1,21 +1,21 @@ #!/usr/bin/env python -import cPickle +import pickle import sys # check there's one argument if (len(sys.argv)!=2): - print "Usage: popheader [pop model]" + print("Usage: popheader [pop model]") sys.exit() if sys.argv[1] == "-h" or sys.argv[1] == "--help": - print "Usage: popheader [pop model]" + print("Usage: popheader [pop model]") sys.exit() # load the model with open(sys.argv[1], 'rb') as f: - pop = cPickle.load(f) + pop = pickle.load(f) -print pop +print(pop) diff --git a/bin/popjoin b/bin/popjoin index d40cfa3..0269a8c 100755 --- a/bin/popjoin +++ b/bin/popjoin @@ -3,7 +3,7 @@ import sys import argparse -import cPickle +import pickle from population import Population @@ -15,7 +15,7 @@ def PopJoin(poplist, outfile): # create "main" population object, using first in list if len(poplist)<2: - print "PopJoin requires at least two populations!" + print("PopJoin requires at least two populations!") sys.exit() # ok, we have multiple file names. @@ -23,9 +23,9 @@ def PopJoin(poplist, outfile): try: f = open(poplist[0], 'rb') except IOError: - print "Could not open file {0}.".format(poplist[0]) + print("Could not open file {0}.".format(poplist[0])) sys.exit() - pop1 = cPickle.load(f) + pop1 = pickle.load(f) f.close() # try to read in the other populations @@ -34,10 +34,10 @@ def PopJoin(poplist, outfile): try: f = open(popfile, 'rb') except IOError: - print "Could not open file {0}.".format(popfile) + print("Could not open file {0}.".format(popfile)) sys.exit() - pop = cPickle.load(f) + pop = pickle.load(f) f.close() populations.append(pop) @@ -46,11 +46,11 @@ def PopJoin(poplist, outfile): newpop = pop1.join(populations) # print out useful info - print newpop + print(newpop) # write out the new file! f = open(outfile, 'wb') - cPickle.dump(newpop, f) + pickle.dump(newpop, f) f.close() diff --git a/bin/wxHist b/bin/wxHist index b575a1a..3d3b94a 100755 --- a/bin/wxHist +++ b/bin/wxHist @@ -6,7 +6,7 @@ import argparse import random import math -import cPickle +import pickle import numpy as np import wx diff --git a/bin/wxView b/bin/wxView index 9e2ac5e..44ca5c4 100755 --- a/bin/wxView +++ b/bin/wxView @@ -5,7 +5,7 @@ import os import argparse import random -import cPickle +import pickle import numpy as np import wx @@ -133,8 +133,8 @@ class VisualizeFrame(wx.Frame): def onRadioClick(self, event): radioBox = event.GetEventObject() - print event.GetId() - print radioBox + print(event.GetId()) + print(radioBox) def onXRadioClick(self, event): # get x index of selected button diff --git a/examples/populate_and_survey.py b/examples/populate_and_survey.py index 2575a98..6df6265 100644 --- a/examples/populate_and_survey.py +++ b/examples/populate_and_survey.py @@ -3,13 +3,14 @@ # import the populate and dosurvey modules import populate import dosurvey +import matplotlib.pyplot as plt -# run the Populate.generate script, setting parameters as necessary. +# run the populate.generate script, setting parameters as necessary. # the resulting model is stored in the variable 'pop' # Any unspecified variables use the default values (see populate # for more details) -pop = populate.generate(1038, +pop = populate.generate(0.1, surveyList=['PMSURV'], radialDistType='lfl06', siDistPars=[-1.41, 0.96], # non-standard SI distribution @@ -25,4 +26,13 @@ dosurvey.write(surveyPopulations, nores=False, summary=False, asc=False) # write out the population model, if required -pop.write(outf="populate.model") +# pop.write(outf="populate.model") + +# lists to store p/pdot +periods = [pulsar.period for pulsar in pop.population if not pulsar.dead] +pdots = [pulsar.pdot for pulsar in pop.population if not pulsar.dead] + +# plot a scatter log-log plot of the p/pdot values +plt.loglog(periods, pdots, 'C0.') +plt.show() +# plt.savefig("test1.png") diff --git a/examples/ppdot.py b/examples/ppdot.py index 8f48e7a..705ccf6 100755 --- a/examples/ppdot.py +++ b/examples/ppdot.py @@ -1,14 +1,23 @@ #!/usr/bin/python +""" +********************************************** +python3 PsrPopPy/examples/ppdot.py +********************************************** +1 argument --> file_name --> *.model file +""" + import sys -import cPickle +import pickle +import matplotlib +matplotlib.use('TkAgg') import matplotlib.pyplot as plt # open file, read in model filename = sys.argv[1] f = open(filename, 'rb') -pop = cPickle.load(f) +pop = pickle.load(f) f.close() # lists to store p/pdot @@ -17,4 +26,5 @@ # plot a scatter log-log plot of the p/pdot values plt.loglog(periods, pdots, 'k.') -plt.show() \ No newline at end of file +plt.plot("test8.png") +plt.show() diff --git a/img/evolve_e5.png b/img/evolve_e5.png new file mode 100644 index 0000000..5ed1a63 Binary files /dev/null and b/img/evolve_e5.png differ diff --git a/img/old_gui.png b/img/old_gui.png new file mode 100644 index 0000000..5f817e3 Binary files /dev/null and b/img/old_gui.png differ diff --git a/img/populate_e5.png b/img/populate_e5.png new file mode 100644 index 0000000..92a0b59 Binary files /dev/null and b/img/populate_e5.png differ diff --git a/img/tkview_hist_gl.png b/img/tkview_hist_gl.png new file mode 100644 index 0000000..916a06f Binary files /dev/null and b/img/tkview_hist_gl.png differ diff --git a/img/tkview_hist_r0.png b/img/tkview_hist_r0.png new file mode 100644 index 0000000..54910ce Binary files /dev/null and b/img/tkview_hist_r0.png differ diff --git a/img/tkview_plot_dm_p.png b/img/tkview_plot_dm_p.png new file mode 100644 index 0000000..526cf42 Binary files /dev/null and b/img/tkview_plot_dm_p.png differ diff --git a/img/visualize_hist_populate_d.png b/img/visualize_hist_populate_d.png new file mode 100644 index 0000000..9e63d5b Binary files /dev/null and b/img/visualize_hist_populate_d.png differ diff --git a/img/visualize_hist_populate_si.png b/img/visualize_hist_populate_si.png new file mode 100644 index 0000000..a359b4d Binary files /dev/null and b/img/visualize_hist_populate_si.png differ diff --git a/img/visualize_plot_evolve_ppdot.png b/img/visualize_plot_evolve_ppdot.png new file mode 100644 index 0000000..0ef85f4 Binary files /dev/null and b/img/visualize_plot_evolve_ppdot.png differ diff --git a/psrpoppy/dataobj.py b/psrpoppy/dataobj.py index 1ba8d70..52d2ff8 100644 --- a/psrpoppy/dataobj.py +++ b/psrpoppy/dataobj.py @@ -2,7 +2,7 @@ import os -import cPickle +import pickle class DataObj: @@ -50,8 +50,8 @@ def readfile_return_dataobj(filename): return try: - pop = cPickle.load(f) - except cPickle.UnpicklingError: + pop = pickle.load(f) + except pickle.UnpicklingError: print "File {0} could not be unpickled!".format(filename) return f.close() diff --git a/psrpoppy/dosurvey.py b/psrpoppy/dosurvey.py index 86ca0bb..8780a6d 100755 --- a/psrpoppy/dosurvey.py +++ b/psrpoppy/dosurvey.py @@ -5,7 +5,7 @@ import math import random -import cPickle +import pickle try: # try and import from the current path (for package usage or use as an uninstalled executable) @@ -38,11 +38,11 @@ def __init__(self, def loadModel(popfile='populate.model', popmodel=None): - """Loads in either a model from disk (popfile, cPickle), + """Loads in either a model from disk (popfile, pickle), or pass in a model from memory (popmodel)""" if popmodel is None: with open(popfile, 'rb') as f: - pop = cPickle.load(f) + pop = pickle.load(f) else: pop = popmodel @@ -97,8 +97,8 @@ def run(pop, # print the population if not nostdout: - print "Running doSurvey on population..." - print pop + print("Running doSurvey on population...") + print(pop) # loop over the surveys we want to run on the pop file surveyPops = [] @@ -106,7 +106,7 @@ def run(pop, s = Survey(surv) s.discoveries = 0 if not nostdout: - print "\nRunning survey {0}".format(surv) + print("\nRunning survey {0}".format(surv)) # create a new population object to store discovered pulsars in survpop = Population() @@ -154,13 +154,13 @@ def run(pop, # report the results if not nostdout: - print "Total pulsars in model = {0}".format(len(pop.population)) - print "Number detected by survey {0} = {1}".format(surv, ndet) - print "Of which are discoveries = {0}".format(s.discoveries) - print "Number too faint = {0}".format(ntf) - print "Number smeared = {0}".format(nsmear) - print "Number out = {0}".format(nout) - print "\n" + print("Total pulsars in model = {0}".format(len(pop.population))) + print("Number detected by survey {0} = {1}".format(surv, ndet)) + print("Of which are discoveries = {0}".format(s.discoveries)) + print("Number too faint = {0}".format(ntf)) + print("Number smeared = {0}".format(nsmear)) + print("Number out = {0}".format(nout)) + print("\n") d = Detections(ndet=ndet, ntf=ntf, @@ -185,7 +185,7 @@ def run(pop, description='Run a survey on your population model') parser.add_argument( '-f', metavar='fname', default='populate.model', - help='file containing population model (def=populate.model') + help='file containing population model (def=populate.model)') parser.add_argument( '-surveys', metavar='S', nargs='+', required=True, diff --git a/psrpoppy/evolve.py b/psrpoppy/evolve.py index 8963747..1de3967 100755 --- a/psrpoppy/evolve.py +++ b/psrpoppy/evolve.py @@ -6,10 +6,14 @@ import random import inspect -import cPickle +import pickle import scipy.integrate import numpy as np +import matplotlib +matplotlib.use('TkAgg') +import matplotlib.pyplot as plt + try: # try and import from the current path (for package usage or use as an uninstalled executable) import distributions as dists @@ -109,30 +113,22 @@ def generate(ngen, pop.zscale = zscale if widthModel == 'kj07': - print "\tLoading KJ07 models...." + print("\tLoading KJ07 models....") kj_p_vals, kj_pdot_vals, kj_dists = beammodels.load_kj2007_models() - print "\tDone\n" + print("\tDone\n") if not nostdout: - print "\tGenerating evolved pulsars with parameters:" - print "\t\tngen = {0}".format(ngen) - print "\t\tUsing electron distn model {0}".format( - pop.electronModel) - print "\n\t\tPeriod mean, sigma = {0}, {1}".format( - pop.pmean, - pop.psigma) - print "\t\tLuminosity mean, sigma = {0}, {1}".format( - pop.lummean, - pop.lumsigma) - print "\t\tSpectral index mean, sigma = {0}, {1}".format( - pop.simean, - pop.sisigma) - print "\t\tGalactic z scale height = {0} kpc".format( - pop.zscale) + print("\tGenerating evolved pulsars with parameters:") + print("\t\tngen = {0}".format(ngen)) + print("\t\tUsing electron distn model {0}".format(pop.electronModel)) + print("\n\t\tPeriod mean, sigma = {0}, {1}".format(pop.pmean, pop.psigma)) + print("\t\tLuminosity mean, sigma = {0}, {1}".format(pop.lummean, pop.lumsigma)) + print("\t\tSpectral index mean, sigma = {0}, {1}".format(pop.simean, pop.sisigma)) + print("\t\tGalactic z scale height = {0} kpc".format(pop.zscale)) if widthModel is None: - print "\t\tWidth {0}% ".format(duty) + print("\t\tWidth {0}% ".format(duty)) else: - print "\t\tUsing Karastergiou & Johnston beam width model" + print("\t\tUsing Karastergiou & Johnston beam width model") # set up progress bar for fun :) prog = ProgressBar(min_value=0, @@ -233,7 +229,7 @@ def generate(ngen, """ else: - print "Undefined width model!" + print("Undefined width model!") sys.exit() # print width # print pulsar.period, width, pulsar.pdot @@ -287,8 +283,9 @@ def generate(ngen, pulsar.spindex = random.gauss(pop.simean, pop.sisigma) # calculate galactic coords and distance - pulsar.gl, pulsar.gb = go.xyz_to_lb(pulsar.galCoords) - pulsar.dtrue = go.calc_dtrue(pulsar.galCoords) + x, y, z = pulsar.galCoords + pulsar.gl, pulsar.gb = go.xyz_to_lb(x, y, z) + pulsar.dtrue = go.calc_dtrue(x, y, z) # then calc DM using fortran libs if pop.electronModel == 'ne2001': @@ -343,7 +340,7 @@ def generate(ngen, # update the counter if not nostdout: prog.increment_amount() - print prog, '\r', + print(prog, '\r',) sys.stdout.flush() else: @@ -353,7 +350,7 @@ def generate(ngen, # update the counter if not nostdout: prog.increment_amount() - print prog, '\r', + print(prog, '\r',) sys.stdout.flush() # pulsar isn't dead, add to population! @@ -367,25 +364,26 @@ def generate(ngen, # update the counter if not nostdout: prog.increment_amount() - print prog, '\r', + print(prog, '\r',) sys.stdout.flush() if not nostdout: - print "\n\n" - print " Total pulsars = {0}".format(len(pop.population)) - print " Total detected = {0}".format(pop.ndet) + print("\n\n") + print(" Total pulsars = {0}".format(len(pop.population))) + print(" Total detected = {0}".format(pop.ndet)) for surv in surveys: - print "\n Results for survey '{0}'".format(surv.surveyName) - print " Number detected = {0}".format(surv.ndet) - print " Number too faint = {0}".format(surv.ntf) - print " Number smeared = {0}".format(surv.nsmear) - print " Number outside survey area = {0}".format(surv.nout) + print("\n Results for survey '{0}'".format(surv.surveyName)) + print(" Number detected = {0}".format(surv.ndet)) + print(" Number too faint = {0}".format(surv.ntf)) + print(" Number smeared = {0}".format(surv.nsmear)) + print(" Number outside survey area = {0}".format(surv.nout)) # save list of arguments into the pop try: - argspec = inspect.getargspec(generate) - key_values = [(arg, locals()[arg]) for arg in argspec.args] + argspec = inspect.getfullargspec(generate) + l = locals() + key_values = [(arg, l[arg]) for arg in argspec.args] pop.arguments = {key: value for (key, value) in key_values} except SyntaxError: pass @@ -435,7 +433,7 @@ def alignpulsar(pulsar, pop): """ # need to be careful to use chi in degrees, but - # consistently remember to take the sins/cosines of radians + # consistently remember to take the sines/cosines of radians if pop.alignModel == 'orthogonal': pulsar.chi = 90.0 pulsar.sinchi_init = 1.0 @@ -770,7 +768,7 @@ def bhattacharya_deathperiod_92(pulsar): help='Output filename for population model \ (def=evolve.model)') - # turn off printing to stdout + # turn off print(ng to stdout) parser.add_argument('--nostdout', nargs='?', const=True, default=False, help='switch off std output') @@ -820,3 +818,12 @@ def bhattacharya_deathperiod_92(pulsar): keepdead=args.keepdead) pop.write(outf=args.o) + # print([i for i in pop.population]) + periods = [pulsar.period for pulsar in pop.population] + pdots = [pulsar.pdot for pulsar in pop.population] + + # plot a scatter log-log plot of the p/pdot values + plt.loglog(periods, pdots, 'C0.') + plt.xlabel("log P") + plt.ylabel(r"log $\dot{P}$") + plt.show() \ No newline at end of file diff --git a/psrpoppy/galacticops.py b/psrpoppy/galacticops.py index b5d775f..9b90108 100755 --- a/psrpoppy/galacticops.py +++ b/psrpoppy/galacticops.py @@ -62,7 +62,7 @@ def vxyz(pulsar): pulsar.vz = vz.value -def calc_dtrue((x, y, z)): +def calc_dtrue(x, y, z): """Calculate true distance to pulsar from the sun.""" rsun = 8.5 # kpc return math.sqrt(x*x + (y-rsun)*(y-rsun) + z*z) @@ -86,7 +86,7 @@ def ne2001_dist_to_dm(dist, gl, gb): dist = C.c_float(dist) gl = C.c_float(gl) gb = C.c_float(gb) - inpath = C.create_string_buffer(fortranpath) + inpath = C.create_string_buffer(str.encode(fortranpath)) linpath = C.c_int(len(fortranpath)) return ne2001lib.dm_(C.byref(dist), C.byref(gl), @@ -105,7 +105,7 @@ def lmt85_dist_to_dm(dist, gl, gb): # passing path to fortran dir and the length of # this path --- removes need to edit getpath.f # during installation - inpath = C.create_string_buffer(fortranpath) + inpath = C.create_string_buffer(str.encode(fortranpath)) linpath = C.c_int(len(fortranpath)) return ne2001lib.dm_(C.byref(dist), @@ -198,7 +198,7 @@ def radec_to_lb(ra, dec): return l.value, b.value -def xyz_to_lb((x, y, z)): +def xyz_to_lb(x, y, z): """ Convert galactic xyz in kpc to l and b in degrees.""" rsun = 8.5 # kpc diff --git a/psrpoppy/orbitalparams.py b/psrpoppy/orbitalparams.py index 70210b1..13952b9 100644 --- a/psrpoppy/orbitalparams.py +++ b/psrpoppy/orbitalparams.py @@ -17,10 +17,10 @@ def test_1802_2124(pulsar): pulsar.inclination_degrees = 78.52 pulsar.pulsar_mass_msolar = 1.24 - print pulsar.gb, pulsar.gl + print(pulsar.gb, pulsar.gl) pulsar.gb = 0.61 pulsar.gl = 4.38 - print pulsar.gb, pulsar.gl + print(pulsar.gb, pulsar.gl) pulsar.galcoords = (0.49, 5.2, 0.04) pulsar.lum_1400 = 8.54 pulsar.t_scatter = 0.0 diff --git a/psrpoppy/populate.py b/psrpoppy/populate.py index 8f4285d..db93713 100755 --- a/psrpoppy/populate.py +++ b/psrpoppy/populate.py @@ -6,7 +6,11 @@ import random import inspect -import cPickle +import pickle + +import matplotlib +matplotlib.use('TkAgg') +import matplotlib.pyplot as plt try: # try and import from the current path (for package usage or use as an uninstalled executable) @@ -39,7 +43,7 @@ class PopulateException(Exception): pass -def generate(ngen, +def generate(ngen=1038, surveyList=None, pDistType='lnorm', radialDistType='lfl06', @@ -84,23 +88,23 @@ def generate(ngen, # check that the distribution types are supported.... if lumDistType not in ['lnorm', 'pow']: - print "Unsupported luminosity distribution: {0}".format(lumDistType) + print("Unsupported luminosity distribution: {0}".format(lumDistType)) if pDistType not in ['lnorm', 'norm', 'cc97', 'lorimer12']: - print "Unsupported period distribution: {0}".format(pDistType) + print("Unsupported period distribution: {0}".format(pDistType)) if radialDistType not in ['lfl06', 'yk04', 'isotropic', 'slab', 'disk', 'gauss']: - print "Unsupported radial distribution: {0}".format(radialDistType) + print("Unsupported radial distribution: {0}".format(radialDistType)) if electronModel not in ['ne2001', 'lmt85']: - print "Unsupported electron model: {0}".format(electronModel) + print("Unsupported electron model: {0}".format(electronModel)) if pattern not in ['gaussian', 'airy']: - print "Unsupported gain pattern: {0}".format(pattern) + print("Unsupported gain pattern: {0}".format(pattern)) if duty_percent < 0.: - print "Unsupported value of duty cycle: {0}".format(duty_percent) + print("Unsupported value of duty cycle: {0}".format(duty_percent)) # need to use properties in this class so they're get/set-type props pop.pDistType = pDistType @@ -130,14 +134,15 @@ def generate(ngen, # store the dict of arguments inside the model. Could be useful. try: - argspec = inspect.getargspec(generate) - key_values = [(arg, locals()[arg]) for arg in argspec.args] + argspec = inspect.getfullargspec(generate) + l = locals() + key_values = [(arg, l[arg]) for arg in argspec.args] pop.arguments = {key: value for (key, value) in key_values} except SyntaxError: pass if not nostdout: - print "\tGenerating pulsars with parameters:" + print("\tGenerating pulsars with parameters:") param_string_list = [] for key, value in key_values: s = ": ".join([key, str(value)]) @@ -145,10 +150,10 @@ def generate(ngen, # join this list of strings, and print it s = "\n\t\t".join(param_string_list) - print "\t\t{0}".format(s) + print("\t\t{0}".format(s)) # set up progress bar for fun :) - prog = ProgressBar(min_value=0, + prog = ProgressBar(min_value=0, max_value=ngen, width=65, mode='dynamic') @@ -182,11 +187,13 @@ def generate(ngen, elif pop.pDistType == 'cc97': p.period = _cc97() elif pop.pDistType == 'gamma': - print "Gamma function not yet supported" + print("Gamma function not yet supported") sys.exit() elif pop.pDistType == 'lorimer12': p.period = _lorimer2012_msp_periods() + p.pdot = 9.768E-40 * random.gauss(12.65, 0.55)**2 * (p.period/1000)**(2.0 - 2.5 + 0.5 * random.random()) + if duty_percent > 0.: # use a simple duty cycle for each pulsar # with a log-normal scatter @@ -213,7 +220,7 @@ def generate(ngen, # suppose it might be nice to be able to have GPS sources # AND double spectra. But for now I assume only have one or # none of these types. - if random.random() > pop.gpsFrac: + if pop.gpsFrac and random.random() > pop.gpsFrac: # This will evaluate true when gpsArgs[0] is NoneType # might have to change in future p.gpsFlag = 0 @@ -221,7 +228,7 @@ def generate(ngen, p.gpsFlag = 1 p.gpsA = pop.gpsA - if random.random() > pop.brokenFrac: + if pop.brokenFrac and random.random() > pop.brokenFrac: p.brokenFlag = 0 else: p.brokenFlag = 1 @@ -268,9 +275,9 @@ def generate(ngen, zheight = random.gauss(0., zscale) gx, gy = go.calcXY(p.r0) p.galCoords = gx, gy, zheight - p.gl, p.gb = go.xyz_to_lb(p.galCoords) + p.gl, p.gb = go.xyz_to_lb(gx, gy, zheight) - p.dtrue = go.calc_dtrue(p.galCoords) + p.dtrue = go.calc_dtrue(gx, gy, zheight) # then calc DM using fortran libs if pop.electronModel == 'ne2001': @@ -293,7 +300,7 @@ def generate(ngen, # add in orbital parameters if orbits: orbitalparams.test_1802_2124(p) - print p.gb, p.gl + print(p.gb, p.gl) # if no surveys, just generate ngen pulsars if surveyList is None: @@ -301,7 +308,7 @@ def generate(ngen, pop.ndet += 1 if not nostdout: prog.increment_amount() - print prog, '\r', + print(prog, '\r',) sys.stdout.flush() # if surveys are given, check if pulsar detected or not # in ANY of the surveys @@ -345,22 +352,22 @@ def generate(ngen, pop.ndet += 1 if not nostdout: prog.increment_amount() - print prog, '\r', + print(prog, '\r',) sys.stdout.flush() # print info to stdout if not nostdout: - print "\n" - print " Total pulsars = {0}".format(len(pop.population)) - print " Total detected = {0}".format(pop.ndet) + print("\n") + print(" Total pulsars = {0}".format(len(pop.population))) + print(" Total detected = {0}".format(pop.ndet)) # print " Number not beaming = {0}".format(surv.nnb) for surv in surveys: - print "\n Results for survey '{0}'".format(surv.surveyName) - print " Number detected = {0}".format(surv.ndet) - print " Number too faint = {0}".format(surv.ntf) - print " Number smeared = {0}".format(surv.nsmear) - print " Number outside survey area = {0}".format(surv.nout) + print("\n Results for survey '{0}'".format(surv.surveyName)) + print(" Number detected = {0}".format(surv.ndet)) + print(" Number too faint = {0}".format(surv.ntf)) + print(" Number smeared = {0}".format(surv.nsmear)) + print(" Number outside survey area = {0}".format(surv.nout)) return pop @@ -376,13 +383,12 @@ def _lorimer2012_msp_periods(): dist = [1., 3., 5., 16., 9., 5., 5., 3., 2.] # calculate which bin to take value of - # bin_num = dists.draw1d(dist) # assume linear distn inside the bins - logp = logpmin + (logpmax-logpmin)*(bin_num+random.random())/len(dist) + logp = logpmin + (logpmax - logpmin) * (bin_num + random.random()) / len(dist) - return 10.**logp + return 10 ** logp def _cc97(): @@ -561,9 +567,9 @@ def _sindegree(angle): f.write(' '.join(sys.argv)) f.write('\n') - # run the code and write out a cPickle population class + # run the code and write out a Pickle population class - pop = generate(args.n, + pop = generate(ngen=args.n, surveyList=args.surveys, pDistType=args.pdist[0], lumDistType=args.ldist[0], @@ -584,3 +590,14 @@ def _sindegree(angle): ) pop.write(outf=args.o) + + # print([(i.period, i.pdot) for i in pop.population]) + # lists to store p/pdot + periods = [pulsar.period for pulsar in pop.population] + pdots = [pulsar.pdot for pulsar in pop.population] + + # plot a scatter log-log plot of the p/pdot values + plt.loglog(periods, pdots, "C0.") + plt.xlabel("log P") + plt.ylabel(r"log $\dot{P}$") + plt.show() diff --git a/psrpoppy/population.py b/psrpoppy/population.py index ef28306..283877e 100755 --- a/psrpoppy/population.py +++ b/psrpoppy/population.py @@ -1,7 +1,7 @@ #!/usr/bin/python import copy -import cPickle +import pickle import numpy as np @@ -102,7 +102,7 @@ def write(self, outf): """Write the population object to a file""" output = open(outf, 'wb') - cPickle.dump(self, output, 2) + pickle.dump(self, output, 2) output.close() def write_asc(self, outf): diff --git a/psrpoppy/progressbar.py b/psrpoppy/progressbar.py index c5dc54c..628d5a9 100644 --- a/psrpoppy/progressbar.py +++ b/psrpoppy/progressbar.py @@ -98,26 +98,26 @@ def main(): print limit = 1000000 - print 'Example 1: Fixed Bar' + print('Example 1: Fixed Bar') prog = ProgressBar(0, limit, 77, mode='fixed') oldprog = str(prog) for i in xrange(limit+1): prog.update_amount(i) if oldprog != str(prog): - print prog, "\r", + print(prog, "\r",) sys.stdout.flush() oldprog = str(prog) - print '\n\n' + print('\n\n') - print 'Example 2: Dynamic Bar' + print('Example 2: Dynamic Bar') prog = ProgressBar(0, limit, 77, mode='dynamic', char='-') oldprog = str(prog) for i in xrange(limit+1): prog.increment_amount() if oldprog != str(prog): - print prog, "\r", + print(prog, "\r",) sys.stdout.flush() oldprog = str(prog) - print '\n\n' + print('\n\n') diff --git a/psrpoppy/pulsar.py b/psrpoppy/pulsar.py index f40cd48..000c438 100644 --- a/psrpoppy/pulsar.py +++ b/psrpoppy/pulsar.py @@ -47,7 +47,7 @@ def __init__(self, self.dm = dm # convert to -180->+180 range - if gl > 180.: + if gl and gl > 180.: gl -= 360. self.gl = gl diff --git a/psrpoppy/survey.py b/psrpoppy/survey.py index ab09d73..aad5780 100755 --- a/psrpoppy/survey.py +++ b/psrpoppy/survey.py @@ -274,7 +274,7 @@ def __init__(self, surveyName, pattern='gaussian'): # turn on AA self.AA = True else: - print "Parameter '", a[1].strip(), "' not recognized!" + print("Parameter '", a[1].strip(), "' not recognized!") f.close() @@ -455,37 +455,37 @@ def SNRcalc(self, if pulsar.is_binary: # print "the pulsar is a binary!" if jerksearch: - print "jerk" + print("jerk") gamma = degradation.gamma3(pulsar, self.tobs, 1) elif accelsearch: - print "accel" + print("accel") gamma = degradation.gamma2(pulsar, self.tobs, 1) else: - print "norm" + print("norm") gamma = degradation.gamma1(pulsar, self.tobs, 1) - print "gamma harm1 = ", gamma + print("gamma harm1 = ", gamma) gamma = degradation.gamma1(pulsar, self.tobs, 2) - print "gamma harm2 = ", gamma + print("gamma harm2 = ", gamma) gamma = degradation.gamma1(pulsar, self.tobs, 3) - print "gamma harm3 = ", gamma + print("gamma harm3 = ", gamma) gamma = degradation.gamma1(pulsar, self.tobs, 4) - print "gamma harm4 = ", gamma + print ("gamma harm4 = ", gamma) # return the S/N accounting for beam offset return sig_to_noise * degfac @@ -518,9 +518,9 @@ def _gpsFlux(self, psr, ref_freq): """Calculate the flux assuming GPS spectrum shape, spindex===b""" log_nu_1 = math.log10(ref_freq/1000.) log_nu_2 = math.log10(self.freq/1000.) - gpsC = math.log10(psr.s_1400()) - (psr.gpsA * log_nu_1**2) \ + gpsC = math.log10(psr.s_1400()) - (psr.gpsA if psr.gpsA else 0 * log_nu_1**2) \ - psr.spindex * log_nu_1 - return 10.**(psr.gpsA * log_nu_2**2 + psr.spindex * log_nu_2 + gpsC) + return 10.**(psr.gpsA if psr.gpsA else 0 * log_nu_2**2 + psr.spindex * log_nu_2 + gpsC) def _dmsmear(self, psr): """Calculate the smearing across a channel due to the pulsar DM""" diff --git a/tests/testPopulate.py b/tests/testPopulate.py index aea93bc..32b5fed 100644 --- a/tests/testPopulate.py +++ b/tests/testPopulate.py @@ -11,11 +11,11 @@ sys.path.append(parentdir) # import the module to be tested -from populate import Populate +from populate import generate class testPopulate(unittest.TestCase): def setUp(self): - self.populate = Populate() + self.populate = generate(1038) def test_double_sided_exp_zero(self): result = self.populate. _double_sided_exp(0.0, origin=0.0) @@ -29,5 +29,4 @@ def test_double_sided_exp_scale(self): # this one is more tricky since there are randoms involved # simple test - check in right range result = self.populate. _double_sided_exp(10.0, origin=0.0) - self.assertTrue(math.fabs(result) < - and math.fabs(result) > ) + self.assertTrue(math.fabs(result) < 10**9 and math.fabs(result) > 0)