In [1]:

    

%load_ext fortranmagic


    

In [2]:

    

%%fortran
SUBROUTINE fortran_lstsqr(ary_x, ary_y, slope, y_interc)
    ! Computes the least-squares solution to a linear matrix equation. """
    IMPLICIT NONE
    REAL(8), INTENT(in), DIMENSION(:) :: ary_x, ary_y
    REAL(8), INTENT(out) :: slope, y_interc
    REAL(8) :: x_avg, y_avg, var_x, cov_xy, temp
    INTEGER(8) :: N, i
    
    N = SIZE(ary_x)

    x_avg = SUM(ary_x) / N
    y_avg = SUM(ary_y) / N
    var_x = 0
    cov_xy = 0
    
    DO i = 1, N
        temp = ary_x(i) - x_avg
        var_x = var_x + temp**2
        cov_xy = cov_xy + (temp*(ary_y(i) - y_avg))
    END DO
    
    slope = cov_xy / var_x
    y_interc = y_avg - slope*x_avg

END SUBROUTINE fortran_lstsqr


    

In [3]:

    

%%fortran
subroutine linregress(x, y, slope, intercept)
implicit none
! Calculates simple linear regression of (x, y)
real(8), intent(in) :: x(:), y(:)
real(8), intent(out) :: slope ! y = intercept + slope*x
real(8), intent(out) :: intercept ! y = intercept + slope*x
real(8) :: xmean, ymean, varx, covxy
integer :: N
N = size(x)
xmean = sum(x)/N
ymean = sum(y)/N
varx  = dot_product(x, x) - N*xmean**2    ! = dot_product(x-xmean, x-xmean)
covxy = dot_product(x, y) - N*xmean*ymean ! = dot_product(x-xmean, y-ymean)
slope = covxy / varx
intercept = ymean - slope*xmean
end subroutine


    

In [4]:

    

import platform
import multiprocessing
from cython import __version__ as cython__version__

def print_sysinfo():
    
    print('\nPython version  :', platform.python_version())
    print('compiler        :', platform.python_compiler())
    print('Cython version  :', cython__version__)
    print('NumPy version   :', np.__version__)
    
    print('\nsystem     :', platform.system())
    print('release    :', platform.release())
    print('machine    :', platform.machine())
    print('processor  :', platform.processor())
    print('CPU count  :', multiprocessing.cpu_count())
    print('interpreter:', platform.architecture()[0])
    print('\n\n')


    

In [7]:

    

import numpy as np

print_sysinfo()

x = ([x_i*np.random.randint(8,12)/10 for x_i in range(10000)])
y = ([y_i*np.random.randint(10,14)/10 for y_i in range(10000)])
xary = np.asarray(x)
yary = np.asarray(y)


%timeit fortran_lstsqr(xary, yary)
%timeit linregress(xary, yary)


    

    




Python version  : 3.4.1
compiler        : GCC 4.2.1 (Apple Inc. build 5577)
Cython version  : 0.20.1
NumPy version   : 1.8.1

system     : Darwin
release    : 13.2.0
machine    : x86_64
processor  : i386
CPU count  : 2
interpreter: 64bit



10000 loops, best of 3: 44.4 µs per loop
10000 loops, best of 3: 52.7 µs per loop

In [5]:

    

import numpy as np

print_sysinfo()

x = ([x_i*np.random.randint(8,12)/10 for x_i in range(10000)])
y = ([y_i*np.random.randint(10,14)/10 for y_i in range(10000)])
xary = np.asarray(x)
yary = np.asarray(y)


%timeit fortran_lstsqr(xary, yary)
%timeit linregress(xary, yary)


    

    




('\nPython version  :', '2.7.7')
('compiler        :', 'GCC 4.0.1 (Apple Inc. build 5493)')
('Cython version  :', '0.20.1')
('NumPy version   :', '1.7.1')
('\nsystem     :', 'Darwin')
('release    :', '13.2.0')
('machine    :', 'x86_64')
('processor  :', 'i386')
('CPU count  :', 2)
('interpreter:', '64bit')



10000 loops, best of 3: 74.3 µs per loop
10000 loops, best of 3: 81 µs per loop

In [ ]: