Skip to content

Instantly share code, notes, and snippets.

@nbren12

nbren12/Benchmark AR(1) process.ipynb

Last active August 29, 2015 14:06
Show Gist options
  • Save nbren12/d8fec9a9cdeb917258d0 to your computer and use it in GitHub Desktop.
Save nbren12/d8fec9a9cdeb917258d0 to your computer and use it in GitHub Desktop.
Download ZIP
Benchmark AR(1) process
Display the source blob
Display the rendered blob
Raw
Benchmark AR(1) process.ipynb
{
"metadata": {
"name": "",
"signature": "sha256:4b17418321a59b5ae5d24d32c5a05e2825723d3f39e79ecca9dbe1f833dc51f2"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, I wanted to benchmark the performance of a simple dynamical system using julia, python, and fortran. The system is a simple complex-valued AR(1) process described by \n",
"$$ u_{n+1} = F u_{n} + \\sigma_{n+1} $$ \n",
"where $\\sigma_{n+1}$ is iid complex normal noise, and F is a complex number. I find that fortran is 10 times faster than julia, which is 10 times faster than python."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Write and compile the fortran OU code"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get a random number generator from netlib"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%file ziggurat.f90\n",
"\n",
"! Marsaglia & Tsang generator for random normals & random exponentials.\n",
"! Translated from C by Alan Miller (amiller@bigpond.net.au)\n",
"\n",
"! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating\n",
"! random variables', J. Statist. Software, v5(8).\n",
"\n",
"! This is an electronic journal which can be downloaded from:\n",
"! http://www.jstatsoft.org/v05/i08\n",
"\n",
"! N.B. It is assumed that all integers are 32-bit.\n",
"! N.B. The value of M2 has been halved to compensate for the lack of\n",
"! unsigned integers in Fortran.\n",
"\n",
"! Latest version - 1 January 2001\n",
"MODULE Ziggurat\n",
" IMPLICIT NONE\n",
"\n",
" PRIVATE\n",
"\n",
" INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND( 12, 60 )\n",
" REAL(DP), PARAMETER :: m1=2147483648.0_DP, m2=2147483648.0_DP, &\n",
" half=0.5_DP\n",
" REAL(DP) :: dn=3.442619855899_DP, tn=3.442619855899_DP, &\n",
" vn=0.00991256303526217_DP, &\n",
" q, de=7.697117470131487_DP, &\n",
" te=7.697117470131487_DP, &\n",
" ve=0.003949659822581572_DP\n",
" INTEGER, SAVE :: iz, jz, jsr=123456789, kn(0:127), &\n",
" ke(0:255), hz\n",
" REAL(DP), SAVE :: wn(0:127), fn(0:127), we(0:255), fe(0:255)\n",
" LOGICAL, SAVE :: initialized=.FALSE.\n",
"\n",
" PUBLIC :: zigset, shr3, uni, rnor, rexp\n",
"\n",
"\n",
"CONTAINS\n",
"\n",
"\n",
"SUBROUTINE zigset( jsrseed )\n",
"\n",
" INTEGER, INTENT(IN) :: jsrseed\n",
"\n",
" INTEGER :: i\n",
"\n",
" ! Set the seed\n",
" jsr = jsrseed\n",
"\n",
" ! Tables for RNOR\n",
" q = vn*EXP(half*dn*dn)\n",
" kn(0) = (dn/q)*m1\n",
" kn(1) = 0\n",
" wn(0) = q/m1\n",
" wn(127) = dn/m1\n",
" fn(0) = 1.0_DP\n",
" fn(127) = EXP( -half*dn*dn )\n",
" DO i = 126, 1, -1\n",
" dn = SQRT( -2.0_DP * LOG( vn/dn + EXP( -half*dn*dn ) ) )\n",
" kn(i+1) = (dn/tn)*m1\n",
" tn = dn\n",
" fn(i) = EXP(-half*dn*dn)\n",
" wn(i) = dn/m1\n",
" END DO\n",
"\n",
" ! Tables for REXP\n",
" q = ve*EXP( de )\n",
" ke(0) = (de/q)*m2\n",
" ke(1) = 0\n",
" we(0) = q/m2\n",
" we(255) = de/m2\n",
" fe(0) = 1.0_DP\n",
" fe(255) = EXP( -de )\n",
" DO i = 254, 1, -1\n",
" de = -LOG( ve/de + EXP( -de ) )\n",
" ke(i+1) = m2 * (de/te)\n",
" te = de\n",
" fe(i) = EXP( -de )\n",
" we(i) = de/m2\n",
" END DO\n",
" initialized = .TRUE.\n",
" RETURN\n",
"END SUBROUTINE zigset\n",
"\n",
"\n",
"\n",
"! Generate random 32-bit integers\n",
"FUNCTION shr3( ) RESULT( ival )\n",
" INTEGER :: ival\n",
"\n",
" jz = jsr\n",
" jsr = IEOR( jsr, ISHFT( jsr, 13 ) )\n",
" jsr = IEOR( jsr, ISHFT( jsr, -17 ) )\n",
" jsr = IEOR( jsr, ISHFT( jsr, 5 ) )\n",
" ival = jz + jsr\n",
" RETURN\n",
"END FUNCTION shr3\n",
"\n",
"\n",
"\n",
"! Generate uniformly distributed random numbers\n",
"FUNCTION uni( ) RESULT( fn_val )\n",
" REAL(DP) :: fn_val\n",
"\n",
" fn_val = half + 0.2328306e-9_DP * shr3( )\n",
" RETURN\n",
"END FUNCTION uni\n",
"\n",
"\n",
"\n",
"! Generate random normals\n",
"FUNCTION rnor( ) RESULT( fn_val )\n",
" REAL(DP) :: fn_val\n",
"\n",
" REAL(DP), PARAMETER :: r = 3.442620_DP\n",
" REAL(DP) :: x, y\n",
"\n",
" IF( .NOT. initialized ) CALL zigset( jsr )\n",
" hz = shr3( )\n",
" iz = IAND( hz, 127 )\n",
" IF( ABS( hz ) < kn(iz) ) THEN\n",
" fn_val = hz * wn(iz)\n",
" ELSE\n",
" DO\n",
" IF( iz == 0 ) THEN\n",
" DO\n",
" x = -0.2904764_DP* LOG( uni( ) )\n",
" y = -LOG( uni( ) )\n",
" IF( y+y >= x*x ) EXIT\n",
" END DO\n",
" fn_val = r+x\n",
" IF( hz <= 0 ) fn_val = -fn_val\n",
" RETURN\n",
" END IF\n",
" x = hz * wn(iz)\n",
" IF( fn(iz) + uni( )*(fn(iz-1)-fn(iz)) < EXP(-half*x*x) ) THEN\n",
" fn_val = x\n",
" RETURN\n",
" END IF\n",
" hz = shr3( )\n",
" iz = IAND( hz, 127 )\n",
" IF( ABS( hz ) < kn(iz) ) THEN\n",
" fn_val = hz * wn(iz)\n",
" RETURN\n",
" END IF\n",
" END DO\n",
" END IF\n",
" RETURN\n",
"END FUNCTION rnor\n",
"\n",
"\n",
"\n",
"! Generate random exponentials\n",
"FUNCTION rexp( ) RESULT( fn_val )\n",
" REAL(DP) :: fn_val\n",
"\n",
" REAL(DP) :: x\n",
"\n",
" IF( .NOT. initialized ) CALL Zigset( jsr )\n",
" jz = shr3( )\n",
" iz = IAND( jz, 255 )\n",
" IF( ABS( jz ) < ke(iz) ) THEN\n",
" fn_val = ABS(jz) * we(iz)\n",
" RETURN\n",
" END IF\n",
" DO\n",
" IF( iz == 0 ) THEN\n",
" fn_val = 7.69711 - LOG( uni( ) )\n",
" RETURN\n",
" END IF\n",
" x = ABS( jz ) * we(iz)\n",
" IF( fe(iz) + uni( )*(fe(iz-1) - fe(iz)) < EXP( -x ) ) THEN\n",
" fn_val = x\n",
" RETURN\n",
" END IF\n",
" jz = shr3( )\n",
" iz = IAND( jz, 255 )\n",
" IF( ABS( jz ) < ke(iz) ) THEN\n",
" fn_val = ABS( jz ) * we(iz)\n",
" RETURN\n",
" END IF\n",
" END DO\n",
" RETURN\n",
"END FUNCTION rexp\n",
"\n",
"END MODULE ziggurat"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Writing ziggurat.f90\n"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The discrete OU Process code:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%file genOU.f90\n",
"\n",
"module main_mod\n",
" use ziggurat, only: rnor\n",
"\n",
"contains\n",
" function crandn() result(y)\n",
" complex(kind=8) y\n",
"\n",
" y = dcmplx(rnor() , rnor()) / dsqrt(2.0d0)\n",
"\n",
" end function crandn\n",
"\n",
" subroutine discreteprocess(uarr, N, u0, sig, F)\n",
" implicit none\n",
" integer :: N, i\n",
" complex(kind=8):: u0, sig, F, uarr(N)\n",
" complex(kind=8) :: u\n",
" intent(out) uarr\n",
"\n",
" u = u0\n",
"\n",
" do i = 1, N\n",
"\n",
" u = F * u + sig * crandn()\n",
" uarr(i) = u\n",
"\n",
" end do\n",
"\n",
" end subroutine discreteprocess\n",
"end module"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Writing genOU.f90\n"
]
}
],
"prompt_number": 3
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%file Makefile\n",
"\n",
"FC = gfortran\n",
"FFLAGS = -O3\n",
"\n",
"genOU.so : ziggurat.o genOU.f90\n",
"\tf2py -c -m genOU $^\n",
"\n",
"%.o : %.f90\n",
"\t$(FC) $(FFLAGS) -c $<\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Writing Makefile\n"
]
}
],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"!make 2> /dev/null > /dev/null\n",
"from genOU import main_mod as mm\n",
"ou_fortran = mm.discreteprocess"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Python Solver"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%pylab\n",
"\n",
"from random import random, gauss\n",
"\n",
"\n",
"def crandn():\n",
" return (1.0j * gauss(0,1) + gauss(0,1)) / 2.0**.5\n",
" \n",
"def ou_python(N, u0, sig, F):\n",
" ret = zeros(N, dtype=np.complex128)\n",
" u = u0\n",
" for i in xrange(N):\n",
" u = F * u + sig * crandn()\n",
" ret[i] = 1.0"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Using matplotlib backend: MacOSX\n",
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"prompt_number": 6
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Timing for fortran and python"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"N, u0, sig, F = 100000, 1.0, 1.0, 1.0\n",
"\n",
"\n",
"%timeit ou_fortran(N, u0, sig, F)\n",
"%timeit ou_python(N, u0, sig, F)\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"100 loops, best of 3: 5.41 ms per loop\n",
"1 loops, best of 3: 522 ms per loop"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
}
],
"prompt_number": 7
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Timing for julia"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%bash\n",
"\n",
"julia << EOF\n",
"\n",
"\n",
"\n",
"crandn() = (randn() + 1.0im * randn()) /2^.5\n",
"\n",
"function genOUJulia(N, u0, sig, F)\n",
" ret = Array(Complex128, N)\n",
" u = complex128(u0) # cast this type explicitly\n",
" for i=1:N\n",
" u = F * u + sig * crandn()\n",
" ret[i] = u\n",
" end\n",
" ret\n",
"end\n",
"\n",
"N = 100000\n",
"\n",
"a = genOUJulia(N, 1.0, 1.0, 1.0); # do jit\n",
"@time genOUJulia(N, 1.0, 1.0, 1.0);\n",
"\n",
"EOF"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"elapsed time: 0.011788437 seconds (1613896 bytes allocated)\n"
]
}
],
"prompt_number": 9
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Summary"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Python takes > 500 ms\n",
"\n",
"Julia takes 11 ms\n",
"\n",
"Fortran takes 5 ms"
]
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment