text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

#

#

#

#
Computational Seismology

#
 SBP-SAT finite difference method for the 1D elastic wave equation in first order form

#

#

#

# This notebook is based on the paper [Dynamic earthquake rupture simulations on nonplanar faults embedded in 3D geometrically complex, heterogeneous Earth models](https://pangea.stanford.edu/~edunham/publications/Duru_Dunham_FD3d_JCP16.pdf), and on the theory of summation-by-parts (SBP) finite difference methods and weak implementation of boundary conditions using the simultaneous-approximation-term (SAT).
#
#
# ##### Authors:
# * Kenneth Duru
#
# ---

# ## Basic Equations ##
#
# We consider in this notebook a system of hyperbolic PDE, the elastic wave equation, in 1D. The source-free elastic wave equation in a heterogeneous 1D medium is
#
# \begin{align}
# \rho(x)\frac{\partial v(x,t)}{\partial t} -\frac{\partial\sigma(x,t)}{\partial x} & = 0\\
# \frac{1}{\mu(x)}\frac{\partial\sigma(x,t)}{\partial t} -\frac{\partial v(x,t)}{\partial x} & = 0
# \end{align}
#
# with $\rho(x)$ the density, $\mu(x)$ the shear modulus and $x \in [0, L]$. At the boundaries $ x = 0, x = L$ we pose the general well-posed linear boundary conditions
#
# \begin{equation}
# \begin{split}
# B_0(v, \sigma, Z_{s}, r_0): =\frac{Z_{s}}{2}\left({1-r_0}\right){v} -\frac{1+r_0}{2} {\sigma} = 0,  \quad \text{at} \quad x = 0, \\
#  B_L(v, \sigma, Z_{s}, r_n): =\frac{Z_{s}}{2} \left({1-r_n}\right){v} + \frac{1+r_n}{2}{\sigma} = 0, \quad \text{at} \quad  x = L.
#  \end{split}
# \end{equation}
#
# with the reflection coefficients $r_0$, $r_n$ being real numbers and  $|r_0|, |r_n| \le 1$.
#
# Note that at $x = 0$,  while  $r_0 = -1$ yields a clamped wall, $r_0 = 0$  yields  an absorbing boundary, and  with $r_0 = 1$  we have a free-surface boundary condition. Similarly, at $x = L$, $r_n = -1$ yields a clamped wall, $r_n = 0$ yields an absorbing boundary, and  $r_n = 1$  gives a free-surface boundary condition.
#
# We introduce the mechanical energy defined by
# \begin{equation}
# E(t) = \int_0^L{\left(\frac{\rho(x)}{2} v^2(x, t) + \frac{1}{2\mu(x)}\sigma^2(x, t)\right) dx},
# \end{equation}
#
# where $E(t)$ is the sum of the kinetic energy and the strain energy.
# We have
#
# \begin{equation}
# \frac{d E(t)}{dt} = -v(0, t)\sigma(0, t) + v(L, t)\sigma(L, t) \le 0.
# \end{equation}
#
# From the boundary conditions, it is easy to check that  $v(0, t)\sigma(0, t) \ge 0$ and $v(L, t)\sigma(L, t) \le 0$, for all $|r_0|, |r_n| \le 1$. This energy loss through the boundaries is what the numerical method should emulate.
#
# 1) Discretize the spatial domain $x$ into $N$ discrete nodes with the uniform spatial step $\Delta{x} = L/(N-1)$, denote the unknown fields at the nodes: $\mathbf{v}\left(t\right) = [v_1\left(t\right), v_2\left(t\right), \cdots, v_N\left(t\right)]$, and $\boldsymbol{\sigma}\left(t\right) = [\sigma_1\left(t\right), \sigma_2\left(t\right), \cdots, \sigma_N\left(t\right)]$.
#
#
# 2) At the grid-point $x_j = (j-1)\Delta{x}$: Approximate the spatial derivative by a finite difference operator $\partial v/\partial x\Big|_{x = x_j} \approx \left(\mathbf{D}\mathbf{v}\right)_j $. Here $\mathbf{D}$ is a finite difference matrix satisfying the summation-by-parts property:
#
# \begin{align}
# \mathbf{H}\mathbf{D} = \mathbf{Q}, \quad \mathbf{Q} + \mathbf{Q} = \left(\boldsymbol{e}_{N}\boldsymbol{e}_{N}^T -\boldsymbol{e}_{1}\boldsymbol{e}_{1}^T\right), \quad \mathbf{H}^T = \mathbf{H} > 0,
# \end{align}
#
# where, $\boldsymbol{e}_{0} = [1, 0, \dots,  0 ]^T, \quad  \boldsymbol{e}_{L} = [ 0, 0, \dots,  1 ]^T$ and $\mathbf{H}$ defines a dicrete norm. We consider only diagonal norm SBP operators with $H_{jj} = h_j > 0$, and define the quadrature rule
#
# \begin{equation}
#  \sum_{i = 1}^{N} f(x_j)h_j \approx \int_{0}^{L}f(x) dx.
# \end{equation}
#
# The second order accurate SBP operator for the first derivative is:
# \begin{align}
# \left(\mathbf{D}\mathbf{v}\right)_j = \frac{v_{j+1}-v_{j-1}}{2 \Delta{x}}, \quad j = 2, 3, \cdots N-1, \quad
# \left(\mathbf{D}\mathbf{v}\right)_1 = \frac{v_{2}-v_{1}}{\Delta{x}},\quad
# \left(\mathbf{D}\mathbf{v}\right)_N = \frac{v_{N}-v_{N-1}}{\Delta{x}}, \quad j = N.
# \end{align}
#
# Note that the interior stencils are centered, with second order accuracy and the boundary stencils are one-sided and first order accurate.
#
# Higher order SBP operators can be found in the book: High Order Difference Methods for Time Dependent PDE, by B. Gustafsson. In this notebook we implement SBP operators with interior accuracy 2, 4 and 6. The implementation of the spatial derivative operators can be found in the file first_derivative_sbp_operators.py
#
# To construct a stable  semi-discrete approximation we replace the spatial derivatives by the SBP operators, and add the boundary conditions as SAT-terms with special penalty weights having:
#
# \begin{align}
# \frac{d \mathbf{v}(t)}{d t} = {\boldsymbol{\rho}}^{-1}
# \left(\mathbf{D} \boldsymbol{\sigma}(t) - \underbrace{\mathbf{H}^{-1}\left(\tau_{11}\boldsymbol{e}_{1}B_0\left(v_1, \sigma_1, Z_{s}, r_0\right) + \tau_{12}\boldsymbol{e}_{N}B_L\left(v_N, \sigma_N, Z_{s}, r_n\right)\right)}_{SAT \to 0}\right),
# \end{align}
#
# \begin{align}
# \frac{d \boldsymbol{\sigma}(t)}{d t} = \boldsymbol{\mu}
# \left(\mathbf{D} \mathbf{v}(t) + \underbrace{\mathbf{H}^{-1}\left(\tau_{21}\frac{\boldsymbol{e}_{1}}{Z_{s}}B_0\left(v_1, \sigma_1, Z_{s}, r_0\right) - \tau_{22}\frac{\boldsymbol{e}_{N}}{Z_{s}}B_L\left(v_N, \sigma_N, Z_{s}, r_n\right)\right)}_{SAT \to 0}\right).
# \end{align}
# Here $\tau_{ij}$ are penalty parameters determined by requiring stability.
#
# Approximate the mechanical energy by the above quadrature rule, having
# \begin{align}
# \mathcal{E}( t) = \sum_{j}^{N}\frac{1}{2}\left(\rho_jv_j^2 + \frac{1}{\mu_j}\sigma_j^2\right)h_j > 0.
# \end{align}
#
# By chosing the penalty parameters $\tau_{ij} = 1$, the semi-discrete approximation satisfies the energy estimate:
# \begin{align}
# \frac{d \mathcal{E}( t)}{d t} = -\frac{1}{2}\left(\left(1-r_0\right)Zv_1^2 + \frac{\left(1+r_0\right)}{Z}\sigma_1^2 +
# \left(1-r_n\right)Zv_N^2 + \frac{\left(1+r_n\right)}{Z}\sigma_N^2\right) \le 0.
# \end{align}
#
#
# 3) The discrete mechanical energy can never grow in time, and thus the semidiscrete numerical approximation is asymptotically stable.
#
# 4)  Time integration can be performed using any stable time stepping scheme. This notebook implements the fourth order accurate Runge-Kutta method.
#
# To keep the problem simple, we use as spatial initial condition a Gauss function with half-width $\delta$
#
# \begin{equation}
# v(x,t=0)  = e^{-1/\delta^2 (x - x_{o})^2}, \quad \sigma(x,t=0) = 0
# \end{equation}

# **** Exercises****
#
#

# + {"code_folding": [0]}
# Parameters initialization and plotting the simulation
# Import necessary routines
import numpy as np
import matplotlib.pyplot as plt
import time_integrator
import rate
import utils
import timeit

# plt.switch_backend("TkAgg")          # plots in external window
plt.switch_backend("nbagg")  # plots within this notebook

# + {"code_folding": [0]}
# Initializations
L = 10.0  # length of the domain (km)
t = 0.0  # initial time
tend = 1.45  # final time
nx = 501  # grid points in x
dx = L / (nx - 1)  # grid increment in x
cs = 3.464  # velocity (km/s) (can be an array)
iplot = 5  # snapshot frequency
rho = 2.6702  # density [g/cm^3]
mu = rho * cs ** 2  # shear modulus [GPa]
Zs = rho * cs  # shear impedance

order = 6  # order of accuracy

# Initialize the domain
y = np.zeros((nx, 1))

# Initial particle velocity perturbation and discretize the domain
for j in range(0, nx):
    y[j, :] = j * dx  # discrete domain

# Time stepping parameters
cfl = 1.0  # CFL number
dt = (cfl / cs) * dx  # Time step
nt = int(round(tend / dt))  # number of time steps
n = 0  # counter

# Boundary condition reflection coefficients
r0 = 1  # r=0:absorbing, r=1:free-surface, r=-1: clamped
r1 = 1  # r=0:absorbing, r=1:free-surface, r=-1: clamped

# penalty parameters
tau_11 = 1
tau_12 = 1
tau_21 = 1
tau_22 = 1

# Initialize: particle velocity (v); and shear stress (s)
v = np.zeros((nx, 1))
s = np.zeros((nx, 1))

U = np.zeros((nx, 1))
V = np.zeros((nx, 1))
U_t = np.zeros((nx, 1))
V_t = np.zeros((nx, 1))
U_x = np.zeros((nx, 1))
V_x = np.zeros((nx, 1))

# Difference between analyticla and numerical solutions
EV = [0]  # initialize errors in V (velocity)
EU = [0]  # initialize errors in U (stress)
T = [0]  # later append every time steps to this

# + {"code_folding": [0]}
# Computation and plotting

# Initialize animated plot for velocity and stress
fig1 = plt.figure(figsize=(10, 10))
ax1 = fig1.add_subplot(4, 1, 1)
line1 = ax1.plot(y, v, "r", y, U, "k--")
plt.title("numerical vs exact")
plt.xlabel("x [km]")
plt.ylabel("velocity [m/s]")

ax2 = fig1.add_subplot(4, 1, 2)
line2 = ax2.plot(y, s, "r", y, V, "k--")
plt.title("numerical vs exact")
plt.xlabel("x[km]")
plt.ylabel("stress [MPa]")

# Initialize error plot (for velocity and stress)
ax3 = fig1.add_subplot(4, 1, 3)
line3 = ax3.plot(T, EV, "r")
plt.title("relative error in particle velocity")
plt.xlabel("time [s]")
ax3.set_ylim([10 ** -5, 1])
plt.ylabel("error")

ax4 = fig1.add_subplot(4, 1, 4)
line4 = ax4.plot(T, EU, "r")
plt.ylabel("error")
plt.xlabel("time[t]")
ax4.set_ylim([10 ** -5, 1])
plt.title("relative error in stress")

plt.tight_layout()
plt.ion()
plt.show()

t = 0  # initial time

forcing = (
    1.0
)  # forcing function, forcing = 1,  and no forcing function, forcing = 0

# type of initial data: Gaussian or Sinusoidal
type_0 = "Gaussian"
# type_0 = 'Sinusoidal'

if type_0 in ("Sinusoidal"):
    forcing = 1.0  # we must use forcing for Sinusoidal initial condition

# L2-norm normalizer
# Generate  conditions for normalization
rate.mms(v, s, U_t, V_t, U_x, V_x, y, 0.65, type_0)
A = np.linalg.norm(v)
B = np.linalg.norm(s)

# Loop through time and evolve the wave-fields using ADER time-stepping scheme of N+1 order of accuracy
start = timeit.default_timer()

# Generate initial conditions
rate.mms(v, s, U_t, V_t, U_x, V_x, y, t, type_0)

for t in utils.drange(0.0, tend + dt, dt):
    n = n + 1

    # compute numerical solution
    time_integrator.elastic_RK4(
        v,
        s,
        v,
        s,
        rho,
        mu,
        nx,
        dx,
        order,
        y,
        t,
        dt,
        r0,
        r1,
        tau_11,
        tau_21,
        tau_12,
        tau_22,
        type_0,
        forcing,
    )

    # Analytical solution
    rate.mms(U, V, U_t, V_t, U_x, V_x, y, t + dt, type_0)

    # compute error and append to the error array
    EU.append(np.linalg.norm(U - v) / A)
    EV.append(np.linalg.norm(V - s) / B)

    T.append(t)

    # Updating plots
    if n % iplot == 0:
        for l in line1:
            l.remove()
            del l
        for l in line2:
            l.remove()
            del l
        for l in line3:
            l.remove()
            del l
        for l in line4:
            l.remove()
            del l

        # Display lines
        line1 = ax1.plot(y, v, "r", y, U, "k--")
        ax1.legend(iter(line1), ("Numerical", "Analytical"))
        line2 = ax2.plot(y, s, "r", y, V, "k--")
        ax2.legend(iter(line2), ("Numerical", "Analytical"))
        line3 = ax3.plot(T, EU, "k--")
        ax3.set_yscale("log")  # , nonposx='clip')
        line4 = ax4.plot(T, EV, "k--")
        ax4.set_yscale("log")  # , nonposx='clip')
        plt.gcf().canvas.draw()

plt.ioff()
plt.show()

# Simulation end time
stop = timeit.default_timer()
print(
    "total simulation time = ", stop - start
)  # print the time required for simulation
print(
    "spatial order  of accuracy = ", order
)  # print the polynomial degree used
print("number of grid points = ", nx)  # print the degree of freedom
print(
    "maximum relative error in particle velocity = ", max(EU)
)  # max. relative error in particle velocity
print(
    "maximum relative error in stress = ", max(EV)
)  # max. relative error in stress
# -


text:
import datetime
from datetime import date
from Model.UserInputs import *
import numpy as np

# black scholes variables
r = 0.01
T = 0
S = df["Adj Close"][-1]
K = []
sigma = []
current_prices = []
open_interest = []
value_type = []

class CleaningData:

    @staticmethod
    def callOrPut():
        global contract_type
        if 'c' in option_type.lower():
            contract_type = chain['calls']
        else:
            contract_type = chain['puts']

    @staticmethod
    def gettingStrike():
        global K
        vector = np.vectorize(np.float)
        K = list(vector(contract_type[['Strike']].values))
        K = [float(x) for x in K]

    @staticmethod
    def gettingLastPrice():
        global current_prices
        vector = np.vectorize(np.float)
        current_prices = list(vector(contract_type[['Last Price']].values))
        current_prices = [float(x) for x in current_prices]

    @staticmethod
    def gettingSigma():
        "this is getting implied volatility"
        global sigma
        new_sig = []
        sigma = list(contract_type[['Implied Volatility']].values)

        for i in sigma:
            iv = float(i[-1][:-2]) / 100
            new_sig.append(iv)
        sigma = new_sig

    @staticmethod
    def gettingT():
        global T
        contract_name = str(contract_type[['Contract Name']].values[0][0])
        date_str = contract_name[len(stock):]
        year = 2000 + int(date_str[:2])
        month = int(date_str[2:4])
        day = int(date_str[4:6])

        today = datetime.today().date()
        expiry_date = date(year, month, day)
        diff = expiry_date - today
        T = diff.days / 365

    @staticmethod
    def getOpenInterest():
        global open_interest
        open_interest = list(contract_type['Open Interest'])

    @staticmethod
    def getValueType():
        for strike in K:
            if strike >= S:
                value_type.append("---Out---")
            else:
                value_type.append("***In***")

    @staticmethod
    def printers():
        print(f'Strike Prices: {K}\n',
              f'Implied Volatility: {sigma}\n',
              f'Value Type: {value_type}')

CleaningData.callOrPut()
CleaningData.gettingT()
CleaningData.gettingStrike()
CleaningData.gettingSigma()
CleaningData.gettingLastPrice()
CleaningData.getValueType()

ext: py

max_stars_repo_name: megz-mugz/OptionPricing-

max_issues_repo_name: megz-mugz/OptionPricing-

max_forks_repo_name: megz-mugz/OptionPricing-

reason: import numpy


text:
# -*- coding: utf-8 -*-

# %% IMPORTS
# Package imports
import numpy as np
import pytest

# e13Tools imports
from e13tools.core import ShapeError
from e13tools.sampling import lhd

# %% PYTEST HELPER FUNCTIONS
# Create fixture that resets the NumPy random seed for every test
@pytest.fixture()
def set_numpy_random_seed():
    np.random.seed(0)

# %% PYTEST CLASSES AND FUNCTIONS
# Pytest class for lhd()-function
@pytest.mark.usefixtures('set_numpy_random_seed')
class Test_lhd(object):
    # Test 3-2 LHD with random samples
    def test_3_2_random(self):
        assert np.allclose(lhd(3, 2, method='random'),
                           np.array([[0.18293783, 0.47919574],
                                     [0.86758779, 0.96392433],
                                     [0.57172979, 0.21529804]]))

    # Test 3-2 LHD with fixed samples
    def test_3_2_fixed(self):
        assert np.allclose(lhd(3, 2, method='fixed'),
                           np.array([[1, 1], [0.5, 0], [0, 0.5]]))

    # Test 3-2 LHD with centered samples
    def test_3_2_center(self):
        assert np.allclose(lhd(3, 2, method='center'),
                           np.array([[0.83333333, 0.83333333],
                                     [0.5, 0.16666667], [0.16666667, 0.5]]))

    # Test 3-2 LHD with specified value range
    def test_3_2_custom_rng(self):
        val_rng = [[0, 2], [1, 4]]
        assert np.allclose(lhd(3, 2, val_rng=val_rng),
                           np.array([[0.36587567, 2.43758721],
                                     [1.73517558, 3.89177300],
                                     [1.14345958, 1.64589411]]))

    # Test 3-3 LHD with invalid val_rng shape raises an error
    def test_3_3_invalid_val_rng_shape(self):
        val_rng = [[0, 2], [0, 1]]
        with pytest.raises(ShapeError):
            lhd(3, 3, val_rng=val_rng)

    # Test 3-2 LHD with constraints given
    def test_3_2_constraints(self):
        constraints = np.random.rand(2, 2)
        assert np.allclose(lhd(3, 2, constraints=constraints, criterion=0),
                           np.array([[1., 1], [0.5, 0], [0, 0.5]]))

    # Test 3-2 LHD with constraints given same sample as LHD
    def test_3_2_same_constraints(self):
        constraints = [0, 0]
        assert np.allclose(lhd(3, 2, constraints=constraints, criterion=0),
                           np.array([[0, 1], [1, 0], [0.5, 0.5]]))

    # Test 3-2 LHD with constraints given outside of val_rng
    def test_3_2_constraints_outside(self):
        constraints = np.random.rand(2, 2)
        val_rng = [[1, 2], [2, 3]]
        assert np.allclose(lhd(3, 2, val_rng=val_rng,
                               constraints=constraints),
                           np.array([[1.54863137, 2.46114717],
                                     [1.14121827, 2.32122092],
                                     [1.81252907, 2.93057501]]))

    # Test 3-2 LHD with no constraints given
    def test_3_2_no_constraints(self):
        assert np.allclose(lhd(3, 2, constraints=[]),
                           np.array([[0.18293783, 0.47919574],
                                     [0.86758779, 0.96392433],
                                     [0.57172979, 0.21529804]]))

    # Test 3-2 LHD with 3D constraints raises an error
    def test_3_2_3D_constraints(self):
        constraints = np.random.rand(2, 2, 2)
        with pytest.raises(ShapeError):
            lhd(3, 2, constraints=constraints)

    # Test 3-2 LHD with invalid constraints raises an error
    def test_3_2_invalid_constraints(self):
        constraints = np.random.rand(2, 3)
        with pytest.raises(ShapeError):
            lhd(3, 2, constraints=constraints)

    # Test 3-2 LHD with maximin criterion, string
    def test_3_2_maximin_str(self):
        assert np.allclose(lhd(3, 2, method='random', criterion='maximin'),
                           np.array([[1, 1], [0.5, 0], [0, 0.5]]))

    # Test 3-2 LHD with maximin criterion, float
    def test_3_2_maximin_float(self):
        assert np.allclose(lhd(3, 2, method='fixed', criterion=0),
                           np.array([[1, 1], [0.5, 0], [0, 0.5]]))

    # Test 3-2 LHD with correlation criterion, string
    def test_3_2_correlation_str(self):
        assert np.allclose(lhd(3, 2, method='center', criterion='correlation'),
                           np.array([[0.83333333, 0.83333333],
                                     [0.5, 0.16666667], [0.16666667, 0.5]]))

    # Test 3-2 LHD with correlation criterion, float
    def test_3_2_correlation_float(self):
        assert np.allclose(lhd(3, 2, method='center', criterion=1),
                           np.array([[0.83333333, 0.83333333],
                                     [0.5, 0.16666667], [0.16666667, 0.5]]))

    # Test 3-2 LHD with multi criterion, string
    def test_3_2_multi_str(self):
        assert np.allclose(lhd(3, 2, method='center', criterion='multi'),
                           np.array([[0.83333333, 0.83333333],
                                     [0.5, 0.16666667], [0.16666667, 0.5]]))

    # Test 3-2 LHD with multi criterion, float
    def test_3_2_multi_float(self):
        assert np.allclose(lhd(3, 2, method='fixed', criterion=0.5),
                           np.array([[1, 1], [0.5, 0], [0, 0.5]]))

    # Test 1-2 LHD with fixed samples
    def test_1_2_fixed(self):
        assert np.allclose(lhd(1, 2, method='fixed', criterion=0),
                           np.array([[0.5, 0.5]]))

    # Test 2-2 LHD with fixed samples
    def test_2_2_fixed(self):
        assert np.allclose(lhd(2, 2, method='fixed', criterion=0),
                           np.array([[1, 0], [0, 1]]))

    # Test 3-1 LHD with fixed samples
    def test_3_1_fixed(self):
        assert np.allclose(lhd(3, 1, method='fixed', criterion=0),
                           np.array([[1], [0.5], [0]]))

    # Test 3-2 LHD with invalid criterion raises an error, string
    def test_3_2_invalid_criterion_str(self):
        with pytest.raises(ValueError):
            lhd(3, 2, criterion='test')

    # Test 3-2 LHD with invalid criterion raises an error, float
    def test_3_2_invalid_criterion_float(self):
        with pytest.raises(ValueError):
            lhd(3, 2, criterion=1.5)

    # Test 3-2 LHD and request the score
    def test_3_2_score(self):
        results = lhd(3, 2, criterion=0, get_score=True)
        assert np.allclose(results[0], np.array([[1, 1], [0.5, 0], [0, 0.5]]))
        assert np.allclose(results[1], np.array([0.80444958, 0.5, 0.80444958]))

    # Test 4-2 LHD and do not request quick-scan
    def test_4_2_no_quick_scan(self):
        assert np.allclose(lhd(4, 2, criterion=0, iterations=100,
                               quickscan=False),
                           np.array([[0.33333333, 0], [0, 0.66666667],
                                     [1, 0.33333333], [0.66666667, 1]]))


text:
import json
import networkx
import networkx as nx
from collections import Counter

G = nx.Graph()

file = ''

with open(file, 'r', encoding='utf8') as jsonfile:
     data = json.load(jsonfile)
     messages = data['messages']
     print('File', file, 'loaded.')

def FindAllUsers():
        users = set()
        for node in messages:
            users.add(node.get('from'))
            print('User', node.get('from'), 'discovered!')
        users.remove(None)
        return users

def FindUserNeighbours(user):
    ids = set()
    neighbours = set()
    for i in messages:
        if i.get('from') == user:
           ids.add(i.get('reply_to_message_id'))

    for i in messages:
        if i.get('id') in ids:
            if i.get('from') in neighbours:
                pass
            else:
                neighbours.add(i.get('from'))
                print('Neighbour',i.get('from'),'for', user, 'discovered!')
        else:
            pass

    return neighbours

def CreateNeighboursDict():
    for i in FindAllUsers():
        if i is None:
            pass
        else:
            for r in FindUserNeighbours(i):
                if r is None:
                    pass
                else:
                    if r == i:
                        pass
                    else:
                        w = CountUserRepliesTo(i, r)
                        print('Adding edges:', i, r, w)
                        G.add_edge(i, r, weight=w)
                        # G.add_edge(i, r)

def CountUserMessages(user):
    counter = set()
    for i in messages:
        if i.get('from') == user:
            if i.get('from') is None:
                pass
            else:
                counter.add(i.get('from'))
    return counter

def CountAllUsersMessages():
    counted = set()
    for i in FindAllUsers():
        counted.add(CountUserMessages(i))
    return counted

def CountUserRepliesTo(user, ruser):
    ids = set()
    counted = []
    for i in messages:
        if i.get('from') == user:
            if i.get('reply_to_message_id') in ids:
                pass
            else:
                ids.add(i.get('reply_to_message_id'))

    for i in messages:
        if i.get('id') in ids:
                counted.append(i.get('from'))
        else:
            pass
    out = Counter(counted)
    ruser = out.get(ruser)
    return ruser

CreateNeighboursDict()
# remove_nodes = ['user1', 'user2']
# G.remove_nodes_from(remove_nodes)
# nx.draw_networkx(G, node_size=300, node_color='#73a7fa', edge_color='#e8e8e8', with_labels=True, font_size=5)
print('------------------ FINISHED ------------------')
networkx.write_gexf(G, 'result.gexf')
print('Exported')

reason: import networkx