from __future__ import division

import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt

import sys
print(sys.prefix)

!echo $PYTHONPATH

%%file circle.f
PROGRAM CIRCLE
REAL R, AREA
! This program reads a real number r and prints
! the area of a circle with radius r.
WRITE (*,*) 'Enter radius, r = '
READ (*,*) R
AREA = 3.14159*R*R
WRITE (*,*) 'Area, A = ', AREA
STOP
END

%%file circle.f
program circle
real r, area, pi
parameter (pi = 3.14159)
! This program reads a real number r and prints
! the area of a circle with radius r.
write (*,*) 'Enter radius, r = '
read (*,*) R
area = 3.14159*r*r
write (*,*) 'Area, A = ', area
stop
end

%%file hellofortran.f
! File  hellofortran.f
        subroutine hellofortran (n)
        integer n
       
        do 100 i=0, n
            print *, "Fortran says hello"
100     continue
        end

!f2py -c -m hellofortran hellofortran.f

%%file hello.py
import hellofortran

hellofortran.hellofortran(5)

# run the script
!python hello.py

%%file dprod.f

       subroutine dprod(x, y, n)
    
       double precision x(n), y
       y = 1.0
    
       do 100 i=1, n
           y = y * x(i)
100    continue
       end


!rm -f dprod.pyf
!f2py -m dprod -h dprod.pyf dprod.f

!cat dprod.pyf

%%file dprod.pyf
python module dprod ! in 
    interface  ! in :dprod
        subroutine dprod(x,y,n) ! in :dprod:dprod.f
            double precision dimension(n), intent(in) :: x
            double precision, intent(out) :: y
            integer, optional,check(len(x)>=n),depend(x),intent(in) :: n=len(x)
        end subroutine dprod
    end interface 
end python module dprod

!f2py -c dprod.pyf dprod.f  #Compile the code which can be included in Python

import dprod #Using the module in Python

help (dprod)

dprod.dprod(np.arange(1,50))

np.prod(np.arange(1.0,50.0))

xvec = np.random.rand(1000)

timeit dprod.dprod(xvec)

timeit xvec.prod()

# pure Python 
def py_dcumsum(a):
    b = np.empty_like(a)
    b[0] = a[0]
    for n in range(1,len(a)):
        b[n] = b[n-1]+a[n]
    return b

%%file dcumsum.f
c File dcumsum.f
       subroutine dcumsum(a, b, n)
       double precision a(n)
       double precision b(n)
       integer n
cf2py  intent(in) :: a
cf2py  intent(out) :: b
cf2py  intent(hide) :: n

       b(1) = a(1)
       do 100 i=2, n
           b(i) = b(i-1) + a(i)
100    continue
       end

!f2py -c dcumsum.f -m dcumsum

import dcumsum

a = np.array([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0])

py_dcumsum(a)

dcumsum.dcumsum(a)

a = np.random.rand(10000)

timeit py_dcumsum(a)

timeit dcumsum.dcumsum(a)

%%file test.c
#include <stdio.h>
int main() {
    printf( "Happy Halloween! \n" );
    getchar();
    return 0;
}

!gcc -o test test.c

!./test

import ctypes
ctypes.cdll.LoadLibrary("libc.so.6")
libc = ctypes.CDLL("libc.so.6")
libc.printf

print(libc.time(None)) 

%%file functions.c

#include <stdio.h>

void hello(int n);

double dprod(double *x, int n);

void dcumsum(double *a, double *b, int n);

void
hello(int n)
{
    int i;
    
    for (i = 0; i < n; i++)
    {
        printf("C says hello\n");
    }
}


double 
dprod(double *x, int n)
{
    int i;
    double y = 1.0;
    
    for (i = 0; i < n; i++)
    {
        y *= x[i];
    }

    return y;
}

void
dcumsum(double *a, double *b, int n)
{
    int i;
    
    b[0] = a[0];
    for (i = 1; i < n; i++)
    {
        b[i] = a[i] + b[i-1];
    }
}

!gcc -c -Wall -O2 -Wall -ansi -pedantic -fPIC -o functions.o functions.c
!gcc -o libfunctions.so -shared functions.o

!file libfunctions.so

%%file functions.py

import numpy
import ctypes

_libfunctions = numpy.ctypeslib.load_library('libfunctions', '.')

_libfunctions.hello.argtypes = [ctypes.c_int]
_libfunctions.hello.restype  =  ctypes.c_void_p

_libfunctions.dprod.argtypes = [numpy.ctypeslib.ndpointer(dtype=numpy.float), ctypes.c_int]
_libfunctions.dprod.restype  = ctypes.c_double

_libfunctions.dcumsum.argtypes = [numpy.ctypeslib.ndpointer(dtype=numpy.float), numpy.ctypeslib.ndpointer(dtype=numpy.float), ctypes.c_int]
_libfunctions.dcumsum.restype  = ctypes.c_void_p

def hello(n):
    return _libfunctions.hello(int(n))

def dprod(x, n=None):
    if n is None:
        n = len(x)
    x = numpy.asarray(x, dtype=numpy.float)
    return _libfunctions.dprod(x, int(n))

def dcumsum(a, n):
    a = numpy.asarray(a, dtype=numpy.float)
    b = numpy.empty(len(a), dtype=numpy.float)
    _libfunctions.dcumsum(a, b, int(n))
    return b

%%file run_hello_c.py

import functions

functions.hello(3)

!python run_hello_c.py

import functions

functions.dprod([1,2,3,4,5]) #product

a = np.random.rand(100000)

res_c = functions.dcumsum(a, len(a)) 

res_fortran = dcumsum.dcumsum(a)

res_c - res_fortran

timeit functions.dcumsum(a, len(a))

timeit dcumsum.dcumsum(a)

timeit a.cumsum()

from numpy import zeros
from scipy import weave

dx = 0.1
dy = 0.1
dx2 = dx*dx
dy2 = dy*dy

def py_update(u):
    nx, ny = u.shape
    for i in xrange(1,nx-1):
        for j in xrange(1, ny-1):
            u[i,j] = ((u[i+1, j] + u[i-1, j]) * dy2 +
                      (u[i, j+1] + u[i, j-1]) * dx2) / (2*(dx2+dy2))

def calc(N, Niter=100, func=py_update, args=()):
    u = zeros([N, N])
    u[0] = 1
    for i in range(Niter):
        func(u,*args)
    return u

%timeit calc(100, 200, func=py_update)

def num_update(u):
    u[1:-1,1:-1] = ((u[2:,1:-1]+u[:-2,1:-1])*dy2 + 
                    (u[1:-1,2:] + u[1:-1,:-2])*dx2) / (2*(dx2+dy2))

%timeit calc(100, 200, func=num_update)

def weave_update(u):
    code = """
    int i, j;
    for (i=1; i<Nu[0]-1; i++) {
       for (j=1; j<Nu[1]-1; j++) {
           U2(i,j) = ((U2(i+1, j) + U2(i-1, j))*dy2 + \
                       (U2(i, j+1) + U2(i, j-1))*dx2) / (2*(dx2+dy2));
       }
    }
    """
    weave.inline(code, ['u', 'dx2', 'dy2'])

%timeit calc(100, 200, func=weave_update)

%%file factorial.pyx
def factorial(n):
    return 1 if (n < 1) else n * factorial(n-1)

!cython factorial.pyx

#!gcc -shared -fPIC -o factorial.so factorial.c
!gcc -I /software/canopy-1.4.1/installer/appdata/canopy-1.4.1.1975.rh5-x86_64/include/python2.7/ -shared -fPIC -o factorial.so factorial.c

import factorial
print(factorial.factorial(6))

from distutils.core import setup
from Cython.Build import cythonize

setup(
    ext_modules=cythonize("factorial.pyx"),
)
# problem is clang ñ gcc here (-fopenmp)

%tb

!python setup.py build_ext --inplace

%load_ext cythonmagic

%%cython
cimport numpy

def cy_dcumsum2(numpy.ndarray[numpy.float64_t, ndim=1] a, numpy.ndarray[numpy.float64_t, ndim=1] b):
    cdef int i, n = len(a)
    b[0] = a[0]
    for i from 1 <= i < n:
        b[i] = b[i-1] + a[i]
    return b

import numpy as np

x = np.ones((32))
y = np.zeros(x.shape)

cy_dcumsum2(x,y)

print(y)