import numpy as np
dt = 0.01
t = np.arange(0, 30, dt)
nse1 = np.random.randn(len(t))                 # white noise 1
nse2 = np.random.randn(len(t))                 # white noise 2
r = np.exp(-t/0.05)

cnse1 = np.convolve(nse1, r, mode='same')*dt   # colored noise 1
cnse2 = np.convolve(nse2, r, mode='same')*dt   # colored noise 2

# two signals with a coherent part and a random part
s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1
s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2
plt.subplot(211)
plt.xlim(0,5)
plt.plot(t, s1, 'b-', t, s2, 'g-')
plt.xlabel('time')
plt.ylabel('s1 and s2')
plt.grid(True)

plt.subplot(212)
cxy, f = plt.csd(s1, s2, 256, 1./dt)
plt.ylabel('CSD (db)')
plt.show()

from bokeh.plotting import output_notebook
output_notebook()

from bokeh.plotting import line, show, figure, xaxis, hold, yaxis
from bokeh.objects import Range1d
import numpy as np
dt = 0.01
t = np.arange(0, 30, dt)
nse1 = np.random.randn(len(t))                 # white noise 1
nse2 = np.random.randn(len(t))                 # white noise 2
r = np.exp(-t/0.05)

cnse1 = np.convolve(nse1, r, mode='same')*dt   # colored noise 1
cnse2 = np.convolve(nse2, r, mode='same')*dt   # colored noise 2

# two signals with a coherent part and a random part
s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1
s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2
hold()
figure(tools="pan,wheel_zoom,box_zoom,reset,previewsave", title="Signals", plot_height=400)
line(t[:500], s1[:500], color="blue", legend='s1', x_range = Range1d(start=0, end=5))
line(t[:500], s2[:500], color='green', legend='s2', x_range = Range1d(start=0, end=5))
yaxis().axis_label='s1 and s2'
yaxis().axis_label_text_font_size = "12pt"
hold(False)
show()
hold()
figure(tools="pan,wheel_zoom,box_zoom,reset,previewsave", title="CSD", plot_height=200)
csd = matplotlib.mlab.csd(s1, s2, 256, 1./dt)
line(csd[1], np.log10(np.abs(csd[0]))*10, x_range = Range1d(start=0, end=50))
yaxis().axis_label='CSD(db)'
xaxis().axis_label='Frequency'
xaxis().axis_label_text_font_size = "12pt"
yaxis().axis_label_text_font_size = "12pt"
show()

import numpy as np

def mandel(x, y, max_iters):
    """
    Given the real and imaginary parts of a complex number,
    determine if it is a candidate for membership in the Mandelbrot
    set given a fixed number of iterations.
    """
    i = 0
    c = complex(x,y)
    z = 0.0j
    for i in range(max_iters):
        z = z*z + c
        if (z.real*z.real + z.imag*z.imag) >= 4:
            return i

    return 255

def create_fractal(min_x, max_x, min_y, max_y, image, iters):
    height = image.shape[0]
    width = image.shape[1]

    pixel_size_x = (max_x - min_x) / width
    pixel_size_y = (max_y - min_y) / height
    for x in range(width):
        real = min_x + x * pixel_size_x
        for y in range(height):
            imag = min_y + y * pixel_size_y
            color = mandel(real, imag, iters)
            image[y, x] = color

    return image

image = np.zeros((1000, 1300), dtype=np.uint8)
%timeit create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
ret = create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
imshow(ret)

import numba
mandel = numba.jit(mandel)
create_fractal = numba.jit(create_fractal)

image = np.zeros((1000, 1300), dtype=np.uint8)
%timeit create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
ret = create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)

@numba.jit("f8[:,:](i8,i8,i8,i8,f8[:,:],i8)")
def create_fractal2(min_x, max_x, min_y, max_y, image, iters):
    height = image.shape[0]
    width = image.shape[1]

    pixel_size_x = (max_x - min_x) / width
    pixel_size_y = (max_y - min_y) / height
    for x in range(width):
        real = min_x + x * pixel_size_x
        for y in range(height):
            imag = min_y + y * pixel_size_y
            color = mandel(real, imag, iters)
            image[y, x] = color

    return image


create_fractal2.inspect_types()

%load_ext cythonmagic

%%cython_inline
cdef int a = 10
cdef int b = 5
return a*b

%%cython_pyximport cython_mandelbrot

cdef unsigned char mandel(double x, double y, int max_iters):
    """
    Given the real and imaginary parts of a complex number,
    determine if it is a candidate for membership in the Mandelbrot
    set given a fixed number of iterations.
    """
    cdef unsigned char i = 0
    cdef double complex c = complex(x,y)
    cdef double complex z = 0.0j
    cdef double zreal, zimag
    for i in range(max_iters):
        z = z*z + c
        if (z.real*z.real + z.imag*z.imag) >= 4:
            return i
    return 255

def create_fractal(double min_x, double max_x, double min_y, double max_y, unsigned char[:, ::1] image, int iters):
    cdef int height = image.shape[0]
    cdef int width = image.shape[1]

    cdef double pixel_size_x = (max_x - min_x) / width
    cdef double pixel_size_y = (max_y - min_y) / height
    cdef int x, y, 
    cdef unsigned char color
    cdef double real, imag
    for x in range(width):
        real = min_x + x * pixel_size_x
        for y in range(height):
            imag = min_y + y * pixel_size_y
            color = mandel(real, imag, iters)
            image[y, x] = color
    return image

from cython_mandelbrot import create_fractal
image = np.zeros((1000, 1300), dtype=np.uint8)
%timeit create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
ret = create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
imshow(ret)