!ls /usr/lib/

from ctypes import CDLL

libc_name = 'libc.dylib'  # OSX
# libc_name = 'libc.so.6'  # Linux
# libc_name = 'libc.dll'  # Windows

libc = CDLL(libc_name)

libc.time

print "Seconds since January 1, 1970:"
print libc.time()

%%file my_sum.c
#include <stdio.h>

// sum all the values in the array x
// x is a pointer to a memory block 
// of length n
int sum(int *x, int n)
{
  int i, counter;
  counter = 0;
  for(i=0; i<n; i++)
    {
      counter += x[i];
    }
  return counter;
}

%%bash
gcc -c my_sum.c
gcc -shared -o my_sum.so my_sum.o

!ls my_sum.*

from ctypes import CDLL, c_void_p
import numpy as np

my_sum=CDLL('my_sum.so')

a = np.arange(10, dtype=np.int32)
adata = a.ctypes.data_as(c_void_p)
asize = a.size

print my_sum.sum(adata, asize)
print a.sum()

%%file my_sum2.c
#include <stdio.h>

// Define a simple array struct, containing
// a pointer and a length
struct Array{
   int *x;
   int n;
};

// Sum the values in the array struct
int sum(struct Array a)
{
  int counter, i;
  counter = 0;
  for(i=0; i<a.n; i++)
    {
      counter += a.x[i];
    }
  return counter;
}
 

%%bash
gcc -c my_sum2.c
gcc -shared -o my_sum2.so my_sum2.o

from ctypes import Structure, POINTER, c_int

# Here's our C structure defined in Python
class Array(Structure):
    _fields_ = ("x", c_void_p), ("n", c_int)
    
my_sum2=CDLL('my_sum2.so')
    
a = np.arange(10, dtype=np.int32)
arr = Array(a.ctypes.data_as(c_void_p), a.size)

print my_sum2.sum(arr)
print a.sum()

from ctypes import CDLL, c_int, c_void_p, byref
import numpy as np

def solve(M, b):
    """Solve the linear equation M x = b using LAPACK's DGESV"""
    # Interface to the LAPACK system library
    # This may be different depending on your operating system
    lapack = CDLL('liblapack.dylib')
    
    # Make sure M is double precision fortran-ordered
    M = np.asarray(M, dtype=np.float64, order='F')
    
    # Make a copy of b.  It will be overwritten with the result
    b = np.array(b, dtype=np.float64, copy=True)
    n = b.size
    
    # We could generalize this, but we'll do a simple problem for now
    assert M.shape == (n, n)
    assert b.shape == (n,)
    
    # prepare the variables for dgesv.
    size = c_int(n)
    one = c_int(1)
    info = c_int(0)
    lapack.dgesv(byref(size), # size of the problem
                 byref(one), # number of columns of b
                 M.ctypes.data_as(c_void_p), # M data array
                 byref(size), # size of the problem
                 (n * c_int)(), # integer work array
                 b.ctypes.data_as(c_void_p), # b data array
                 byref(size), # size of b
                 byref(info)) # exit flag
    return b

X = np.random.random((5, 5))
b = np.random.random(5)

print solve(X, b)

print np.linalg.solve(X, b)

import numpy as np

def fib_py(N):
    x = np.zeros(N, dtype=float)
    for i in range(N):
        if i == 0:
            x[i] = 0
        elif i == 1:
            x[i] = 1
        else:
            x[i] = x[i - 1] + x[i - 2]
    return x

print fib_py(10)

%timeit fib_py(10000)

%%file fib.f

      subroutine fib(A,N)
      integer N
      double precision A(N)
      do i=1,N
         if (i.EQ.1) then
            A(i) = 0.0D0
         elseif (i.EQ.2) then
            A(i) = 1.0D0
         else
            A(i) = A(i-1) + A(i-2)
         endif
      enddo
      end

!f2py -c -m fib fib.f

import numpy as np
import fib
a = np.zeros(10, dtype='d')
fib.fib(a) 
print a 

a = np.zeros(10000, dtype='d')
%timeit fib.fib(a)

!rm -f _fib2.pyf
!f2py fib.f -m fib2 -h _fib2.pyf

!cat _fib2.pyf

from fib import fib
a = np.zeros(10)
fib(a, 10)  # specify the optional keyword
print a

a = np.zeros(10)
fib(a, 5)  # specify only part of the array
print a

a = np.zeros(10)
fib(a, 20)  # specify a length which is too long

%%file fib2.pyf
!    -*- f90 -*-
! Note: the context of this file is case sensitive.

python module fib2 ! in 
    interface  ! in :fib2
        subroutine fib(a,n) ! in :fib2:fib.f
            double precision dimension(n), intent(out), depend(n) :: a
            integer intent(in) :: n
        end subroutine fib
    end interface 
end python module fib2

!f2py -c -m fib2 fib2.pyf fib.f

import fib2
print fib2.fib(10)

fib2.fib?

%timeit fib2.fib(10000)

def nth_fib(n):
    a, b = 0, 1
    for i in range(n):
        b, a = a + b, b
    return a

[nth_fib(i) for i in range(10)]

%timeit nth_fib(10000)

%load_ext cythonmagic

%%cython
def nth_fib2(n):
    a, b = 0, 1
    for i in range(n):
        b, a = a + b, b
    return a

[nth_fib2(i) for i in range(10)]

%timeit nth_fib2(10000)

%%cython
def nth_fib3(int n):
    cdef int a = 0
    cdef int b = 1
    cdef int tmp
    for i in range(n):
        tmp = b
        b = a + b
        a = tmp
    return a

[nth_fib3(i) for i in range(10)]

print "Python only:"
%timeit nth_fib(10000)
print "Bare Cython:"
%timeit nth_fib(10000)
print "Typed Cython:"
%timeit nth_fib3(10000)

%%file nth_fib.pyx
# Slow version:
def nth_fib2(n):
    a, b = 0, 1
    for i in range(n):
        b, a = a + b, b
    return a

# Fast Version:
def nth_fib3(int n):
    cdef int a = 0
    cdef int b = 1
    cdef int tmp
    for i in range(n):
        tmp = b
        b = a + b
        a = tmp
    return a

!cython -a nth_fib.pyx

%%cython
import numpy as np

def fib_cy(N):
    x = np.zeros(N, dtype=float)
    for i in range(N):
        if i == 0:
            x[i] = 0
        elif i == 1:
            x[i] = 1
        else:
            x[i] = x[i - 1] + x[i - 2]
    return x

print np.allclose(fib_py(10000), fib_cy(10000))
%timeit fib_py(10000)
%timeit fib_cy(10000)

%%cython
import numpy as np
from numpy cimport float_t

def fib_cy2(int N):
    cdef int i
    cdef float_t[::1] x = np.zeros(N, dtype=float)
    for i in range(N):
        if i == 0:
            x[i] = 0
        elif i == 1:
            x[i] = 1
        else:
            x[i] = x[i - 1] + x[i - 2]
    return x

%timeit fib_cy(10000)
%timeit fib_cy2(10000)

def fib(N):
    x = np.zeros(N, dtype=np.float64)
    for i in range(N):
        if i == 0:
            x[i] = 0
        elif i == 1:
            x[i] = 1
        else:
            x[i] = x[i - 1] + x[i - 2]
    return x

from numba import autojit
fib_nb = autojit(fib)

fib_nb(10)

%timeit fib(10000)
%timeit fib_nb(10000)

@autojit
def fib_nb(N):
    x = np.zeros(N, dtype=np.float64)
    for i in range(N):
        if i == 0:
            x[i] = 0
        elif i == 1:
            x[i] = 1
        else:
            x[i] = x[i - 1] + x[i - 2]
    return x

%timeit fib_nb(10000)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def mandel(x, y, max_iter):
    """
    Given z = x + iy and max_iter, determine whether the candidate
    is in the mandelbrot set for the given number of iterations
    """
    c = complex(x, y)
    z = 0.0j
    for i in range(max_iter):
        z = z*z + c
        if (z.real*z.real + z.imag*z.imag) >= 4:
            return i

    return max_iter

def create_fractal(Nx, Ny, xmin, xmax, ymin, ymax, max_iter):
    """Create and return a fractal image"""
    image = np.zeros((Ny, Nx), dtype=float)
    dx = (xmax - xmin) * 1. / Nx
    dy = (ymax - ymin) * 1. / Ny
    
    for x in range(Nx):
        rpart = xmin + x * dx
        for y in range(Ny):
            ipart = ymin + y * dy
            color = mandel(rpart, ipart, max_iter)
            image[y, x] = color
    return image

image = create_fractal(300, 200, -2, 1, -1, 1, 20)
plt.imshow(image, cmap=plt.cm.jet)