Skip to content

adding jitted scalar maximization routine, first build #416

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 12 commits into from
Jun 29, 2018
Merged
31 changes: 29 additions & 2 deletions docs/qe_apidoc.py
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,15 @@
:show-inheritance:
"""

optimize_module_template = """{mod_name}
{equals}

.. automodule:: quantecon.optimize.{mod_name}
:members:
:undoc-members:
:show-inheritance:
"""

random_module_template = """{mod_name}
{equals}

Expand Down Expand Up @@ -133,6 +142,7 @@

game_theory
markov
optimize
random
tools
util
Expand Down Expand Up @@ -211,6 +221,12 @@ def model_tool():
# Alphabetize
markov.sort()

# list file names with optimize
optimize_files = glob("../quantecon/optimize/[a-z0-9]*.py")
optimize = list(map(lambda x: x.split('/')[-1][:-3], optimize_files))
# Alphabetize
optimize.sort()

# list file names with random
random_files = glob("../quantecon/random/[a-z0-9]*.py")
random = list(map(lambda x: x.split('/')[-1][:-3], random_files))
Expand All @@ -230,7 +246,7 @@ def model_tool():
# Alphabetize
util.sort()

for folder in ["game_theory", "markov", "random", "tools", "util"]:
for folder in ["game_theory", "markov", "optimize", "random", "tools", "util"]:
if not os.path.exists(source_join(folder)):
os.makedirs(source_join(folder))

Expand All @@ -257,6 +273,13 @@ def model_tool():
equals = "=" * len(mod)
f.write(markov_module_template.format(mod_name=mod, equals=equals))

# Write file for each optimize file
for mod in optimize:
new_path = os.path.join("source", "optimize", mod + ".rst")
with open(new_path, "w") as f:
equals = "=" * len(mod)
f.write(optimize_module_template.format(mod_name=mod, equals=equals))

# Write file for each random file
for mod in random:
new_path = os.path.join("source", "random", mod + ".rst")
Expand Down Expand Up @@ -284,24 +307,28 @@ def model_tool():

gt = "game_theory/" + "\n game_theory/".join(game_theory)
mark = "markov/" + "\n markov/".join(markov)
opti = "optimize/" + "\n optimize/".join(optimize)
rand = "random/" + "\n random/".join(random)
tlz = "tools/" + "\n tools/".join(tools)
utls = "util/" + "\n util/".join(util)
#-TocTree-#
toc_tree_list = {"game_theory": gt,
"markov": mark,
"optimize" : opti,
"tools": tlz,
"random": rand,
"util": utls,
}

for f_name in ("game_theory", "markov", "random", "tools", "util"):
for f_name in ("game_theory", "markov", "optimize", "random", "tools", "util"):
with open(source_join(f_name + ".rst"), "w") as f:
m_name = f_name
if f_name == "game_theory":
f_name = "Game Theory" #Produce Nicer Title for Game Theory Module
if f_name == "util":
f_name = "Utilities" #Produce Nicer Title for Utilities Module
if f_name == "optimize":
f_name = "Optimize"
temp = split_file_template.format(name=f_name.capitalize(),
equals="="*len(f_name),
files=toc_tree_list[m_name])
Expand Down
1 change: 1 addition & 0 deletions docs/source/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ mainly used by developers internal to the package.

game_theory
markov
optimize
random
tools
util
Expand Down
7 changes: 7 additions & 0 deletions docs/source/optimize.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
Optimize
========

.. toctree::
:maxdepth: 2

optimize/scalar_maximization
7 changes: 7 additions & 0 deletions docs/source/optimize/scalar_maximization.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
scalar_maximization
===================

.. automodule:: quantecon.optimize.scalar_maximization
:members:
:undoc-members:
:show-inheritance:
1 change: 1 addition & 0 deletions docs/source/tools.rst
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ Tools
tools/distributions
tools/ecdf
tools/estspec
tools/filter
tools/graph_tools
tools/gridtools
tools/ivp
Expand Down
7 changes: 7 additions & 0 deletions docs/source/tools/filter.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
filter
======

.. automodule:: quantecon.filter
:members:
:undoc-members:
:show-inheritance:
1 change: 1 addition & 0 deletions quantecon/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
from . import game_theory
from . import quad
from . import random
from . import optimize

#-Objects-#
from .compute_fp import compute_fixed_point
Expand Down
6 changes: 6 additions & 0 deletions quantecon/optimize/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
"""
Initialization of the optimize subpackage
"""

from .scalar_maximization import brent_max

144 changes: 144 additions & 0 deletions quantecon/optimize/scalar_maximization.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,144 @@
import numpy as np
from numba import jit, njit

@njit
def brent_max(func, a, b, args=(), xtol=1e-5, maxiter=500):
"""
Uses a jitted version of the maximization routine from SciPy's fminbound.
The algorithm is identical except that it's been switched to maximization
rather than minimization, and the tests for convergence have been stripped
out to allow for jit compilation.

Note that the input function `func` must be jitted or the call will fail.

Parameters
----------
func : jitted function
a : scalar
Lower bound for search
b : scalar
Upper bound for search
args : tuple, optional
Extra arguments passed to the objective function.
maxiter : int, optional
Maximum number of iterations to perform.
xtol : float, optional
Absolute error in solution `xopt` acceptable for convergence.

Returns
-------
xf : float
The maximizer
fval : float
The maximum value attained
info : tuple
A tuple of the form (status_flag, num_iter). Here status_flag
indicates whether or not the maximum number of function calls was
attained. A value of 0 implies that the maximum was not hit.
The value `num_iter` is the number of function calls.

Example
-------

```
@njit
def f(x):
return -(x + 2.0)**2 + 1.0

xf, fval, info = maximize_scalar(f, -2, 2)
```

"""

maxfun = maxiter
status_flag = 0

sqrt_eps = np.sqrt(2.2e-16)
golden_mean = 0.5 * (3.0 - np.sqrt(5.0))

fulc = a + golden_mean * (b - a)
nfc, xf = fulc, fulc
rat = e = 0.0
x = xf
fx = -func(x, *args)
num = 1

ffulc = fnfc = fx
xm = 0.5 * (a + b)
tol1 = sqrt_eps * np.abs(xf) + xtol / 3.0
tol2 = 2.0 * tol1

while (np.abs(xf - xm) > (tol2 - 0.5 * (b - a))):
golden = 1
# Check for parabolic fit
if np.abs(e) > tol1:
golden = 0
r = (xf - nfc) * (fx - ffulc)
q = (xf - fulc) * (fx - fnfc)
p = (xf - fulc) * q - (xf - nfc) * r
q = 2.0 * (q - r)
if q > 0.0:
p = -p
q = np.abs(q)
r = e
e = rat

# Check for acceptability of parabola
if ((np.abs(p) < np.abs(0.5*q*r)) and (p > q*(a - xf)) and
(p < q * (b - xf))):
rat = (p + 0.0) / q
x = xf + rat

if ((x - a) < tol2) or ((b - x) < tol2):
si = np.sign(xm - xf) + ((xm - xf) == 0)
rat = tol1 * si
else: # do a golden section step
golden = 1

if golden: # Do a golden-section step
if xf >= xm:
e = a - xf
else:
e = b - xf
rat = golden_mean*e

if rat == 0:
si = np.sign(rat) + 1
else:
si = np.sign(rat)

x = xf + si * np.maximum(np.abs(rat), tol1)
fu = -func(x, *args)
num += 1

if fu <= fx:
if x >= xf:
a = xf
else:
b = xf
fulc, ffulc = nfc, fnfc
nfc, fnfc = xf, fx
xf, fx = x, fu
else:
if x < xf:
a = x
else:
b = x
if (fu <= fnfc) or (nfc == xf):
fulc, ffulc = nfc, fnfc
nfc, fnfc = x, fu
elif (fu <= ffulc) or (fulc == xf) or (fulc == nfc):
fulc, ffulc = x, fu

xm = 0.5 * (a + b)
tol1 = sqrt_eps * np.abs(xf) + xtol / 3.0
tol2 = 2.0 * tol1

if num >= maxfun:
status_flag = 1
break

fval = -fx
info = status_flag, num

return xf, fval, info
58 changes: 58 additions & 0 deletions quantecon/optimize/tests/test_scalar_max.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
"""
Tests for scalar maximization.

"""
import numpy as np
from numpy.testing import assert_almost_equal
from numba import njit

from quantecon.optimize import brent_max

@njit
def f(x):
"""
A function for testing on.
"""
return -(x + 2.0)**2 + 1.0

def test_brent_max():
"""
Uses the function f defined above to test the scalar maximization
routine.
"""
true_fval = 1.0
true_xf = -2.0
xf, fval, info = brent_max(f, -2, 2)
assert_almost_equal(true_fval, fval, decimal=4)
assert_almost_equal(true_xf, xf, decimal=4)

@njit
def g(x, y):
"""
A multivariate function for testing on.
"""
return -x**2 + y

def test_brent_max():
"""
Uses the function f defined above to test the scalar maximization
routine.
"""
y = 5
true_fval = 5.0
true_xf = -0.0
xf, fval, info = brent_max(g, -10, 10, args=(y,))
assert_almost_equal(true_fval, fval, decimal=4)
assert_almost_equal(true_xf, xf, decimal=4)


if __name__ == '__main__':
import sys
import nose

argv = sys.argv[:]
argv.append('--verbose')
argv.append('--nocapture')
nose.main(argv=argv, defaultTest=__file__)


2 changes: 1 addition & 1 deletion quantecon/version.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
"""
This is a VERSION file and should NOT be manually altered
"""
version = '0.3.7'
version = '0.3.8'
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this will be set in setup.py
I will update this tomorrow morning.

3 changes: 2 additions & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ def write_version_py(filename=None):
'quantecon.game_theory',
'quantecon.game_theory.game_generators',
'quantecon.markov',
'quantecon.optimize',
'quantecon.random',
'quantecon.tests',
'quantecon.util',
Expand All @@ -113,7 +114,7 @@ def write_version_py(filename=None):
download_url='https://github.com/QuantEcon/QuantEcon.py/tarball/' + VERSION,
keywords=['quantitative', 'economics'],
install_requires=[
'numba>=0.36.2',
'numba>=0.38',
'numpy',
'scipy>=1.0.0',
'sympy',
Expand Down