In [54]:
import quantumpropagator as qp
import matplotlib.pyplot as plt
%matplotlib ipympl
plt.rcParams.update({'font.size': 8})
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as pltfrom
from ipywidgets import interact,fixed #, interactive, fixed, interact_manual
import ipywidgets as widgets
from matplotlib import cm
import pickle
# name_data_file = '/home/alessio/n-Propagation/newExtrapolated_allCorrection.pickle'
name_data_file = '/home/alessio/n-Propagation/newExtrapolated_gammaExtrExag.pickle'
# dataDict = np.load('/home/alessio/n-Propagation/datanewoneWithNACnow.npy')[()]
# # name_data_file = '/home/alessio/n-Propagation/NAC_2_1_little_exagerated.pickle'
with open(name_data_file, "rb") as input_file:
data = pickle.load(input_file)
%load_ext Cython
# data.keys()
In [55]:
# name_data_file2 = 'NAC_2_1_little_exagerated.pickle'
# with open(name_data_file2, "rb") as input_file:
# data2 = pickle.load(input_file)
# name_data_file3 = 'newExtrapolated_gammaExtrExag.pickle'
# with open(name_data_file3, "rb") as input_file:
# data3 = pickle.load(input_file)
In [56]:
# pot = data['potCube']
# pot2= data2['potCube']
# pot3 = data3['potCube']
# np.all(pot == pot3)
In [63]:
pot = data['potCube']
data['potCube'].shape
print(pot.shape)
qp.find_numpy_index_minumum(pot), pot[29, 28, 55, 0]
Out[63]:
In [67]:
%matplotlib ipympl
pot_difference_AU = pot[15:-15,15:-15,30:-30,2] - pot[15:-15,15:-15,30:-30,3]
#phiC, gamC, theC, phiD, gamD, theD = (27,25,85, 4, 7, 24)
#pot_difference_AU = pot[phiC-phiD:phiC+phiD,gamC-gamD:gamC+gamD,theC-theD:theC+theD,2] - pot[phiC-phiD:phiC+phiD,gamC-gamD:gamC+gamD,theC-theD:theC+theD,3]
pot_difference = qp.fromHartoEv(pot_difference_AU)
print(qp.find_numpy_index_minumum(pot_difference))
b = pd.Series(pot_difference.flatten())
b.describe()
b.hist(bins=100)
Out[67]:
In [70]:
plt.close('all')
phiC, gamC, theC, phiD, gamD, theD = (27,26,85, 2, 2, 24)
pot_difference_AU = pot[phiC-phiD:phiC+phiD,gamC-gamD:gamC+gamD,theC-theD:theC+theD,2] - pot[phiC-phiD:phiC+phiD,gamC-gamD:gamC+gamD,theC-theD:theC+theD,3]
pot_difference = qp.fromHartoEv(pot_difference_AU)
print(qp.find_numpy_index_minumum(pot_difference))
b = pd.Series(pot_difference.flatten())
b.describe()
b.hist(bins=100)
Out[70]:
In [59]:
%matplotlib ipympl
dp = 7
dg = 7
dt = 20
mask = pot_difference[22-dp:22+dp,22-dg:22+dg,110-dt:110+dt]
c = pd.Series(mask.flatten())
c.describe()
#c.hist()
Out[59]:
In [60]:
diff_0_1_all = pot[:,:,:,1]-pot[:,:,:,0]
diff_2_3_all = pot[:,:,:,3]-pot[:,:,:,2]
In [62]:
diff_0_1 = np.zeros_like(diff_0_1_all) + 999
diff_0_1[15:-15,15:-15,30:-30] = diff_0_1_all[15:-15,15:-15,30:-30]
diff_2_3 = np.zeros_like(diff_2_3_all) + 999
diff_2_3[15:-15,15:-15,30:-30] = diff_2_3_all[15:-15,15:-15,30:-30]
save_pot_diff = True
dictio = {}
a = 0
if save_pot_diff:
filename = '/home/alessio/IMPORTANTS/VISUALIZE_ENERGY_DIFFERENCE/PotDiff{:04}.h5'.format(a)
# I reput the zeros out.
dictio['diff'] = qp.fromHartoEv(diff_0_1)
dictio['lab'] = 'Diff 0 1'
qp.writeH5fileDict(filename, dictio)
a = a + 1
if save_pot_diff:
filename = '/home/alessio/IMPORTANTS/VISUALIZE_ENERGY_DIFFERENCE/PotDiff{:04}.h5'.format(a)
# I reput the zeros out.
dictio['diff'] = qp.fromHartoEv(diff_2_3)
dictio['lab'] = 'Diff 2 3'
qp.writeH5fileDict(filename, dictio)
a = 0
for i in range(8):
a = a + 1
thisState = np.zeros_like(diff_0_1_all) + 999
thisState[15:-15,15:-15,30:-30] = pot[15:-15,15:-15,30:-30,i]
dictio = {}
if save_pot_diff:
filename = '/home/alessio/IMPORTANTS/VISUALIZE_ENERGY_DIFFERENCE/Energy{:04}.h5'.format(a)
dictio['diff'] = qp.fromHartoEv(thisState)
dictio['lab'] = i
qp.writeH5fileDict(filename, dictio)
In [ ]:
In [9]:
from quantumpropagator import fromLabelsToFloats, labTranformA
phis_ext = labTranformA(data['phis'])
gams_ext = labTranformA(data['gams'])
thes_ext = labTranformA(data['thes'])
phiV_ext, gamV_ext, theV_ext = fromLabelsToFloats(data)
# take step
dphi = phis_ext[0] - phis_ext[1]
dgam = gams_ext[0] - gams_ext[1]
dthe = thes_ext[0] - thes_ext[1]
# take range
range_phi = phis_ext[-1] - phis_ext[0]
range_gam = gams_ext[-1] - gams_ext[0]
range_the = thes_ext[-1] - thes_ext[0]
phis = phis_ext[15:-15]
gams = gams_ext[15:-15]
thes = thes_ext[30:-30]
phiV = phiV_ext[15:-15]
gamV = gamV_ext[15:-15]
theV = theV_ext[30:-30]
header = ' Labels extr. internal extr. dq range\n'
string = 'Phi -> {:8.4f} {:8.4f} {:8.4f} {:8.4f} {:8.4f} {:8.4f}\nGam -> {:8.4f} {:8.4f} {:8.4f} {:8.4f} {:8.4f} {:8.4f}\nThe -> {:8.4f} {:8.4f} {:8.4f} {:8.4f} {:8.4f} {:8.4f}'
out = (header + string).format(phiV_ext[-1],phiV_ext[0],phis_ext[-1],phis_ext[0],dphi,range_phi,
gamV_ext[-1],gamV_ext[0],gams_ext[-1],gams_ext[0],dgam,range_gam,
theV_ext[-1],theV_ext[0],thes_ext[-1],thes_ext[0],dthe,range_the)
print(out)
In [10]:
nacs = data['smoCube']
# take out zeros
NACS = nacs[15:-15,15:-15,30:-30]
# select the two states
print(NACS.shape, nacs.shape)
pL, gL, tL, sL, dL, coorL = NACS.shape
In [17]:
#%%time
n=10
makeGraph = True
states_to_consider = 2
if makeGraph:
for s1 in range(states_to_consider):
for s2 in range(s1):
a = np.abs(NACS[:,:,:,s1,s2,0].flatten())
binZ = [0.0000000000000001, 0.0000001, 0.000001, 0.00001,0.0001,0.001,0.01,0.1]
# thing here is the integer where I plot the bar (x position)
thing = np.arange(len(binZ)-1)
label_names = [ '{}'.format(x) for x in binZ ]
counts, bins = np.histogram(a,bins=binZ)
fig, ax0 = plt.subplots(1,1)
ax0.bar(thing,counts)
plt.xticks(thing,label_names)
plt.title('Nacs values between states {} {}'.format(s1,s2))
for xy in zip(thing, counts):
ax0.annotate('{}'.format(xy[1]), xy=xy)
In [18]:
cart = 0
s1 = 5
s2 = 4
p = 22
g=5
t=77
elem = np.abs(NACS[p,g,t,s1,s2,cart])
neighbors = np.abs(np.array([NACS[p+1,g,t,s1,s2,cart],
NACS[p-1,g,t,s1,s2,cart],
NACS[p,g+1,t,s1,s2,cart],
NACS[p,g-1,t,s1,s2,cart],
NACS[p,g,t+1,s1,s2,cart],
NACS[p,g,t-1,s1,s2,cart]]))
lol = neighbors - elem
differences = np.amin(lol)
print('{} {} {} {}'.format(elem, neighbors, lol, differences))
print('States({},{}) -> Cube({:2},{:2},{:2}): {:5.3e}'.format(s1,s2,p,g,t,differences))
In [19]:
NACS
cart = 0
for s1 in range(sL):
for s2 in range(s1):
#for p in qp.log_progress(range(pL),every=1,size=(pL)):
for p in range(1,pL-1):
for g in range(1,gL-1):
for t in range(1,tL-1):
elem = np.abs(NACS[p,g,t,s1,s2,cart])
neighbors = np.abs(np.array([NACS[p+1,g,t,s1,s2,cart],
NACS[p-1,g,t,s1,s2,cart],
NACS[p,g+1,t,s1,s2,cart],
NACS[p,g-1,t,s1,s2,cart],
NACS[p,g,t+1,s1,s2,cart],
NACS[p,g,t-1,s1,s2,cart]]))
differences = neighbors - elem
#print('{} {} {}'.format(elem, neighbors, differences))
if np.all(differences > 0.0001):
print('States({},{}) -> Cube({:2},{:2},{:2}): {:5.3e}'.format(s1,s2,p,g,t,differences))
In [20]:
# AAA is the plane at which I want to study the "only" point
AAA = NACS[10,:,:,1,2,2]
gam_That, the_That = np.unravel_index(AAA.argmin(), AAA.shape)
10, gam_That, the_That
phis[10],gams[gam_That],thes[the_That]
Out[20]:
In [21]:
%%cython --annotate --compile-args=-fopenmp --link-args=-fopenmp --force
### #%%cython
### #%%cython --annotate
import numpy as np
cimport numpy as np
cimport cython
from cython.parallel import prange
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@cython.nonecheck(False)
cdef void neigh(double [:,:,:,:] nacs_2_1, double [:,:,:,::1] bigger_2_1_nacs):
cdef:
int pL = 25, gL = 26, tL = 100,coorL=3
int p1,g1,t1,p2,g2,t2,coor,tuplL,pg
double thresh
thresh = 0.0001
tuplL = pL*gL
for coor in range(coorL):
#for pg in prange(tuplL, nogil=True, schedule='dynamic',num_threads=16):
for pg in range(tuplL):
p1 = pg // gL
g1 = pg % gL
for t1 in range(tL):
# for p2 in range(pL):
# for g2 in range(gL):
# for t2 in range(tL):
# if abs(nacs_2_1[p1,g1,t1,coor]) < 0.000001:
# bigger_2_1_nacs[p1,g1,t1,coor] = nacs_2_1[p1,g1,t1,coor]*100
# elif abs(nacs_2_1[p1,g1,t1,coor]) < 0.00001:
bigger_2_1_nacs[p1,g1,t1,coor] = nacs_2_1[p1,g1,t1,coor]*100
# return(bigger_2_1_nacs)
def neighbor(nacs_2_1,bigger_2_1_nacs):
return np.asarray(neigh(nacs_2_1,bigger_2_1_nacs))
print('done')
Out[21]:
In [22]:
%%time
state1 = 2
state2 = 1
nacs_2_1 = NACS[:,:,:,state1,state2,:]
nacs_other = NACS[:,:,:,state2,state1,:]
print(np.all(nacs_2_1 == -nacs_other))
print(nacs_2_1.shape)
bigger_2_1_nacs = np.zeros_like(nacs_2_1)
neighbor(nacs_2_1,bigger_2_1_nacs)
saveFile = False
dictio = {}
a=0
if saveFile:
for coord in range(3):
filename = '/home/alessio/k-nokick/IMPORTANTS/VISUALIZE_NACS/newNacSmoother/Nac{:04}.h5'.format(a)
# I reput the zeros out.
external = np.pad(vector[:,:,:,coord], ((15,15),(15,15),(30,30)), 'constant')
dictio['NACS'] = external
dictio['state1'] = state1
dictio['state2'] = state2
dictio['cart'] = coord
qp.writeH5fileDict(filename, dictio)
a += 1
filename = '/home/alessio/k-nokick/IMPORTANTS/VISUALIZE_NACS/newNacSmoother/Nac{:04}.h5'.format(a)
dictio['NACS'] = np.pad(nacs_2_1[:,:,:,coord], ((15,15),(15,15),(30,30)), 'constant')
dictio['state1'] = state1
dictio['state2'] = state2
dictio['cart'] = coord
qp.writeH5fileDict(filename, dictio)
a += 1
In [23]:
# PUT TRUE IF YOU WANT TO EXAGERATE AND CHANGE THE NACS
do_it = False
nacs_2_1 = nacs[:,:,:,1,2,:]
bigger_2_1_nacs = np.empty_like(nacs_2_1)
#print(bigger_2_1_nacs.shape)
pL,gL,tL,coorL = bigger_2_1_nacs.shape
for p in qp.log_progress(pL,every=1,size=(len(pL))):
for g in range(gL):
for t in range(tL):
for coor in range(coorL):
elem = nacs_2_1[p,g,t,coor]
if np.abs(elem) > 0.0001:
first = 2
secon = 4
# proximity(6)
bigger_2_1_nacs[p+1,g,t,coor] = elem/first
bigger_2_1_nacs[p-1,g,t,coor] = elem/first
bigger_2_1_nacs[p,g+1,t,coor] = elem/first
bigger_2_1_nacs[p,g-1,t,coor] = elem/first
bigger_2_1_nacs[p,g,t+1,coor] = elem/first
bigger_2_1_nacs[p,g,t-1,coor] = elem/first
# Corners (8)
bigger_2_1_nacs[p+1,g+1,t+1,coor] = elem/secon # 000
bigger_2_1_nacs[p+1,g+1,t-1,coor] = elem/secon
bigger_2_1_nacs[p+1,g-1,t+1,coor] = elem/secon
bigger_2_1_nacs[p+1,g-1,t-1,coor] = elem/secon # 011
bigger_2_1_nacs[p-1,g+1,t+1,coor] = elem/secon # 000
bigger_2_1_nacs[p-1,g+1,t-1,coor] = elem/secon
bigger_2_1_nacs[p-1,g-1,t+1,coor] = elem/secon
bigger_2_1_nacs[p-1,g-1,t-1,coor] = elem/secon # 011
# Half sides (12)
bigger_2_1_nacs[p+1,g,t+1,coor] = elem/secon
bigger_2_1_nacs[p+1,g,t-1,coor] = elem/secon
bigger_2_1_nacs[p-1,g,t+1,coor] = elem/secon
bigger_2_1_nacs[p-1,g,t-1,coor] = elem/secon
bigger_2_1_nacs[p+1,g+1,t,coor] = elem/secon
bigger_2_1_nacs[p+1,g-1,t,coor] = elem/secon
bigger_2_1_nacs[p-1,g+1,t,coor] = elem/secon
bigger_2_1_nacs[p-1,g-1,t,coor] = elem/secon
bigger_2_1_nacs[p,g+1,t+1,coor] = elem/secon
bigger_2_1_nacs[p,g+1,t-1,coor] = elem/secon
bigger_2_1_nacs[p,g-1,t+1,coor] = elem/secon
bigger_2_1_nacs[p,g-1,t-1,coor] = elem/secon
# 2 distant (6)
bigger_2_1_nacs[p+2,g,t,coor] = elem/secon
bigger_2_1_nacs[p-2,g,t,coor] = elem/secon
bigger_2_1_nacs[p,g+2,t,coor] = elem/secon
bigger_2_1_nacs[p,g-2,t,coor] = elem/secon
bigger_2_1_nacs[p,g,t+2,coor] = elem/secon
bigger_2_1_nacs[p,g,t-2,coor] = elem/secon
#print('{} {} {} {} {}'.format(p,g,t,coor,elem))
else:
bigger_2_1_nacs[p,g,t,coor] = elem
if do_it:
data_new = data
name_data_file_new = 'NAC_2_1_little_exagerated.pickle'
print(data_new.keys())
nacs[:,:,:,1,2,:] = bigger_2_1_nacs
nacs[:,:,:,2,1,:] = -bigger_2_1_nacs
data_new['smoCube'] = nacs
pickle.dump( data_new, open( name_data_file_new, "wb" ) )
In [11]:
dipo = data['dipCUBE']
DIPO = dipo[15:-15,15:-15,30:-30]
dipo.shape, DIPO.shape
Out[11]:
In [12]:
plt.close('all')
def do3dplot(xs,ys,zss):
'with mesh function'
fig = plt.figure(figsize=(9,9))
ax = fig.add_subplot(111, projection='3d')
X,Y = np.meshgrid(ys,xs)
#ax.set_zlim(-1, 1)
#ax.scatter(X, Y, zss)
ax.plot_surface(X, Y, zss,cmap=cm.coolwarm, linewidth=1, antialiased=False)
fig.canvas.layout.height = '800px'
fig.tight_layout()
def visualize_this_thing(thing,state1,state2,cart,kind,dim):
along = ['X','Y','Z']
print('DIPOLE between state ({},{}) along {} - Doing cut in {} with value ({:8.4f},{:8.4f}) - shape: {}'.format(state1,
state2,
along[cart],
kind,
dimV[kind][dim],
dims[kind][dim],
thing.shape))
if kind == 'Phi':
pot = thing[dim,:,:,cart,state1,state2]
print('Looking at DIPOLE with indexes [{},:,:,{},{},{}]'.format(dim,cart,state1,state2))
do3dplot(gams,thes,pot)
elif kind == 'Gam':
print('Looking at DIPOLE with indexes [:,{},:,{},{},{}]'.format(dim,cart,state1,state2))
pot = thing[:,dim,:,cart,state1,state2]
do3dplot(phis,thes,pot)
elif kind == 'The':
print('Looking at DIPOLE with indexes [:,:,{},{},{},{}]'.format(dim,cart,state1,state2))
pot = thing[:,:,dim,cart,state1,state2]
do3dplot(phis,gams,pot)
dimV = { 'Phi': phiV, 'Gam': gamV, 'The': theV } # real values
dims = { 'Phi': phis, 'Gam': gams, 'The': thes } # for labels
kinds = ['Phi','Gam','The']
def fun_pot2D(kind,state1, state2, cart,dim):
visualize_this_thing(DIPO, state1, state2, cart, kind, dim)
def nested(kinds):
dimensionV = dimV[kinds]
interact(fun_pot2D, kind=fixed(kinds),
state1 = widgets.IntSlider(min=0,max=7,step=1,value=0,continuous_update=False),
state2 = widgets.IntSlider(min=0,max=7,step=1,value=1,continuous_update=False),
cart = widgets.IntSlider(min=0,max=2,step=1,value=2,continuous_update=False),
dim = widgets.IntSlider(min=0,max=(len(dimensionV)-1),step=1,value=0,continuous_update=False))
interact(nested, kinds = ['Gam','Phi','The']);
In [49]:
import ipyvolume as ipv
def do3dplot2(xs,ys,zss):
X,Y = np.meshgrid(ys,xs)
ipv.figure()
ipv.plot_surface(X, zss, Y, color="orange")
ipv.plot_wireframe(X, zss, Y, color="red")
ipv.show()
In [13]:
pot = data['potCube'] - data['potCube'].min()
A = pot
# find the minimum index having the shape
phi_min, gam_min, the_min, state_min = np.unravel_index(A.argmin(), A.shape)
phi_min, gam_min, the_min, state_min
Out[13]:
In [14]:
nacs.shape
Out[14]:
In [15]:
B = nacs[:,:,:,:,1,0]
In [16]:
# this should be absolute value
phi_ci, gam_ci, the_ci, cart_ci = np.unravel_index(B.argmax(), B.shape)
np.unravel_index(B.argmax(), B.shape)
Out[16]:
In [17]:
phis_ext[16],gams_ext[15],thes_ext[112]
Out[17]:
In [18]:
phi_ci, gam_ci, the_ci = [16,15,112]
In [19]:
# I start by making a cube of 0.
# boolean that creates (and overwrite) the file
create_mask_file = False
ZERO = np.zeros_like(pot[:,:,:,0])
print(ZERO.shape)
region_a = np.zeros_like(ZERO)
region_b = np.zeros_like(ZERO)
region_c = np.zeros_like(ZERO)
# for pi, p in qp.log_progress(enumerate(phis_n),every=1,size=(len(phis_n))):
# p_linea_a = 18
# g_linea_S = 30
# t_linea_b = 112
p_linea_a, g_linea_S, t_linea_b = 18, 30, 112
m_coeff,q_coeff = 1.7,75
for p,phi in qp.log_progress(enumerate(phis_ext),every=1,size=(len(phis_ext))):
for g,gam in enumerate(gams_ext):
lineValue_theta = m_coeff * g + q_coeff
for t,the in enumerate(thes_ext):
if p > p_linea_a:
region_a[p,g,t] = 1
if p <= p_linea_a:
if t > lineValue_theta:
region_c[p,g,t] = 1
else:
region_b[p,g,t] = 1
# if t > t_linea_b and g < g_linea_S:
# region_c[p,g,t] = 1
# else:
# region_b[p,g,t] = 1
# to paint cubes on the verge I make zero the sides values
if p==0 or g == 0 or t == 0 or p == len(phis_ext)-1 or g == len(gams_ext)-1 or t == len(thes_ext)-1:
region_a[p,g,t] = 0
region_b[p,g,t] = 0
region_c[p,g,t] = 0
regions = [{'label' : 'FC', 'cube': region_a},{'label' : 'reactants', 'cube': region_b},{'label' : 'products', 'cube': region_c}]
if create_mask_file:
print('I created the regions pickle file')
pickle.dump(regions, open('regions.pickle', "wb" ) )
else:
qp.warning("file region NOT written, check the variable 'create_mask_file' if you want to write region file")
In [21]:
# I start by making a cube of 0.
# boolean that creates (and overwrite) the file
create_adv_mask_files = False
ZERO = np.zeros_like(pot[:,:,:,0])
ones = np.ones_like(pot[:,:,:,0])
print(ZERO.shape)
mask_CI = np.zeros_like(ZERO)
mask_AFC = np.zeros_like(ZERO)
# 18, 30, 112
#p_cube1, p_cube2, g_cube1, g_cube2, t_cube1, t_cube2 = 12, 29, 15, 32, 82, 142
p_CI_cube1, p_CI_cube2, g_CI_cube1, g_CI_cube2, t_CI_cube1, t_CI_cube2 = 12, 29, 15, 32, 82, 142
p_AFC_cube1, p_AFC_cube2, g_AFC_cube1, g_AFC_cube2, t_AFC_cube1, t_AFC_cube2 = 0, 54, 0, 55, 85, 159
for p,phi in qp.log_progress(enumerate(phis_ext),every=1,size=(len(phis_ext))):
for g,gam in enumerate(gams_ext):
for t,the in enumerate(thes_ext):
if p > p_CI_cube1 and p < p_CI_cube2 and g > g_CI_cube1 and g < g_CI_cube2 and t > t_CI_cube1 and t < t_CI_cube2:
mask_CI[p,g,t] = 1
if p > p_AFC_cube1 and p < p_AFC_cube2 and g > g_AFC_cube1 and g < g_AFC_cube2 and t > t_AFC_cube1 and t < t_AFC_cube2:
mask_AFC[p,g,t] = 1
# to paint cubes on the verge I make zero the sides values
if p==0 or g == 0 or t == 0 or p == len(phis_ext)-1 or g == len(gams_ext)-1 or t == len(thes_ext)-1:
ones[p,g,t] = 0
mask_AFC[p,g,t] = 0
masks_only_one = [{'label' : 'All', 'cube': ones, 'states':[0,1,2,3,4,5,6,7], 'show' : False}]
masks = [
{'label' : 'All', 'cube': ones, 'states':[0,1,2,3,4,5,6,7], 'show' : False},
{'label' : 'AfterFC', 'cube': mask_AFC, 'states':[0,1,2,3,4,5,6,7], 'show' : True},
{'label' : 'CI', 'cube': mask_CI, 'states':[0,1,2,3,4,5,6,7], 'show' : True},
]
# {'label' : 'Mask CI', 'cube': mask_a, 'states':[0,1]}]
#masks = [{'label' : 'Mask', 'cube': ones, 'states':[0,1,2,3,4,5,6,7]}]
if create_adv_mask_files:
print('I created the regions pickle file')
pickle.dump(masks_only_one, open('advanced_masks_onlyONE.pickle', "wb" ))
pickle.dump(masks, open('advanced_masks.pickle', "wb" ))
else:
qp.warning("file advanced_masks NOT written, check the variable 'create_adv_mask_files' if you want to write region file")
In [22]:
# dipo_min = dipo[phi_min, gam_min, the_min]
# dipo_ci = dipo[phi_ci, gam_ci, the_ci]
# difference_dipo = dipo_ci - dipo_min
# for i in range(8):
# permanent = difference_dipo[:,i,i]
# print('S_{} -> {}'.format(i,permanent))
In [23]:
# dipo_min[:,1,2],dipo_min[:,0,1],dipo_min[:,0,6],dipo_min[:,0,3],dipo_min[:,0,2],dipo_min[:,0,7]
In [24]:
# npy = '/home/alessio/Desktop/NAC_CORRECTION_NOVEMBER2018/dataprova.npy'
# dictio = np.load(npy)[()]
# dictio.keys()
In [25]:
# NACS2 = dictio['nacCUBE']
# NACS2.shape
In [26]:
# def do3dplot(xs,ys,zss):
# 'with mesh function'
# fig = plt.figure(figsize=(9,9))
# ax = fig.add_subplot(111, projection='3d')
# X,Y = np.meshgrid(ys,xs)
# #ax.set_zlim(-1, 1)
# #ax.scatter(X, Y, zss)
# ax.plot_wireframe(X, Y, zss)
# fig.tight_layout()
# def visualize_this_thing(thing,state1,state2,cart,kind,dim):
# print(thing.shape)
# print('\nWARNING, this is not fully correct!!! Not SMO and not really what you think\n')
# along = ['Phi','Gam','The']
# print('NAC between state ({},{}) along {}\nDoing cut in {} with value ({:8.4f},{:8.4f})'.format(state1,
# state2,
# along[cart],
# kind,
# dimV[kind][dim],
# dims[kind][dim]))
# if kind == 'Phi':
# pot = thing[dim,:,:,state1,state2,0,cart]
# print('\nLooking at SMO with indexes [{},:,:,{},{},{}]'.format(dim, state1,state2,cart))
# do3dplot(gams,thes,pot)
# elif kind == 'Gam':
# print('\nLooking at SMO with indexes [:,{},:,{},{},{}]'.format(dim, state1,state2,cart))
# pot = thing[:,dim,:,state1,state2,0,cart]
# do3dplot(phis,thes,pot)
# elif kind == 'The':
# print('\nLooking at SMO with indexes [:,:,{},{},{},{}]'.format(dim, state1,state2,cart))
# pot = thing[:,:,dim,state1,state2,0,cart]
# do3dplot(phis,gams,pot)
# dimV = { 'Phi': phiV, 'Gam': gamV, 'The': theV } # real values
# dims = { 'Phi': phis, 'Gam': gams, 'The': thes } # for labels
# kinds = ['Phi','Gam','The']
# def fun_pot2D(kind,state1, state2, cart,dim):
# visualize_this_thing(NACS2, state1, state2, cart, kind, dim)
# def nested(kinds):
# dimensionV = dimV[kinds]
# interact(fun_pot2D, kind=fixed(kinds), state1 = widgets.IntSlider(min=0,max=7,step=1,value=0), state2 = widgets.IntSlider(min=0,max=7,step=1,value=0), cart = widgets.IntSlider(min=0,max=2,step=1,value=0), dim = widgets.IntSlider(min=0,max=(len(dimensionV)-1),step=1,value=0))
# interact(nested, kinds = ['Phi','Gam','The']);
In [27]:
# print(data.keys())
# data_new = data
# nacs_new = data['smoCube']
# NACS_new = nacs_new[15:-15,15:-15,30:-30]
# print(NACS_new.shape,nacs_new.shape)
# phi_ext_000_000 = 29
# phi_prev = 28
# phi_next = 30
# new_nacs = np.copy(nacs)
# for g in range(56):
# for t in range(160):
# not_correct = nacs[phi_ext_000_000,g,t]
# correct_prev = nacs[phi_prev ,g,t]
# correct_next = nacs[phi_next ,g,t]
# #if np.linalg.norm(not_correct) > 0.001:
# # print('{} {}\nThis {} \nMiddle {}\n After {}'.format(g,t,correct_prev[:,:,1], not_correct[:,:,1],correct_next[:,:,1]))
# for state1 in range(8):
# for state2 in range(8):
# for cart in range(3):
# value_prev = correct_prev[state1,state2,cart]
# value_this = not_correct [state1,state2,cart]
# value_next = correct_next[state1,state2,cart]
# average = (value_prev + value_next)/2
# if np.sign(average) == np.sign(value_this):
# new_value = value_this
# else:
# new_value = -value_this
# new_nacs[phi_ext_000_000,g,t,state1,state2,cart] = new_value
In [28]:
# def do3dplot(xs,ys,zss):
# 'with mesh function'
# fig = plt.figure(figsize=(9,9))
# ax = fig.add_subplot(111, projection='3d')
# X,Y = np.meshgrid(ys,xs)
# #ax.set_zlim(-1, 1)
# #ax.scatter(X, Y, zss)
# ax.plot_wireframe(X, Y, zss)
# fig.tight_layout()
# def visualize_this_thing(thing,state1,state2,cart,kind,dim):
# print(thing.shape)
# along = ['Phi','Gam','The']
# print('NAC between state ({},{}) along {}\nDoing cut in {} with value ({:8.4f},{:8.4f})'.format(state1,
# state2,
# along[cart],
# kind,
# dimV[kind][dim],
# dims[kind][dim]))
# if kind == 'Phi':
# pot = thing[dim,:,:,state1,state2,cart]
# print('\nLooking at SMO with indexes [{},:,:,{},{},{}]'.format(dim, state1,state2,cart))
# do3dplot(gams_ext,thes_ext,pot)
# elif kind == 'Gam':
# print('\nLooking at SMO with indexes [:,{},:,{},{},{}]'.format(dim, state1,state2,cart))
# pot = thing[:,dim,:,state1,state2,cart]
# do3dplot(phis_ext,thes_ext,pot)
# elif kind == 'The':
# print('\nLooking at SMO with indexes [:,:,{},{},{},{}]'.format(dim, state1,state2,cart))
# pot = thing[:,:,dim,state1,state2,cart]
# do3dplot(phis_ext,gams_ext,pot)
# dimV = { 'Phi': phiV_ext, 'Gam': gamV_ext, 'The': theV_ext } # real values
# dims = { 'Phi': phis_ext, 'Gam': gams_ext, 'The': thes_ext } # for labels
# kinds = ['Phi','Gam','The']
# def fun_pot2D(kind,state1, state2, cart,dim):
# visualize_this_thing(new_nacs, state1, state2, cart, kind, dim)
# def nested(kinds):
# dimensionV = dimV[kinds]
# interact(fun_pot2D, kind=fixed(kinds), state1 = widgets.IntSlider(min=0,max=7,step=1,value=0), state2 = widgets.IntSlider(min=0,max=7,step=1,value=0), cart = widgets.IntSlider(min=0,max=2,step=1,value=0), dim = widgets.IntSlider(min=0,max=(len(dimensionV)-1),step=1,value=0))
# interact(nested, kinds = ['Phi','Gam','The']);
In [29]:
#name_data_file_new = 'newExtrapolated_allCorrectionSECOND.pickle'
#data_new.keys()
In [30]:
# data_new['smoCube'] = new_nacs
# pickle.dump( data_new, open( name_data_file_new, "wb" ) )
In [31]:
folder = '.'
a=0
saveFile = False
for state1 in range(8):
for state2 in range(state1):
for cart in range(3):
dictio = {}
cartL = ['X','Y','Z']
print('Nacs ({},{}) along {} -> {:04}'.format(state1,state2,cartL[cart],a))
a+=1
if saveFile:
filename = 'Nac{:04}.h5'.format(a)
dictio['NACS'] = nacs[:,:,:,state1,state2,cart]
dictio['state1'] = state1
dictio['state2'] = state2
dictio['cart'] = cart
qp.writeH5fileDict(filename, dictio)
In [32]:
# lol2 is the new function to be added
phi_index = 16
theta_index = 81
state_index = 0
lol = pot[phi_index,:,theta_index,state_index]
num = 15
constant = 0.001
lol2 = np.zeros_like(lol)
for i in range(num):
lol2[i] = constant * (i-num)**2
#print('{} {} {}'.format(i,num,i-num))
In [33]:
fig = plt.figure()
plt.title('Gamma wall')
plt.xlabel('Gamma')
plt.ylabel('Energy')
plt.plot(lol)
plt.plot(lol2+lol);
In [33]:
newpot = np.zeros_like(pot)
for p in range(55):
for t in range(160):
for s in range(8):
newpot[p,:,t,s] = pot[p,:,t,s] + lol2
In [34]:
do_it = False
if do_it:
data_new = data
name_data_file_new = 'newExtrapolated_gammaExtrExag.pickle'
data_new.keys()
data_new['potCube'] = newpot
pickle.dump( data_new, open( name_data_file_new, "wb" ) )
else:
qp.warning('Here it is set to false, new file is NOT created')
In [35]:
fig = plt.figure()
phi = 20
the = 100
plt.plot(pot[phi,:,the,1])
plt.plot(newpot[phi,:,the,1]);
In [ ]: