AdityaSavara / PEUQSE

Parameter estimation for complex physical problems often suffers from finding ‘solutions’ that are not physically realistic. The PEUQSE software provides tools for finding physically realistic parameter estimates, graphs of the parameter estimate positions within parameter space, and plots of the final simulation results.
13 stars 5 forks source link

mumpce plots failing due to numerical errors causing linalg error #73

Open AdityaSavara opened 4 years ago

AdityaSavara commented 4 years ago

To run the example attached to this post in the following zip file, use runfile_Example8_optimization_mcmc.py. It should take under a minute to run. 200323_2.2a-MiniExample.zip

I received this error:

Traceback (most recent call last):

  File "<ipython-input-25-550d60b1b48a>", line 1, in <module>
    runfile('C:/Users/fvs/Desktop/Temp/Junk/b/BayesKin/200323_2.2a/runfile_Example8_optimization_mcmc.py', wdir='C:/Users/fvs/Desktop/Temp/Junk/b/BayesKin/200323_2.2a')

  File "C:\Users\fvs\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 678, in runfile
    execfile(filename, namespace)

  File "C:\Users\fvs\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 106, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/fvs/Desktop/Temp/Junk/b/BayesKin/200323_2.2a/runfile_Example8_optimization_mcmc.py", line 71, in <module>
    PE_object.createAllPlots() #This function calls each of the below functions.

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a\CheKiPEUQ\InverseProblem.py", line 683, in createAllPlots
    self.createMumpcePlots()

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a\CheKiPEUQ\InverseProblem.py", line 670, in createMumpcePlots
    figureObject_beta.mumpce_plots(model_parameter_info = self.UserInput.model_parameter_info, active_parameters = active_parameters, pairs_of_parameter_indices = pairs_of_parameter_indices, posterior_mu_vector = posterior_mu_vector, posterior_cov_matrix = posterior_cov_matrix, prior_mu_vector = np.array(self.UserInput.mu_prior), prior_cov_matrix = self.UserInput.covmat_prior, contour_settings_custom = self.UserInput.contour_settings_custom)

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a\CheKiPEUQ\plotting_functions.py", line 49, in mumpce_plots
    mumpceSolutionsObject = mumpceSolution.Solution(posterior_mu_vector, posterior_cov_matrix, initial_x=prior_mu_vector, initial_covariance=prior_cov_matrix)

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a\CheKiPEUQ\mumpce\solution.py", line 18, in __init__
    self.alpha = np.linalg.cholesky(covariance_x) #: The lower triangular decomposition :math:`\alpha` where :math:`\Sigma = \alpha \alpha^{\text{T}}`

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\numpy\linalg\linalg.py", line 759, in cholesky
    r = gufunc(a, signature=signature, extobj=extobj)

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\numpy\linalg\linalg.py", line 100, in _raise_linalgerror_nonposdef
    raise LinAlgError("Matrix is not positive definite")

LinAlgError: Matrix is not positive definite

I think it is because a matrix (probably posterior cov) had negative or zero values due to numerical errors or insufficient sampling. https://stackoverflow.com/questions/21604498/numpy-cholesky-decomposition-linalgerror

I think that what we should do is give people a warning while also forcing the matrix to become positive definite.

For forcing: There is some r code, don't now if it's in python https://rdrr.io/cran/stsm/man/force-defpos.html here is a small python script: https://stackoverflow.com/questions/43238173/python-convert-matrix-to-positive-semi-definite

For warning: I found that for my case increasing hte mcmc step length and increasing the mcmc_modulate_accept_probability to 1000 solved my problem. More details of what I think should be put in the warning will be in the comments below.

AdityaSavara commented 4 years ago

So what we need to do is... a) Check for positive definite in posterior_cov_matrix, and if it fails then give user warning to use mcmc_modulate_accept_probability (they actually might need to go much higher than what the instructions provided say, I had to go as high as 1000 with eight parameters, so we need to revise our instructions). Also, to consider to increase/decrease the mcmc step size and increase the sample length, change the starting point. b) later, put code to force positive definite. I think that this may be happening due to covariances of zero. One thing we can do is to put in an arbitrarily small covariance that is something like 1/100 * mu_var that way it is effectively zero on the scale of the mumpce plot. c) even later, give user suggestions of how to change mcmc step size. In future, I want to make mcmc step size an array instead of a single value (separate step sizes for each variable). I think pymc does that automatically. So another option for c is to wait until the pymc engine is inluded as an option.

AdityaSavara commented 4 years ago

Need to put a try and except mumpce plots colorbars. for a totally different reason, I got this (I believe this was related to not having enough contour lines in the range displayed, for the settings I had).


Traceback (most recent call last):

  File "<ipython-input-22-04d16fe900a9>", line 1, in <module>
    runfile('C:/Users/fvs/Desktop/Temp/Junk/b/BayesKin/200323_2.2a - Copy3reducedParamOptResult/runfile_Example8_optimization_mcmc_posterior_plot.py', wdir='C:/Users/fvs/Desktop/Temp/Junk/b/BayesKin/200323_2.2a - Copy3reducedParamOptResult')

  File "C:\Users\fvs\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 678, in runfile
    execfile(filename, namespace)

  File "C:\Users\fvs\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 106, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/fvs/Desktop/Temp/Junk/b/BayesKin/200323_2.2a - Copy3reducedParamOptResult/runfile_Example8_optimization_mcmc_posterior_plot.py", line 76, in <module>
    PE_object.createAllPlots() #This function calls each of the below functions.

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a - Copy3reducedParamOptResult\CheKiPEUQ\InverseProblem.py", line 687, in createAllPlots
    self.createMumpcePlots()

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a - Copy3reducedParamOptResult\CheKiPEUQ\InverseProblem.py", line 674, in createMumpcePlots
    figureObject_beta.mumpce_plots(model_parameter_info = self.UserInput.model_parameter_info, active_parameters = active_parameters, pairs_of_parameter_indices = pairs_of_parameter_indices, posterior_mu_vector = posterior_mu_vector, posterior_cov_matrix = posterior_cov_matrix, prior_mu_vector = np.array(self.UserInput.mu_prior), prior_cov_matrix = self.UserInput.covmat_prior, contour_settings_custom = self.UserInput.contour_settings_custom)

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a - Copy3reducedParamOptResult\CheKiPEUQ\plotting_functions.py", line 87, in mumpce_plots
    mumpceProjectObject.plot_pdfs(mumpceProjectObject.pairsOfParameterIndices, contour_settings_custom = contour_settings_custom)

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a - Copy3reducedParamOptResult\CheKiPEUQ\mumpce\Project.py", line 1017, in plot_pdfs
    self._single_pdf_plot(factors=factors,ax=ax, fig=fig, contour_settings_custom = contour_settings_custom)

  File "C:\Users\fvs\Desktop\Temp\Junk\b\BayesKin\200323_2.2a - Copy3reducedParamOptResult\CheKiPEUQ\mumpce\Project.py", line 975, in _single_pdf_plot
    poste_colorbar = fig.colorbar(cposte,ax=ax,orientation='horizontal')

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\matplotlib\figure.py", line 2215, in colorbar
    cb = cbar.colorbar_factory(cax, mappable, **cb_kw)

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\matplotlib\colorbar.py", line 1640, in colorbar_factory
    cb = Colorbar(cax, mappable, **kwargs)

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\matplotlib\colorbar.py", line 1173, in __init__
    ColorbarBase.__init__(self, ax, **kw)

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\matplotlib\colorbar.py", line 460, in __init__
    self.draw_all()

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\matplotlib\colorbar.py", line 493, in draw_all
    self._config_axes(X, Y)

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\matplotlib\colorbar.py", line 702, in _config_axes
    xy = self._outline(X, Y)

  File "C:\Users\fvs\AppData\Roaming\Python\Python36\site-packages\matplotlib\colorbar.py", line 749, in _outline
    x = X.T.reshape(-1)[ii]

IndexError: index 2 is out of bounds for axis 0 with size 2

Suffice to say there is a work around by putting a try and except around this block:

if contour_settings_custom["colorbars"] == True: