QML Hyperparameter fitter

fit_qml.py

#!/usr/bin/env python2
import sys
import os

import copy 
import qml
from qml.fchl import get_atomic_kernels_fchl
from qml.fchl import get_atomic_symmetric_kernels_fchl
from qml.fchl import generate_fchl_representation
import random
from time import time
import numpy as np
from copy import deepcopy
from scipy.optimize import minimize
from fml.math import cho_solve
np.set_printoptions(linewidth=999999999)


CUT_DISTANCE = 1e6

class CostFunction(object):

    def __init__(self):

        print "Initializing"
        self.Y = np.load("Y.npy")
        self.X = np.load("X.npy")

        self.n_cross = 10

    def calc_maes(self, Ki, Yi):

        diffs = []

        total = len(Yi)
        test = total // 5
        train = total - test

        random.seed(667)

        maes = []

        for i in range(self.n_cross):

            split = range(total)

            random.shuffle(split)

            training_index = split[:train]
            test_index = split[-test:]

            Y  = Yi[training_index]
            Ys = Yi[test_index]

            C = deepcopy(Ki[0][training_index][:,training_index])
            C[np.diag_indices_from(C)] += 10.0**(-7.0)

            alpha = cho_solve(C, Y)

            Yss = np.dot(Ki[0][test_index][:,training_index], alpha)

            diff = (Yss  - Ys)
            mae = np.mean(np.abs(diff))
            maes.append(mae)


        maes = np.array(maes)
        mae = np.mean(maes)
        stddev = np.std(maes) / np.sqrt(self.n_cross)

        return mae, stddev

    def get_mae(self, parameters):

        R_WIDTH = parameters[0]
        C_WIDTH = parameters[1]
        D_WIDTH = parameters[2]
        T_WIDTH = parameters[3]
        SCALE_DISTANCE = parameters[4]
        SCALE_ANGULAR  = parameters[5]
        SIGMA = 25.0
        SIGMA = parameters[6]

        print "p = np." + repr(parameters)

        sigmas = [SIGMA]

        start_K = time()

        K = get_atomic_symmetric_kernels_fchl(self.X, sigmas, cut_distance=CUT_DISTANCE, 
        t_width=T_WIDTH, d_width=D_WIDTH,  r_width=R_WIDTH, c_width=C_WIDTH, 
        order=1, scale_angular=SCALE_ANGULAR, scale_distance=SCALE_DISTANCE)

        k_time = time() - start_K


        start_cv = time()
        mae, stddev = self.calc_maes(K, self.Y)
        cv_time = time() - start_cv

        print "Cost: %10.4f +/- %10.4f kcal/mol     time = %10.4f" % (mae, stddev, (k_time+cv_time))

        return mae

if __name__ == "__main__":

    # p = np.array([ 0.01, 0.01, 0.15,  np.pi/1.0,  0.1,  0.10, 25.0])
    p = np.array([ 1e-5, 1e-5, 0.05,  np.pi/1.0,  0.1,  0.10, 25.0])

    cost = CostFunction()

    minimize(cost.get_mae, p, method="Nelder-Mead",
        options={"maxiter": 1000, "disp": True})