Replace joblib to cPickle, change min_descriptor_len & train_test_split ratio, refactor training script and folder layout, request only non-disordered

Author: blokhin

mpds_ml_labs/prediction.py

@@ -1,11 +1,10 @@
 
 from __future__ import division
 import os
+import cPickle
 
 import numpy as np
 
-from sklearn.externals import joblib
-
 
 human_names = {
     'z': {
@@ -17,7 +16,7 @@
     'y': {
         'name': 'enthalpy of formation',
         'units': 'kJ g-at.-1',
-        'symbol': 'H',
+        'symbol': 'ΔH',
         'rounding': 0
     },
     'x': {
@@ -26,6 +25,12 @@
         'symbol': 'C<sub>p</sub>',
         'rounding': 0
     },
+    #'w': {
+    #    'name': 'band gap for direct transition',
+    #    'units': 'eV',
+    #    'symbol': 'e<sub>dir.</sub>',
+    #    'rounding': 1
+    #},
     'k': {
         'name': 'Seebeck coefficient',
         'units': 'muV K-1',
@@ -62,21 +67,24 @@
 periodic_numbers_normed = [(i - pmin)/(pmax - pmin) for i in periodic_numbers]
 
 
-def get_descriptor(ase_obj, kappa=18, overreach=False):
+def get_descriptor(ase_obj, kappa=None, overreach=False):
     """
     From ASE object obtain
     a vectorized atomic structure
     populated to a certain fixed (relatively big) volume
     defined by kappa
     """
+    if not kappa: kappa = 18
     if overreach: kappa *= 2
 
     norms = np.array([ np.linalg.norm(vec) for vec in ase_obj.get_cell() ])
     multiple = np.ceil(kappa / norms).astype(int)
     ase_obj = ase_obj.repeat(multiple)
     com = ase_obj.get_center_of_mass() # NB use recent ase version here, because of the new element symbols
     ase_obj.translate(-com)
-    del ase_obj[[atom.index for atom in ase_obj if np.sqrt(np.dot(atom.position, atom.position)) > kappa]]
+    del ase_obj[
+        [atom.index for atom in ase_obj if np.sqrt(np.dot(atom.position, atom.position)) > kappa]
+    ]
 
     ase_obj.center()
     ase_obj.set_pbc((False, False, False))
@@ -99,6 +107,10 @@ def get_descriptor(ase_obj, kappa=18, overreach=False):
 def load_ml_model(prop_model_files):
     ml_model = {}
     for n, file_name in enumerate(prop_model_files, start=1):
+        if not os.path.exists(file_name):
+            print("No file %s" % file_name)
+            continue
+
         basename = file_name.split(os.sep)[-1]
         if basename.startswith('ml') and basename[3:4] == '_' and basename[2:3] in human_names:
             prop_id = basename[2:3]
@@ -107,10 +119,11 @@ def load_ml_model(prop_model_files):
             prop_id = str(n)
             print("No property name detected in file %s" % basename)
 
-        model = joblib.load(file_name)
-        if hasattr(model, 'predict') and hasattr(model, 'metadata'):
-            ml_model[prop_id] = model
-            print("Model metadata: %s" % model.metadata)
+        with open(file_name, 'rb') as f:
+            model = cPickle.load(f)
+            if hasattr(model, 'predict') and hasattr(model, 'metadata'):
+                ml_model[prop_id] = model
+                print("Model metadata: %s" % model.metadata)
 
     print("Loaded property models: %s" % len(ml_model))
     return ml_model
@@ -134,9 +147,7 @@ def ase_to_ml_model(ase_obj, ml_model):
     d_dim = len(descriptor)
 
     if not ml_model: # testing
-
-        test_prop = round(np.sum(descriptor))
-        return {prop_id: {'value': test_prop, 'mae': 0, 'r2': 0} for prop_id in human_names.keys()}, None
+        return {prop_id: {'value': 42, 'mae': 0, 'r2': 0} for prop_id in human_names.keys()}, None
 
     for prop_id, regr in ml_model.items(): # production
 

mpds_ml_labs/ml_mpds.py renamed to train_model.py

@@ -33,7 +33,7 @@ def get_regr(a=None, b=None):
     )
 
 
-def estimate_quality(algo, args, values, attempts=30, nsamples=0.4):
+def estimate_quality(algo, args, values, attempts=30, nsamples=0.33):
     results = []
     for _ in range(attempts):
         X_train, X_test, y_train, y_test = train_test_split(args, values, test_size=nsamples)
@@ -52,13 +52,12 @@ def estimate_quality(algo, args, values, attempts=30, nsamples=0.4):
     return avg_mae, avg_r2
 
 
-def mpds_get_data(prop_id):
+def mpds_get_data(prop_id, descriptor_kappa):
     """
-    NB
-    currently pressure is not taken into account,
-    however must be
+    Fetch, massage, and save dataframe from the MPDS
+    NB currently pressure is not taken into account!
     """
-    print("Getting %s" % human_names[prop_id]['name'])
+    print("Getting %s with descriptor kappa = %s" % (human_names[prop_id]['name'], descriptor_kappa))
     starttime = time.time()
 
     client = MPDSDataRetrieval()
@@ -84,12 +83,14 @@ def mpds_get_data(prop_id):
     # these should be corrected by LPF editors soon
     if prop_id == 'z':
         props = props[props['Value'] < 2000]
-    elif prop_id == 'w':
-        props = props[(props['Value'] > 0) & (props['Value'] < 20)]
+    #elif prop_id == 'w': # NB this requires additional treatment for zero band gaps
+    #    props = props[(props['Value'] > 0) & (props['Value'] < 20)]
     elif prop_id == 'u':
         props = props[props['Value'] > 0]
 
-    to_drop = props[(props['Cname'] == 'Temperature') & (props['Cunits'] == 'K') & ((props['Cvalue'] < 200) | (props['Cvalue'] > 400))]
+    to_drop = props[
+        (props['Cname'] == 'Temperature') & (props['Cunits'] == 'K') & ((props['Cvalue'] < 200) | (props['Cvalue'] > 400))
+    ]
 
     print("Rows to drop by criteria: %s" % len(to_drop))
     props.drop(to_drop.index, inplace=True)
@@ -100,22 +101,26 @@ def mpds_get_data(prop_id):
 
     print("Got %s distinct crystalline phases" % len(phases))
 
-    min_descriptor_len = 220
+    min_descriptor_len = 200
     max_descriptor_len = min_descriptor_len*10
     data_by_phases = {}
 
+    print("Computing descriptors...")
     pbar = ProgressBar()
     for item in pbar(client.get_data(
-        {"props": "atomic structure"},
+        {
+            "props": "atomic structure",
+            "classes": "non-disordered"
+        },
         fields={'S':['phase_id', 'entry', 'chemical_formula', 'cell_abc', 'sg_n', 'setting', 'basis_noneq', 'els_noneq']},
         phases=phases
     )):
         crystal = MPDSDataRetrieval.compile_crystal(item, 'ase')
         if not crystal: continue
-        descriptor = get_descriptor(crystal)
+        descriptor = get_descriptor(crystal, kappa=descriptor_kappa)
 
         if len(descriptor) < min_descriptor_len:
-            descriptor = get_descriptor(crystal, overreach=True)
+            descriptor = get_descriptor(crystal, kappa=descriptor_kappa, overreach=True)
             if len(descriptor) < min_descriptor_len:
                 continue
 
@@ -147,79 +152,89 @@ def mpds_get_data(prop_id):
 
     print("Done %s rows in %1.2f sc" % (len(struct_props), time.time() - starttime))
 
-    export_file = MPDSExport.save_df(struct_props, prop_id)
-    print("Saving %s" % export_file)
+    struct_props.export_file = MPDSExport.save_df(struct_props, prop_id)
+    print("Saving %s" % struct_props.export_file)
 
     return struct_props
 
 
+def tune_model(data_file):
+    """
+    Load saved data and perform simple regressor parameter tuning
+    """
+    basename = data_file.split(os.sep)[-1]
+    if basename.startswith('df') and basename[3:4] == '_' and basename[2:3] in human_names:
+        tag = basename[2:3]
+        print("Detected property %s" % human_names[tag]['name'])
+    else:
+        tag = None
+        print("No property name detected")
+
+    df = pd.read_pickle(data_file)
+
+    X = np.array(df['Descriptor'].tolist())
+    y = df['Avgvalue'].tolist()
+
+    results = []
+    for parameter_a in range(20, 501, 20):
+        avg_mae, avg_r2 = estimate_quality(get_regr(a=parameter_a), X, y)
+        results.append([parameter_a, avg_mae, avg_r2])
+        print("%s\t\t\t%s\t\t\t%s" % (parameter_a, avg_mae, avg_r2))
+    results.sort(key=lambda x: (-x[1], x[2]))
+
+    print("Best result:", results[-1])
+    parameter_a = results[-1][0]
+
+    results = []
+    for parameter_b in range(1, 13):
+        avg_mae, avg_r2 = estimate_quality(get_regr(a=parameter_a, b=parameter_b), X, y)
+        results.append([parameter_b, avg_mae, avg_r2])
+        print("%s\t\t\t%s\t\t\t%s" % (parameter_b, avg_mae, avg_r2))
+    results.sort(key=lambda x: (-x[1], x[2]))
+
+    print("Best result:", results[-1])
+    parameter_b = results[-1][0]
+
+    print("a = %s b = %s" % (parameter_a, parameter_b))
+
+    regr = get_regr(a=parameter_a, b=parameter_b)
+    regr.fit(X, y)
+    regr.metadata = {'mae': avg_mae, 'r2': round(avg_r2, 2)}
+
+    if tag:
+        export_file = MPDSExport.save_model(regr, tag)
+        print("Saving %s" % export_file)
+
+
 if __name__ == "__main__":
     try:
         arg = sys.argv[1]
     except IndexError:
         sys.exit(
-            "What to do?\n"
-            "Please, provide either a *prop_id* letter (%s) for a property data to be downloaded,\n"
-            "or a data *filename* generated after a property data download." % ", ".join(human_names.keys())
+    "What to do?\n"
+    "Please, provide either a *prop_id* letter (%s) for a property data to be downloaded and fitted,\n"
+    "or a data *filename* for tuning the model." % ", ".join(human_names.keys())
         )
+    try:
+        descriptor_kappa = int(sys.argv[2])
+    except:
+        descriptor_kappa = None
 
     if arg in human_names.keys():
 
-        # getting the data from scratch by prop_id
-        struct_props = mpds_get_data(arg)
+        struct_props = mpds_get_data(arg, descriptor_kappa)
 
         X = np.array(struct_props['Descriptor'].tolist())
         y = struct_props['Avgvalue'].tolist()
 
         avg_mae, avg_r2 = estimate_quality(get_regr(), X, y)
+
         print("Avg. MAE: %.2f" % avg_mae)
         print("Avg. R2 score: %.2f" % avg_r2)
 
-    elif os.path.exists(arg):
-
-        # loading saved data
-        basename = arg.split(os.sep)[-1]
-        if basename.startswith('df') and basename[3:4] == '_' and basename[2:3] in human_names:
-            tag = basename[2:3]
-            print("Detected property %s" % human_names[tag]['name'])
-        else:
-            tag = None
-            print("No property name detected")
-
-        df = pd.read_pickle(arg)
-
-        X = np.array(df['Descriptor'].tolist())
-        y = df['Avgvalue'].tolist()
+        tune_model(struct_props.export_file)
 
-        # simple regressor parameter tuning
-        results = []
-        for a in range(60, 501, 20):
-            avg_mae, avg_r2 = estimate_quality(get_regr(a=a), X, y)
-            results.append([a, avg_mae, avg_r2])
-            print("%s\t\t\t%s\t\t\t%s" % (a, avg_mae, avg_r2))
-        results.sort(key=lambda x: (-x[1], x[2]))
-
-        print("Best result:", results[-1])
-        a = results[-1][0]
-
-        results = []
-        for b in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]:
-            avg_mae, avg_r2 = estimate_quality(get_regr(a=a, b=b), X, y)
-            results.append([b, avg_mae, avg_r2])
-            print("%s\t\t\t%s\t\t\t%s" % (b, avg_mae, avg_r2))
-        results.sort(key=lambda x: (-x[1], x[2]))
-
-        print("Best result:", results[-1])
-        b = results[-1][0]
-
-        print("a = %s b = %s" % (a, b))
-
-        regr = get_regr(a=a, b=b)
-        regr.fit(X, y)
-        regr.metadata = {'mae': avg_mae, 'r2': round(avg_r2, 2)}
-
-        if tag:
-            export_file = MPDSExport.save_model(regr, tag)
-            print("Saving %s" % export_file)
+    elif os.path.exists(arg):
+        tune_model(arg)
 
     else: raise RuntimeError("Unrecognized argument: %s" % arg)