In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from numpy import poly1d, polyfit, power
import scipy.optimize
from math import *
from IPython.display import HTML
from IPython.display import Image
import os
import pandas as pd
import PIL as pil
import heapq
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 18, 14
# import seaborn as sns
# sns.set_palette("deep", desat=.6)
try:
from PIL import Image
except:
import Image
incl = 40.
sin_i2 = np.sin(incl*np.pi/180.)**2
cos_i2 = np.cos(incl*np.pi/180.)**2
In [2]:
os.chdir("C:\\science\\2FInstability\\data\\ngc2775")
plt.imshow(np.asarray(Image.open("shapiro_fg3c.png")))
plt.plot()
Out[2]:
In [3]:
plt.imshow(np.asarray(Image.open("shapiro_fg4c.png")))
plt.plot()
Out[3]:
In [4]:
pylab.rcParams['figure.figsize'] = 15, 5
#Velocity
lu = (55, 15)
lu_val = (26, 220)
rd = (417, 135)
rd_val = (64, 80)
data = [(105,71), (104,80), (133,70), (138,67), (157,67), (170,67), (177,64), (203,65),
(235,60), (252,54), (285,57), (318,49), (350,60), (383,52), (399,67)]
def extend_values(data, lu, lu_val, rd, rd_val):
xscale = 1.0*(rd_val[0]-lu_val[0])/(rd[0]-lu[0])
yscale = 1.0*(rd_val[1]-lu_val[1])/(lu[1]-rd[1])
extended = []
for d in data:
extended.append((xscale*(d[0]-lu[0])+lu_val[0], yscale*(rd[1] - d[1])+rd_val[1]))
return extended
e_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*e_data)[0],zip(*e_data)[1], 'or')
vel_fit = lambda l: 283.*np.power(l, 0)
plt.plot(np.linspace(30., 90., 100), map(vel_fit, np.linspace(30., 90., 100)), '--')
plt.ylim(80, 220)
plt.xlim(26, 64)
plt.show()
In [5]:
#Sig_maj
lu = (56, 136)
lu_val = (26, 110)
rd = (417, 256)
rd_val = (64, 30)
data = [(105,153), (105,163), (134,184), (137,177), (157,189),
(169,194), (176,189), (203,198), (203,201), (236,205), (252,198), (285,203), (318,200), (350,197), (383,204), (400,203)]
sig_maj_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*sig_maj_data)[0],zip(*sig_maj_data)[1], 'or')
sigR = lambda l: 197.*np.exp(-l/45.)
sigZ = lambda l: 201.*np.exp(-l/45.)
phi_to_R = lambda l: 0.5*(1 - 0.)
sig_maj = lambda l: sqrt(sigR(l)**2 * (phi_to_R(l) * sin_i2 + sigZ(l)**2 * cos_i2/sigR(l)**2))
plt.plot(np.linspace(30., 90., 100), map(sig_maj, np.linspace(30., 90., 100)), '--')
plt.ylim(30, 110)
plt.xlim(26, 64)
plt.show()
In [6]:
#Sig_min
lu = (56, 256)
lu_val = (26, 110)
rd = (417, 376)
rd_val = (64, 30)
data = [(117,284), (143,281), (160,309), (185,306), (199,315), (228,325), (228,302), (254,315), (271,313), (288,313),
(314,329), (353,332), (357,322), (400,347)]
sig_min_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*sig_min_data)[0], zip(*sig_min_data)[1], 'or')
sig_min = lambda l: sqrt(sigR(l)**2 * sin_i2 + sigZ(l)**2 * cos_i2)
plt.plot(np.linspace(30., 90., 100), map(sig_min, np.linspace(30., 90., 100)), '--')
plt.ylim(30, 110)
plt.xlim(26, 64)
plt.show()
In [7]:
h_kin = 45.
beta = 0.
def sig_maj_exp(R):
global alpha, sig_R_0, h_kin, beta
return np.exp(-R/h_kin)*sig_R_0*sqrt(0.5*(1 - beta) * sin_i2 + alpha**2 * cos_i2)
def sig_min_exp(R):
global alpha, sig_R_0, h_kin
return np.exp(-R/h_kin)*sig_R_0*sqrt(sin_i2 + alpha**2 * cos_i2)
def calc_chi2_normal(obs, obserr, predicted):
return sum([(o-p)**2/err**2 for (o,p,err) in zip(obs, predicted, obserr)])/(len(obs))
alphas = np.arange(0.25, 1., 0.01)
sigmas = np.arange(80.0, 400, 0.25)
def compute_chi2_maps(alph=(), sigm=()):
'''Вычисляем все изображения, чтобы потом только настройки менять'''
image_min = np.random.uniform(size=(len(sigm), len(alph)))
image_maj = np.random.uniform(size=(len(sigm), len(alph)))
image = np.random.uniform(size=(len(sigmas), len(alphas)))
for i,si in enumerate(sigm):
for j,al in enumerate(alph):
global alpha, sig_R_0
alpha = al
sig_R_0 = si
sqerr_maj = sum(power([sig_maj_exp(p[0]) - p[1] for p in sig_maj_data], 2))/len(sig_maj_data)
sqerr_min = sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data)
image_maj[i][j] = sqerr_maj
image_min[i][j] = sqerr_min
return image_maj, image_min
pics_path = "C:\\science\\2FInstability\\data\\ngc2775\\"
if not os.path.exists(pics_path):
os.makedirs(pics_path)
if os.path.isfile(pics_path + 'chi2_map_maj.npy'):
image_maj = np.load(pics_path + "chi2_map_maj.npy")
image_min = np.load(pics_path + "chi2_map_min.npy")
else:
image_maj, image_min = compute_chi2_maps(alph=alphas, sigm=sigmas)
np.save(pics_path + 'chi2_map_maj', image_maj)
np.save(pics_path + 'chi2_map_min', image_min)
In [8]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
def plot_chi2_map(image, ax, log_scale=False, title='$\chi^2$', is_contour=False, vmax=0.):
if image is not None:
im = ax.imshow(image, cmap='jet', vmin=image.min(), vmax=vmax, interpolation='spline16',
origin="lower", extent=[alphas[0], alphas[-1],sigmas[0],sigmas[-1]], aspect="auto")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(im, cax=cax)
ax.set_title(title, size=20.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
ax.set_xlabel(r'$\alpha$', size=20.)
ax.grid(True)
min_sigmas = np.where(image_min < image_min.min() + 10.)
slice_alph, slice_sig = min_sigmas[1], min_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_min[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_min), slice_sig)
poly_slice_min = poly1d(polyfit(slice_alph, slice_sig, deg=3))
maj_sigmas = np.where(image_maj < image_maj.min() + 10.)
slice_alph, slice_sig = maj_sigmas[1], maj_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_maj[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_maj), slice_sig)
poly_slice_maj = poly1d(polyfit(slice_alph, slice_sig, deg=3))
fig, axes = plt.subplots(nrows=1, ncols=2, sharex=False, sharey=True, figsize=[15,5])
plot_chi2_map(image_maj, axes[0], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=300.)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=300.)
axes[0].plot(0.58, 213., 'o', color='m')
axes[0].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')
axes[0].plot(alphas[20:], map(poly_slice_maj, alphas[20:]), '.-', label = 'maj -> maj', color= 'red')
axes[0].plot(alphas[20:], map(poly_slice_min, alphas[20:]), '.-', label = 'min -> maj', color= 'pink')
axes[1].plot(0.58, 213., 'o', color='m')
axes[1].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')
axes[1].plot(alphas[20:], map(poly_slice_min, alphas[20:]), '.-', label = 'min -> min', color='red')
axes[1].plot(alphas[20:], map(poly_slice_maj, alphas[20:]), '.-', label = 'maj -> min', color='pink')
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=2, sharex=False, sharey=False, figsize=[15,5])
plot_chi2_map((image_min + image_maj)/2, axes[0], log_scale=False, title='$\chi^2_{maj}+\chi^2_{min}$', is_contour=False, vmax=150.)
axes[0].plot(0.58, 213., 'o', color='m')
axes[0].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')
err_maj_1, err_maj_2 = [], []
err_min_1, err_min_2 = [], []
for al in alphas:
global alpha, sig_R_0
alpha = al
sig_R_0 = poly_slice_min(alpha)
err_maj_1.append(sum(power([sig_maj_exp(p[0]) - p[1] for p in sig_maj_data], 2))/len(sig_maj_data))
err_min_1.append(sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data))
sig_R_0 = poly_slice_maj(alpha)
err_maj_2.append(sum(power([sig_maj_exp(p[0]) - p[1] for p in sig_maj_data], 2))/len(sig_maj_data))
err_min_2.append(sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data))
axes[1].plot(alphas, err_maj_1, '.-', label = 'min -> maj', color= 'red')
axes[1].plot(alphas, err_min_1, '.-', label = 'min -> min', color='orange')
axes[1].plot(alphas, err_maj_2, '.-', label = 'maj -> maj', color= 'green')
axes[1].plot(alphas, err_min_2, '.-', label = 'maj -> min', color='blue')
axes[1].legend()
plt.show()
In [9]:
alphas = np.arange(0.3, 0.9, 0.12)
sigmas = np.arange(210., 310., 15.)
points = np.arange(26., 64., 0.1)
good_pics = []
def plot_ranges_gers(sigmas_range, alphas_range, good_pics=[], calc_chi=False, best_err=3):
nrows = alphas.size
ncols = sigmas.size
fig, axes = plt.subplots(nrows=nrows, ncols=ncols, sharex=True, sharey=True, figsize=[16,10])
plt_index = 0
# Последнее - среднее геометрическое
sqerr_majs, sqerr_mins, sqerr_mean = [],[],[]
for al in alphas_range:
for si in sigmas_range:
global alpha, sig_R_0
alpha = al
sig_R_0 = si
ax = axes[plt_index/ncols, plt_index % ncols]
ax.set_title(r'$\alpha = %s, \sigma_{R,0}=%s$' % (al,si))
sqerr_maj = sum(power([sig_maj_exp(p[0]) - p[1] for p in sig_maj_data], 2))/len(sig_maj_data)
sqerr_min = sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data)
if sqerr_maj < 120. and sqerr_min < 120.:
ax.plot(80., 105., 'og', ms = 10)
ax.plot(points, [sig_maj_exp(R) for R in points], '--', color='blue')
ax.plot(points, [sig_min_exp(R) for R in points], '--', color='red')
ax.plot(zip(*sig_min_data)[0],zip(*sig_min_data)[1], 's', color='red', ms=3)
ax.plot(zip(*sig_maj_data)[0],zip(*sig_maj_data)[1], 's', color='blue', ms=3)
ax.set_ylim(30, 110)
ax.set_xlim(26, 64)
plt_index = plt_index + 1
plot_ranges_gers(sigmas, alphas, good_pics=good_pics, calc_chi=True)
plt.show()
In [10]:
from IPython.html.widgets import *
def widget_plot_maps(h_k, b, vm, contour_swap):
global alpha, sig_R_0, h_kin, beta, alphas, sigmas
h_kin = h_k
beta = b
alphas = np.arange(0.25, 0.9, 0.05)
sigmas = np.arange(80.0, 400, 5.)
image_maj, image_min = compute_chi2_maps(alph=alphas, sigm=sigmas)
fig, axes = plt.subplots(nrows=1, ncols=2, sharex=False, sharey=True, figsize=[15,5])
plot_chi2_map(image_maj, axes[0], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=vm)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=vm)
if contour_swap:
norm = cm.colors.Normalize(vmax=image_min.max(), vmin=-image_min.max())
cmap = cm.PRGn
levels = np.linspace(start=image_min.min(), stop=vm, num=10)
cset=axes[0].contour(image_min, levels, hold='on', colors = 'k', origin='lower',
extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
norm = cm.colors.Normalize(vmax=image_maj.max(), vmin=-image_maj.max())
cmap = cm.PRGn
levels = np.linspace(start=image_maj.min(), stop=vm, num=10)
cset=axes[1].contour(image_maj, levels, hold='on', colors = 'k', origin='lower',
extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
axes[0].plot(0.58, 213., 'o', color='m')
axes[0].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')
axes[1].plot(0.58, 213., 'o', color='m')
axes[1].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')
plt.show()
interact(widget_plot_maps, h_k=(60, 80, 0.5), b=(0.15, 0.25, 0.01), vm = (70., 350., 10), contour_swap=False);