<jupyter_start><jupyter_text>Imports<jupyter_code>%matplotlib inline
from __future__ import division

import tools.convolution
from tools.readtree import readtree

from tqdm import tqdm

import pandas as pd
import numpy as np
from tqdm import tqdm
import itertools
import matplotlib.pyplot as plt

import matplotlib
from mpl_toolkits.mplot3d import Axes3D
import plotly.offline as py
import plotly.graph_objs as GO

import astropy
import os
from scipy.signal import argrelextrema

from tools.extrema import extrema
from tools.convolution import convolution
import glob

py.init_notebook_mode()
matplotlib.rcParams.update(matplotlib.rcParamsDefault)
matplotlib.style.use(['beamer', 'ggplot'])

def savefig(figure, path):
    figure.savefig(path+'.png')
    figure.savefig(path+'.pdf')
%load_ext autotime
%load_ext line_profiler<jupyter_output><empty_output><jupyter_text>Load some data about halos<jupyter_code>from tools.infos import read_infos
infos = read_infos('/data52/Horizon-AGN/OUTPUT_DIR/output_00002/info_00002.txt')['infos']
inertia_data = pd.read_hdf('data/halo_kin.00002.1e10_m.hdf')
print(str(inertia_data.index.size) + ' halos!')<jupyter_output>280219 halos!
time: 412 ms<jupyter_text>Smoothing tree  Set halo number<jupyter_code>halo_i = 536<jupyter_output>time: 1.07 ms<jupyter_text>Load data<jupyter_code>from scipy.io import FortranFile as FF
ff = FF('smoothingtrees/11.5_logm_12_both/halo_%s-centers.bin' % halo_i)
ndim, noutput = ff.read_ints()
outputs = ff.read_ints()
centers = ff.read_reals().reshape((ndim, noutput))
outputs, centers
df = pd.read_csv('smoothingtrees/11.5_logm_12_both/halo_%s_around' % halo_i, delim_whitespace=True, engine='c', index_col=0)
df[['x', 'y', 'z']].describe()
dfin = pd.read_csv('smoothingtrees/11.5_logm_12_both/halo_%s-00002' % halo_i, delim_whitespace=True, engine='c', index_col=0)
desc = dfin.describe()

center = desc.loc['mean'][['x', 'y', 'z']]
span = desc.loc['max'][['x', 'y', 'z']] - desc.loc['min'][ ['x', 'y', 'z']]

desc
left = center - 1.5*span
right = center + 1.5*span
df = df[df.x > left.x][df.y > left.y][df.z > left.z][df.x < right.x][df.y < right.y][df.z < right.z]
df.describe()<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>plt.figure()
ax = plt.subplot(111, projection='3d')
each = 1000
ax.scatter3D(df.x[::each], df.y[::each], df.z[::each])
plt.show()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter3D(df.x[::each], df.y[::each], df.z[::each], alpha=0.5)
ax.scatter3D(dfin.x, dfin.y, dfin.z,  c='red')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

df[['x', 'y', 'z']].mean(), dfin[['x', 'y', 'z']].mean()
plt.show()
df.quantile(0.01), df.quantile(0.99)
plt.close('all')<jupyter_output>time: 1.27 ms<jupyter_text>Home made!  Convolution<jupyter_code>def f():
    pixs, pos, eigvect, eigvals, peaktype = extrema.extrema_get(ndim, npeak)
    #convert back pixs into
    index = np.array(np.unravel_index(pixs, conv.shape, order='F'))
    return dict(
        index = index,
        pos = pos,
        eigvect = eigvect,
        eigvals = eigvals,
        type = peaktype
    )

def as_df(peaks):
    df = pd.DataFrame(peaks['pos'].T, columns=['x', 'y', 'z'])
    df['i'] = df['j'] = df['k'] = 0
    df[['i', 'j', 'k']] = peaks['index'].T
    df['type'] = peaks['type']
    
    df['xr'] = df['i'] / edges[0].shape[0] * (edges[0].max() - edges[0].min()) + edges[0].min()
    df['yr'] = df['j'] / edges[1].shape[0] * (edges[1].max() - edges[1].min()) + edges[1].min()
    df['zr'] = df['k'] / edges[2].shape[0] * (edges[2].max() - edges[2].min()) + edges[2].min()
    
    return df
from tools.convolution import convolution
bins = 64
dens, edges = convolution.conv_density(df[['x', 'y', 'z']].as_matrix().T, bins)
gaussian = convolution.kernel_gaussian3d(bins, 2)
gfft = convolution.fft(gaussian)
hfft = convolution.fft(dens)
cfft = convolution.conv_prod(hfft, gfft)

conv = convolution.ifft(cfft)
from pymses.utils import constants as C
(infos.loc['unit_l'].value * C.cm).express(C.Mpc)
extrema.npeaks = bins**3
extrema.nbins = np.array([bins, bins, bins])
extrema.npeaks, extrema.nbins

ndim, npeak = extrema.extrema_compute(np.asfortranarray(conv))
ndim, npeak

peaks = f()
dfp = as_df(peaks)
%matplotlib inline
from pymses.utils import constants as C
from matplotlib.colors import LogNorm
with matplotlib.style.context(['seaborn-paper']):
    toMpc = (infos.loc['unit_l'].value * C.cm).express(C.Mpc)
    cmean = conv.mean()
    dens_mean = (conv - cmean) / cmean * 8e7 #Msun/DM particle
    cent = 28
    around = 2
    plt.imshow(dens_mean[:,:,cent-around:cent+around].sum(axis=2).T, origin='lower', 
               extent=(edges[0, 0]*toMpc, edges[0, -1]*toMpc, edges[1, 0]*toMpc, edges[1, -1]*toMpc), cmap='viridis',
               interpolation='none', aspect='auto')
    plt.xlim(edges[0, 0]*toMpc, edges[0, -1]*toMpc)
    plt.ylim(edges[1, 0]*toMpc, edges[1, -1]*toMpc)
    
    plt.xlabel('$x$ [$\mathrm{Mpc}/h$]')
    plt.ylabel('$y$ [$\mathrm{Mpc}/h$]')
    
    plt.title(u'$z=%.2f$, épaisseur $\Delta z = %.2f\ \mathrm{Mpc}/h$' % (1/(infos.loc['aexp']) - 1, around*2 / bins * toMpc))
    
    cb = plt.colorbar()
    cb.set_label(u'Surdensité ' + r'$\delta = \frac{\rho - \bar{\rho}}{\bar{\rho}}$')
        
    a = dfp[dfp.z >= cent-around][dfp.z <= cent+around][dfp.type == 4]
    plt.scatter(a.xr*toMpc, a.yr*toMpc, marker='+', c='red', label='Maximum')
    
    a = dfp[dfp.z >= cent-around][dfp.z <= cent+around][dfp.type == 3]
    plt.scatter(a.xr*toMpc, a.yr*toMpc, marker='.', c='green', label='Point selle')
    
    plt.grid('off')
    plt.legend(loc='best', framealpha=0)
    plt.tight_layout()
    plt.savefig('/home/cadiou/plots/paper/halo_%s_density.pdf'%halo_i)
dfp[dfp.type==4].plot.scatter('x', 'z')<jupyter_output><empty_output><jupyter_text>Find extrema<jupyter_code>print(conv.shape)
extrema.npeaks = bins**3
extrema.nbins = np.array([bins, bins, bins])
extrema.npeaks, extrema.nbins
ndim, npeak = extrema.extrema_compute(np.asfortranarray(conv))
ndim, npeak
peaks = f()
dfp = as_df(peaks)
dfp[dfp.k == 32][dfp.type == 4]
plt.pcolormesh(edges[0], edges[1], conv[:,:,5].T)

colors = ['red', 'blue', 'green', 'purple']
for i, (t, l) in enumerate(zip(range(1, 5), ['min', 'saddle (1 max)', 'saddle (2 max)', 'max'])):
    mask = (dfp.type == t)
    plt.scatter(dfp[mask].xr, dfp[mask].yr, c=colors[i], label=l)
    print('%i of type %s' % (dfp[mask].index.size, l))
plt.legend()

plt.xlim()<jupyter_output>70 of type min
246 of type saddle (1 max)
220 of type saddle (2 max)
79 of type max<jupyter_text>All at once!<jupyter_code>import cython
'''%%cython
# estimate density
from tools.convolution import convolution
from tools.extrema import extrema
cimport numpy as np
import numpy as np
import pandas as pd
from tqdm import tqdm

import cython

@cython.boundscheck(False)
@cython.wraparound(False)'''
def treat(df, bins=64,steps=5):
    #cdef np.ndarray[double,  ndim=3] dens, gaussian
    '''    cdef np.ndarray[double,  ndim=3] dens, gaussian
    cdef np.ndarray[complex, ndim=3] dfft, gfft
    cdef np.ndarray[double,  ndim=2] edges
    cdef np.ndarray[long,    ndim=2] index
    
    cdef np.ndarray[int] peaktype, pixs
    cdef np.ndarray[double, ndim=2] pos, eigvect
    cdef np.ndarray[double] eigvals
    cdef double sigma'''
    
    dens, edges = convolution.conv_density(df[['x', 'y', 'z']].as_matrix().T, bins)
    dfft = convolution.fft(dens)

    sigmas_list = np.linspace(1, bins // 2, steps)

    # peak finding
    extrema.npeaks = bins**3
    extrema.nbins = np.array([bins, bins, bins])
    extrema.nproc = 1
    dfpeak = pd.DataFrame()
    
    for sigma in tqdm(sigmas_list):
        # convolve
        gaussian = convolution.kernel_gaussian3d(bins, sigma)

        gfft = convolution.fft(gaussian)
        cfft = convolution.conv_prod(dfft, gfft)
        conv = convolution.ifft(cfft)

        ndim, npeak = extrema.extrema_compute(np.asfortranarray(conv))

        pixs, pos, eigvect, eigvals, peaktype = extrema.extrema_get(ndim, npeak)
        #convert back pixs into
        index = np.array(np.unravel_index(pixs, conv.shape, order='F'))

        dfpeak = dfpeak.append(pd.DataFrame(dict(
                    i=index[0, :], j=index[1, :], k=index[2, :], type=peaktype, eigval= eigvals,
                    x=pos[0], y=pos[1], z=pos[2],
                    evx= eigvect[0], evy=eigvect[1], evz=eigvect[2],
                    sigma=sigma
                    )))
        
    return dfpeak, dens, edges
bins = 32
dfpeak, dens, edges = treat(df, bins=bins, steps=40)
dfpeak['amp'] = dfpeak.eval('evx**2 + evy**2 + evz**2')
print(str(dfpeak.index.size) + ' elements')
dfpeak.head()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

m = dfpeak.type==4

def mapon(lbound, ubound, maxval, minval, arr):
    maxarr = np.max(arr)
    minarr = np.min(arr)
    return (arr - minval)/(maxval - minval)*ubound + lbound


mmapon= lambda arr: mapon(5, 50, np.max(np.sqrt(1/dfpeak['amp'])), np.min(np.sqrt(1/dfpeak['amp'])), arr)

m = dfpeak.type == 3
ax.scatter3D(dfpeak[m]['x'], dfpeak[m]['y'], dfpeak[m]['sigma'], s=mmapon(np.sqrt(1/dfpeak[m]['amp'])), c='red', alpha=0.3, marker='x',
            label='Saddle point, 2 max')

m = dfpeak.type == 2
ax.scatter3D(dfpeak[m]['x'], dfpeak[m]['y'], dfpeak[m]['sigma'], s=mmapon(np.sqrt(1/dfpeak[m]['amp'])), c='green', alpha=0.1, marker='x',
            label='Saddle point, 1 max')

m = dfpeak.type == 4
ax.scatter3D(dfpeak[m]['x'], dfpeak[m]['y'], dfpeak[m]['sigma'], s=mmapon(np.sqrt(1/dfpeak[m]['amp'])), c='blue', marker='o',
            label='Peak')

ax.set_xlim(0, bins)
ax.set_ylim(0, bins)
ax.set_zlim(0, dfpeak['sigma'].max())

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('$\sigma$')

#ax.legend(loc='best')

fig.tight_layout()

savefig(plt.gcf(), '/home/cadiou/plots/smoothingtree/halo_%s' % halo_i)

ax.set_zlim(0, bins // 2)

savefig(plt.gcf(), '/home/cadiou/plots/smoothingtree/halo_%s_zoom' % halo_i)
plt.show()
fig = plt.figure()
ax = fig.add_subplot(111)

def mapon(lbound, ubound, maxval, minval, arr):
    maxarr = np.max(arr)
    minarr = np.min(arr)
    return (arr - minval)/(maxval - minval)*ubound + lbound


mmapon= lambda arr: mapon(5, 50, np.max(np.sqrt(1/dfpeak['amp'])), np.min(np.sqrt(1/dfpeak['amp'])), arr)

m = dfpeak.type == 3
ax.scatter(dfpeak[m]['x'],
           #dfpeak[m]['y'], 
           dfpeak[m]['sigma'], s=mmapon(np.sqrt(1/dfpeak[m]['amp'])), c='red', alpha=0.3, marker='x',
           label='Saddle point, 2 max')

m = dfpeak.type == 2
ax.scatter(dfpeak[m]['x'], 
           #dfpeak[m]['y'],
           dfpeak[m]['sigma'], s=mmapon(np.sqrt(1/dfpeak[m]['amp'])), c='green', alpha=0.1, marker='x',
           label='Saddle point, 1 max')

m = dfpeak.type == 4
ax.scatter(dfpeak[m]['x'], 
           #dfpeak[m]['y'], 
           dfpeak[m]['sigma'], s=mmapon(np.sqrt(1/dfpeak[m]['amp'])), c='blue', marker='o',
           label='Peak')

ax.set_xlim(0, bins)
ax.set_ylim(0, dfpeak['sigma'].max())

ax.set_xlabel('x')
ax.set_ylabel('$\sigma$')

#ax.legend(loc='best')

fig.tight_layout()

savefig(plt.gcf(), '/home/cadiou/plots/smoothingtree/halo_%s_2D' % halo_i)

ax.set_ylim(0, bins // 2)

savefig(plt.gcf(), '/home/cadiou/plots/smoothingtree/halo_%s_zoom_2D' % halo_i)
plt.close('all')<jupyter_output>time: 1.25 ms
{'fraction_non_alphanumeric': 0.19417647579431807, 'fraction_numerical': 0.022391362067439597, 'mean_word_length': 4.1504102096627165, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '622943', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4969, 'n_tokens_neox': 4683, 'n_words': 840, 'path': 'Smoothing tree.ipynb', 'repo_name': 'cphyc/cosmo_z17to0', 'size': '11299'}