rus_cubic_cmsx4_8modes.py

#%%

import numpy
import time
import scipy
import sympy
import os
os.chdir('/home/bbales2/modal')

import rus
reload(rus)
#%%
#656 samples
#%%

# basis polynomials are x^n * y^m * z^l where n + m + l <= N
N = 12

## Dimensions for TF-2
X = 0.011959#0.007753
Y = 0.013953#0.009057
Z = 0.019976#0.013199

#Sample density
density = 8700.0#4401.695921 #Ti-64-TF2

c110 = 2.0
anisotropic0 = 2.0
c440 = 1.0
c120 = -(c440 * 2.0 / anisotropic0 - c110)

# Standard deviation around each mode prediction
std0 = 5.0

# Measured resonance modes
data = numpy.array([71.25925,
75.75875,
86.478,
89.947375,
111.150125,
112.164125,
120.172125,
127.810375])

#%%

# These are the two HMC parameters
#   L is the number of timesteps to take -- use this if samples in the traceplots don't look random
#   epsilon is the timestep -- make this small enough so that pretty much all the samples are being accepted, but you
#       want it large enough that you can keep L ~ 50 -> 100 and still get independent samples
L = 50
# start epsilon at .0001 and try larger values like .0005 after running for a while
# epsilon is timestep, we want to make as large as possibe, wihtout getting too many rejects
epsilon = 0.0001

# Set this to true to debug the L and eps values
debug = False

#%%
# Initialize model and warm it up

reload(rus)

c11, anisotropic, c44 = sympy.symbols('c11 anisotropic c44')

c12 = sympy.sympify("-(c44 * 2.0 / anisotropic - c11)") # The c11 and c44 and anisotropic are the same as above

C = sympy.Matrix([[c11, c12, c12, 0, 0, 0],
                  [c12, c11, c12, 0, 0, 0],
                  [c12, c12, c11, 0, 0, 0],
                  [0, 0, 0, c44, 0, 0],
                  [0, 0, 0, 0, c44, 0],
                  [0, 0, 0, 0, 0, c44]])

hmc = rus.HMC(tol = 1e-3,
              density = density, X = X, Y = Y, Z = Z,
              resonance_modes = data, # List of resonance modes
              stiffness_matrix = C, # Stiffness matrix
              parameters = { c11 : c110, anisotropic : anisotropic0, c44 : c440, 'std' : std0 }, # Parameters
              rotations = True,
              T = 1.0)

hmc.set_labels({ c11 : 'c11', anisotropic : 'a', c44 : 'c44', 'std' : 'std' })
hmc.set_timestepping(epsilon = epsilon, L = 50)
#hmc.print_current()
#hmc.computeResolutions(1e-3)
hmc.sample(steps = 5, debug = True)
#%%
hmc.derivative_check()
#%%
# Run for a long time
hmc.set_timestepping(epsilon = epsilon * 10.0, L = 50)
hmc.sample(debug = False)#True)#False)#True)
#%%
# Plot traceplots
import matplotlib.pyplot as plt
import seaborn

for name, data1 in zip(*hmc.format_samples()):
    plt.plot(data1)
    plt.title('{0}'.format(name, numpy.mean(data1), numpy.std(data1)), fontsize = 24)
    plt.tick_params(axis='y', which='major', labelsize=16)
    plt.tick_params(axis='x', which='major', labelsize=16)
    fig = plt.gcf()
    fig.set_size_inches((10, 6))
    plt.show()
#%%
# Plot approximate posteriors as histograms
for name, data1 in zip(*hmc.format_samples()):
    data1 = data1[-14000:]
    seaborn.distplot(data1, kde = False, fit = scipy.stats.norm)
    plt.title('{0}, $\mu$ = {1:0.3f}, $\sigma$ = {2:0.3f}'.format(name, numpy.mean(data1), numpy.std(data1)), fontsize = 36)
    plt.tick_params(axis='x', which='major', labelsize=16)
    fig = plt.gcf()
    fig.set_size_inches((10, 6))
    plt.show()

#%%
# Plot posterior predictive
hmc.posterior_predictive(plot = True, lastN = 200)
plt.title('Posterior predictive', fontsize = 72)
plt.xlabel('Mode', fontsize = 48)
plt.ylabel('Computed - Measured (khz)', fontsize = 48)
plt.tick_params(axis='y', which='major', labelsize=48)
plt.tick_params(axis='x', which='major', labelsize=48)
fig = plt.gcf()
fig.set_size_inches((24, 16))
plt.savefig('dec2/cmsxlow/posteriorpredictive.png', dpi = 144)
plt.show()
#%%

hmc.posterior_predictive(200)
#%%

gmm = sklearn.mixture.GMM(2)

gmm.fit(qs[:, :1])

print gmm.weights_
print gmm.means_
print gmm.covars_