diff --git a/docs/qe_apidoc.py b/docs/qe_apidoc.py index ffc345b5e..3c4c87c9d 100644 --- a/docs/qe_apidoc.py +++ b/docs/qe_apidoc.py @@ -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} @@ -133,6 +142,7 @@ game_theory markov + optimize random tools util @@ -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)) @@ -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)) @@ -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") @@ -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]) diff --git a/docs/source/index.rst b/docs/source/index.rst index 9e969cea5..6dcf25317 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -13,6 +13,7 @@ mainly used by developers internal to the package. game_theory markov + optimize random tools util diff --git a/docs/source/optimize.rst b/docs/source/optimize.rst new file mode 100644 index 000000000..7aa11d50d --- /dev/null +++ b/docs/source/optimize.rst @@ -0,0 +1,7 @@ +Optimize +======== + +.. toctree:: + :maxdepth: 2 + + optimize/scalar_maximization diff --git a/docs/source/optimize/scalar_maximization.rst b/docs/source/optimize/scalar_maximization.rst new file mode 100644 index 000000000..af15d0538 --- /dev/null +++ b/docs/source/optimize/scalar_maximization.rst @@ -0,0 +1,7 @@ +scalar_maximization +=================== + +.. automodule:: quantecon.optimize.scalar_maximization + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/tools.rst b/docs/source/tools.rst index ef3923f93..7b8456710 100644 --- a/docs/source/tools.rst +++ b/docs/source/tools.rst @@ -11,6 +11,7 @@ Tools tools/distributions tools/ecdf tools/estspec + tools/filter tools/graph_tools tools/gridtools tools/ivp diff --git a/docs/source/tools/filter.rst b/docs/source/tools/filter.rst new file mode 100644 index 000000000..f102fb19b --- /dev/null +++ b/docs/source/tools/filter.rst @@ -0,0 +1,7 @@ +filter +====== + +.. automodule:: quantecon.filter + :members: + :undoc-members: + :show-inheritance: diff --git a/quantecon/__init__.py b/quantecon/__init__.py index 7600dc6e4..1428a7f35 100644 --- a/quantecon/__init__.py +++ b/quantecon/__init__.py @@ -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 diff --git a/quantecon/optimize/__init__.py b/quantecon/optimize/__init__.py new file mode 100644 index 000000000..4937fb738 --- /dev/null +++ b/quantecon/optimize/__init__.py @@ -0,0 +1,6 @@ +""" +Initialization of the optimize subpackage +""" + +from .scalar_maximization import brent_max + diff --git a/quantecon/optimize/scalar_maximization.py b/quantecon/optimize/scalar_maximization.py new file mode 100644 index 000000000..9eff745fb --- /dev/null +++ b/quantecon/optimize/scalar_maximization.py @@ -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 diff --git a/quantecon/optimize/tests/test_scalar_max.py b/quantecon/optimize/tests/test_scalar_max.py new file mode 100644 index 000000000..b05c00113 --- /dev/null +++ b/quantecon/optimize/tests/test_scalar_max.py @@ -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__) + + diff --git a/quantecon/version.py b/quantecon/version.py index e38e023ca..aad6dede3 100644 --- a/quantecon/version.py +++ b/quantecon/version.py @@ -1,4 +1,4 @@ """ This is a VERSION file and should NOT be manually altered """ -version = '0.3.7' \ No newline at end of file +version = '0.3.8' \ No newline at end of file diff --git a/setup.py b/setup.py index 2df3ab4c8..efbe6f05c 100644 --- a/setup.py +++ b/setup.py @@ -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', @@ -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',