import os, getpass
os.environ['MPDS_KEY'] = getpass.getpass()

!pip install git+https://github.com/mpds-io/pytheos.git#egg=pytheos

import numpy as np
from pytheos import BM3Model

v0, k0, k0p = 10, 200, 4

exp_bm3 = BM3Model()
v_data = v0 * np.linspace(0.6, 1, 20)

fit_params = exp_bm3.make_params(v0=v0, k0=k0, k0p=k0p)
p_bm3 = exp_bm3.eval(fit_params, v=v_data)
p_data = p_bm3 + np.random.normal(0.0, 2, 20)

fit_params['v0'].vary = False
fitresult_bm3 = exp_bm3.fit(p_data, fit_params, v=v_data, weights=None)

print(fitresult_bm3.fit_report())

!pip install mpds_client>=0.0.17

import warnings
import pandas as pd
import numpy as np
from numpy.linalg import det
from ase.geometry import cellpar_to_cell

from mpds_client import MPDSDataRetrieval

client = MPDSDataRetrieval()

dfrm_k0p = client.get_dataframe({"classes": "binary", "elements": "O", "props": "pressure derivative of isothermal bulk modulus"})
dfrm_k0p = dfrm_k0p[np.isfinite(dfrm_k0p['Phase'])] # only data for the existing distinct phases
avg_k0p =  dfrm_k0p.groupby('Phase')['Value'].median().to_frame().reset_index().rename(columns={'Value': 'avg_k0p'})

dfrm_k0 = client.get_dataframe({"props": "isothermal bulk modulus"}, phases=set(dfrm_k0p['Phase'].tolist()))
avg_k0 =  dfrm_k0.groupby('Phase')['Value'].median().to_frame().reset_index().rename(columns={'Value': 'avg_k0'})
avg_k0 =  avg_k0.merge(avg_k0p, how='inner', on='Phase')

def get_cell_volume(a, b, c, alpha, beta, gamma):
    '''
    Calculate V from cell parameters.
    NB ab_normal and a_direction are standard.
    '''
    return abs(
        det(
            cellpar_to_cell([a, b, c, alpha, beta, gamma])
        )
    )

pvts = {}

for matrix in client.get_data(
    {"props": "cell parameters - pressure diagram"},
    phases=set(dfrm_k0p['Phase'].tolist()), # only those phases we have experimental bulk modulus
    fields={}                               # all fields
):
    try: decks = matrix['sample']['measurement'][0]['property']['matrix']
    except KeyError:
        warnings.warn('Error: there is no expected property in the entry %s' % matrix)
        continue

    p_all, v_all, t_all = [], [], []
    for deck in decks:
        p, t, v = deck[0], deck[1], get_cell_volume(*deck[2:])
        p_all.append(p)
        v_all.append(v)
        t_all.append(t)

    # TODO here in principle we should do something smarter
    # than just omitting the data for the same phase
    if matrix['sample']['material']['phase_id'] in pvts and len(pvts[ matrix['sample']['material']['phase_id'] ][0]) > len(p_all):
            warnings.warn('Skipping entry %s' % matrix['sample']['material']['entry'])
            continue

    pvts[ matrix['sample']['material']['phase_id'] ] = [matrix['sample']['material']['entry'], p_all, v_all, t_all]

pvts = [[key] + value for key, value in pvts.items()]
pvts = pd.DataFrame(pvts, columns=['Phase', 'Entry', 'P', 'V', 'T'])
print(pvts)

dfrm = avg_k0.merge(pvts, how='inner', on='Phase')

for n, system in dfrm.iterrows():

    params = exp_bm3.make_params(v0=system['V'][0], k0=system['avg_k0'], k0p=system['avg_k0p'])

    fitresult_bm3 = exp_bm3.fit(system['P'], params, v=system['V'], weights=None)
    if abs(fitresult_bm3.params['k0'].value - system['avg_k0']) > 50:
        # show the discrepancies, if any
        print("*" * 30 + " Distinct phase https://mpds.io/#phase_id/%s " % system['Phase'] + "*" * 30)
        print("BM_fit: %4.1f 	 BM_exp: %4.1f" % (fitresult_bm3.params['k0'].value, system['avg_k0']))
        print("BM0p_fit: %1.2f 	 BM0p_exp: %1.2f" % (fitresult_bm3.params['k0p'].value, system['avg_k0p']))