Skip to content

fixes to find_MAP #2468

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.

Sign up for GitHub

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

Closed
bwengals wants to merge 13 commits into from
Closed

fixes to find_MAP #2468

Changes from 1 commit
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
8d77aa4
fixes to find_MAP
bwengals Aug 1, 2017
4bd6f7e
Merge remote-tracking branch 'bwengals/find_MAP-float32' into mletest
kyleabeauchamp Aug 7, 2017
4647c3e
Add failing test case for 2482
kyleabeauchamp Aug 8, 2017
5b87e36
Merge remote-tracking branch 'upstream/master' into mletest
kyleabeauchamp Aug 12, 2017
f06ab8f
Fix some test failures
kyleabeauchamp Aug 12, 2017
757af19
Add logp_nojac and logp_sum
aseyboldt Aug 12, 2017
537d0bb
Merge pull request #1 from kyleabeauchamp/mletest
bwengals Aug 13, 2017
d28dce0
Resolve conflicts
kyleabeauchamp Aug 14, 2017
3f3a871
Merge master
kyleabeauchamp Aug 16, 2017
87f9ed2
Fix
kyleabeauchamp Aug 16, 2017
a8e42bd
Merge pull request #2 from kyleabeauchamp/mlefix
bwengals Aug 16, 2017
ba18757
Fix rtol
kyleabeauchamp Aug 16, 2017
3105778
Merge pull request #3 from kyleabeauchamp/mlefix
bwengals Aug 16, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
282 changes: 99 additions & 183 deletions pymc3/tuning/starting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,23 +3,27 @@

@author: johnsalvatier
'''
from scipy import optimize
import numpy as np
from numpy import isfinite, nan_to_num, logical_not
import pymc3 as pm
import time
from ..vartypes import discrete_types, typefilter
from ..model import modelcontext, Point
from ..theanof import inputvars
from ..blocking import DictToArrayBijection, ArrayOrdering
from ..util import update_start_vals
from tqdm import tqdm

from pymc3.vartypes import discrete_types, typefilter
from pymc3.model import modelcontext, Point
from pymc3.blocking import ArrayOrdering, DictToArrayBijection
from pymc3.theanof import inputvars, floatX
from pymc3.util import update_start_vals
from scipy.optimize import minimize


from inspect import getargspec

__all__ = ['find_MAP']

def find_MAP(start=None, vars=None, fmin=None,
return_raw=False, model=None, live_disp=False, callback=None, *args, **kwargs):

def find_MAP(start=None, vars=None, method=None, progressbar=True, return_raw=False,
model=None, maxeval=50000, callback=None, *args, **kwargs):
"""
Sets state to the local maximum a posteriori point given a model.
Current default of fmin_Hessian does not deal well with optimizing close
Expand All @@ -30,21 +34,21 @@ def find_MAP(start=None, vars=None, fmin=None,
start : `dict` of parameter values (Defaults to `model.test_point`)
vars : list
List of variables to set to MAP point (Defaults to all continuous).
fmin : function
Optimization algorithm (Defaults to `scipy.optimize.fmin_bfgs` unless
method : string or callable

@kyleabeauchamp kyleabeauchamp Aug 4, 2017

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Is it really that much more convenient to allow strings here? In particular, is this minor convenience more important than preserving the API?

@bwengals bwengals Aug 4, 2017

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I'm not sure. Some benefits I can think of is it matches up a bit better with sample and fit, which specify the method with strings. Also, minimize can take a homebrew optimizer, which may be nice. It also gives a consistent interface for all scipys optimizers. The different fmin_*` functions are sometimes a little bit different. Draw back is breaking peoples code...

@kyleabeauchamp kyleabeauchamp Aug 4, 2017

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Either way SGTM, just wanted to bring up the API stability concern.

@bwengals bwengals Aug 5, 2017

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

yeah I'm not sure how folks feel about it

Optimization algorithm (Defaults to `BFGS` unless
discrete variables are specified in `vars`, then
`scipy.optimize.fmin_powell` which will perform better).
return_raw : Bool
Whether to return extra value returned by fmin (Defaults to `False`)
`Powell` which will perform better).
progressbar : bool
Whether or not to display a progress bar in the command line.
return_raw : bool
Whether to return extra values returned by fmin (Defaults to `False`)
model : Model (optional if in `with` context)
live_disp : Bool
Display table tracking optimization progress when run from within
an IPython notebook.
maxeval : int
The maximum number of times the posterior distribution is evaluated.
callback : callable
Callback function to pass to scipy optimization routine. Overrides
live_disp if callback is given.
Callback function to pass to scipy optimization routine.
*args, **kwargs
Extra args passed to fmin
Extra args passed to fmin.
"""
model = modelcontext(model)
if start is None:
Expand All @@ -69,62 +73,46 @@ def find_MAP(start=None, vars=None, fmin=None,
except AttributeError:
gradient_avail = False

if disc_vars or not gradient_avail :
if disc_vars or not gradient_avail:
pm._log.warning("Warning: gradient not available." +
"(E.g. vars contains discrete variables). MAP " +
"estimates may not be accurate for the default " +
"parameters. Defaulting to non-gradient minimization " +
"fmin_powell.")
fmin = optimize.fmin_powell
"'Powell'.")
method = "Powell"

if fmin is None:
if method is None:
if disc_vars:
fmin = optimize.fmin_powell
method = "Powell"
else:
fmin = optimize.fmin_bfgs
method = "BFGS"

allinmodel(vars, model)

start = Point(start, model=model)
bij = DictToArrayBijection(ArrayOrdering(vars), start)
logp_func = bij.mapf(model.fastlogp)
x0 = bij.map(start)

logp = bij.mapf(model.fastlogp)
def logp_o(point):
return nan_to_high(-logp(point))

# Check to see if minimization function actually uses the gradient
if 'fprime' in getargspec(fmin).args:
dlogp = bij.mapf(model.fastdlogp(vars))
def grad_logp_o(point):
return nan_to_num(-dlogp(point))

if live_disp and callback is None:
callback = Monitor(bij, logp_o, model, grad_logp_o)

r = fmin(logp_o, bij.map(start), fprime=grad_logp_o, callback=callback, *args, **kwargs)
if method in ["CG", "BFGS", "Newton-CG", "L-BFGS-B", "TNC",
"SLSQP", "dogleg", "trust-ncg"]:

@junpenglao junpenglao Aug 2, 2017

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This would mean if I provide a custom method I can not use the gradient?

@bwengals bwengals Aug 2, 2017

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

oh good catch, it would be nice to be able to provide a custom method. I'll try and fix that.

dlogp_func = bij.mapf(model.fastdlogp(vars))
cost_func = CostFuncWrapper(maxeval, progressbar, logp_func, dlogp_func)
compute_gradient = True
else:
if live_disp and callback is None:
callback = Monitor(bij, logp_o, dlogp=None)

# Check to see if minimization function uses a starting value
if 'x0' in getargspec(fmin).args:
r = fmin(logp_o, bij.map(start), callback=callback, *args, **kwargs)
else:
r = fmin(logp_o, callback=callback, *args, **kwargs)
cost_func = CostFuncWrapper(maxeval, progressbar, logp_func)
compute_gradient = False

if isinstance(r, tuple):
mx0 = r[0]
else:
mx0 = r

if live_disp:
try:
callback.update(mx0)
except:
pass

try:
r = minimize(cost_func, x0, method=method, jac=compute_gradient, *args, **kwargs)
mx0 = r["x"]
except (KeyboardInterrupt, StopIteration) as e:
mx0, r = cost_func.previous_x, None
cost_func.progress.close()
if isinstance(e, StopIteration):
pm._log.info(e)
finally:
cost_func.progress.close()
mx = bij.rmap(mx0)

allfinite_mx0 = allfinite(mx0)
Expand Down Expand Up @@ -169,148 +157,76 @@ def message(name, values):
"density. 2) your distribution logp's are " +
"properly specified. Specific issues: \n" +
specific_errors)
mx = {v.name: mx[v.name].astype(v.dtype) for v in model.vars}

vars = model.unobserved_RVs
mx = {var.name: value for var, value in zip(vars, model.fastfn(vars)(mx))}
if return_raw:
return mx, r
else:
return mx


def allfinite(x):
return np.all(isfinite(x))


def nan_to_high(x):
return np.where(isfinite(x), x, 1.0e100)


def allinmodel(vars, model):
notin = [v for v in vars if v not in model.vars]
if notin:
raise ValueError("Some variables not in the model: " + str(notin))



class Monitor(object):
def __init__(self, bij, logp, model, dlogp=None):
try:
from IPython.display import display
from ipywidgets import HTML, VBox, HBox, FlexBox
self.prog_table = HTML(width='100%')
self.param_table = HTML(width='100%')
r_col = VBox(children=[self.param_table], padding=3, width='100%')
l_col = HBox(children=[self.prog_table], padding=3, width='25%')
self.hor_align = FlexBox(children = [l_col, r_col], width='100%', orientation='vertical')
display(self.hor_align)
self.using_notebook = True
self.update_interval = 1
except:
self.using_notebook = False
self.update_interval = 2

self.iters = 0
self.bij = bij
self.model = model
self.fn = model.fastfn(model.unobserved_RVs)
self.logp = logp
self.dlogp = dlogp
self.t_initial = time.time()
self.t0 = self.t_initial
self.paramtable = {}
class CostFuncWrapper(object):
def __init__(self, maxeval=5000, progressbar=True, logp_func=None, dlogp_func=None):
self.t0 = time.time()
self.n_eval = 0
self.maxeval = maxeval
self.logp_func = logp_func
if dlogp_func is None:
self.use_gradient = False
self.desc = 'lp = {:,.5g}'
else:
self.dlogp_func = dlogp_func
self.use_gradient = True
self.desc = 'lp = {:,.5g}, ||grad|| = {:,.5g}'
self.previous_x = None
self.progress = tqdm(total=maxeval, disable=not progressbar)

def __call__(self, x):
self.iters += 1
if time.time() - self.t0 > self.update_interval or self.iters == 1:
self.update(x)

def update(self, x):
self._update_progtable(x)
self._update_paramtable(x)
if self.using_notebook:
self._display_notebook()
self.t0 = time.time()

def _update_progtable(self, x):
s = time.time() - self.t_initial
hours, remainder = divmod(int(s), 3600)
minutes, seconds = divmod(remainder, 60)
self.t_elapsed = "{:2d}h{:2d}m{:2d}s".format(hours, minutes, seconds)
self.logpost = -1.0*np.float(self.logp(x))
self.dlogpost = np.linalg.norm(self.dlogp(x))

def _update_paramtable(self, x):
var_state = self.fn(self.bij.rmap(x))
for var, val in zip(self.model.unobserved_RVs, var_state):
if not var.name.endswith("_"):
valstr = format_values(val)
self.paramtable[var.name] = {"size": val.size, "valstr": valstr}

def _display_notebook(self):
## Progress table
html = r"""<style type="text/css">
table { border-collapse:collapse }
.tg {border-collapse:collapse;border-spacing:0;border:none;}
.tg td{font-family:Arial, sans-serif;font-size:14px;padding:3px 3px;border-style:solid;border-width:0px;overflow:hidden;word-break:normal;}
.tg th{Impact, Charcoal, sans-serif;font-size:13px;font-weight:bold;padding:3px 3px;border-style:solid;border-width:0px;overflow:hidden;word-break:normal; background-color:#0E688A;color:#ffffff;}
.tg .tg-vkoh{white-space:pre;font-weight:normal;font-family:"Lucida Console", Monaco, monospace !important; background-color:#ffffff;color:#000000}
.tg .tg-suao{font-weight:bold;font-family:"Lucida Console", Monaco, monospace !important;background-color:#0E688A;color:#ffffff;}
"""
html += r"""
</style>
<table class="tg" style="undefined;">
<col width="400px" />
<tr>
<th class= "tg-vkoh">Time Elapsed: {:s}</th>
</tr>
<tr>
<th class= "tg-vkoh">Iteration: {:d}</th>
</tr>
<tr>
<th class= "tg-vkoh">Log Posterior: {:.3f}</th>
</tr>
""".format(self.t_elapsed, self.iters, self.logpost)
if self.dlogp is not None:
html += r"""
<tr>
<th class= "tg-vkoh">||grad||: {:.3f}</th>
</tr>""".format(self.dlogpost)
html += "</table>"
self.prog_table.value = html
## Parameter table
html = r"""<style type="text/css">
.tg .tg-bgft{font-weight:normal;font-family:"Lucida Console", Monaco, monospace !important;background-color:#0E688A;color:#ffffff;}
.tg td{font-family:Arial, sans-serif;font-size:12px;padding:3px 3px;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#504A4E;color:#333;background-color:#fff;word-wrap: break-word;}
.tg th{Impact, Charcoal, sans-serif;font-size:13px;font-weight:bold;padding:3px 3px;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;border-color:#504A4E;background-color:#0E688A;color:#ffffff;}
</style>
<table class="tg" style="undefined;">
<col width="130px" />
<col width="50px" />
<col width="600px" />
<tr>
<th class="tg">Parameter</th>
<th class="tg">Size</th>
<th class="tg">Current Value</th>
</tr>
"""
for var, values in self.paramtable.items():
html += r"""
<tr>
<td class="tg-bgft">{:s}</td>
<td class="tg-vkoh">{:d}</td>
<td class="tg-vkoh">{:s}</td>
</tr>
""".format(var, values["size"], values["valstr"])
html += "</table>"
self.param_table.value = html


def format_values(val):
fmt = "{:8.3f}"
if val.size == 1:
return fmt.format(np.float(val))
elif val.size < 9:
return "[" + ", ".join([fmt.format(v) for v in val]) + "]"
else:
start = "[" + ", ".join([fmt.format(v) for v in val[:4]])
end = ", ".join([fmt.format(v) for v in val[-4:]]) +"]"
return start + ", ... , " + end
neg_value = np.float64(self.logp_func(pm.floatX(x)))
value = -1.0 * nan_to_high(neg_value)
if self.use_gradient:
neg_grad = self.dlogp_func(pm.floatX(x))
if np.all(np.isfinite(neg_grad)):
self.previous_x = x
grad = nan_to_num(-1.0*neg_grad)
grad = grad.astype(np.float64)
else:
self.previous_x = x
grad = None

if self.n_eval % 10 == 0:
self.update_progress_desc(neg_value, grad)

if self.n_eval > self.maxeval:
self.update_progress_desc(neg_value, grad)
self.progress.close()
raise StopIteration

self.n_eval += 1
self.progress.update(1)

if self.use_gradient:
return value, grad
else:
return value

def update_progress_desc(self, neg_value, grad=None):
if grad is None:
self.progress.set_description(self.desc.format(neg_value))
else:
norm_grad = np.linalg.norm(grad)
self.progress.set_description(self.desc.format(neg_value, norm_grad))