#!/usr/bin/env python
################################################################################
#
#  CDDL HEADER START
#
#  The contents of this file are subject to the terms of the Common Development
#  and Distribution License Version 1.0 (the "License").
#
#  You can obtain a copy of the license at
#  http:# www.opensource.org/licenses/CDDL-1.0.  See the License for the
#  specific language governing permissions and limitations under the License.
#
#  When distributing Covered Code, include this CDDL HEADER in each file and
#  include the License file in a prominent location with the name LICENSE.CDDL.
#  If applicable, add the following below this CDDL HEADER, with the fields
#  enclosed by brackets "[]" replaced with your own identifying information:
#
#  Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
#
#  CDDL HEADER END
#
#  Copyright (c) 2017, Regents of the University of Minnesota.
#  All rights reserved.
#
#  Contributor(s):
#     Ellad B. Tadmor
#
################################################################################

# The docstring below is vc_description
'''Comparison of the analytical forces obtained from the model with forces
computed by numerical differentiation using Richardson extrapolation for a
randomly distorted non-periodic structure based off the equilibrium
face-centered cubic (fcc) structure. Separate configurations are tested for
each species supported by the model, as well as one containing a random
distribution of all species. Configurations used for testing are provided
as auxiliary files.'''

# TODO:
#
# This program does not deal with the issue of discontinuities in the
# model energy.  Such discontinuities could lead to incorrect numerical
# derivatives and to incorrect conclusions about the accuracy of the
# analytical forces.  A useful extension would be to check whether
# energy discontinuities are responsible for the component with the
# largest error, and if yes, ignore it and move on to the next largest,
# and so on.  Discontinuities could be found using the "detectedge"
# algorithm developed for the "Chebfun" Matlab package. See the
# discussion in "Piecewise Smooth Chebfuns", Pachon, Platte and Trefethen,
# https://www.cs.ox.ac.uk/files/717/NA-08-07.pdf
# and the code here:
# https://github.com/chebfun/chebfun/blob/development/%40fun/detectEdge.m

# Python 2-3 compatible code issues
from __future__ import print_function
try:
    input = raw_input
except NameError:
    pass

from ase.lattice.cubic import FaceCenteredCubic
from ase.calculators.kim.kim import KIM, KIM_get_supported_species_list
import kimvc
import sys
import random
import numpy as np
import numdifftools as nd
from numdifftools.step_generators import MaxStepGenerator
import math
import scipy.optimize

__version__ = "002"
__author__ = "Ellad Tadmor"

################################################################################
#
#   FUNCTIONS
#
################################################################################

################################################################################
def negpot(p, at=0, dof=0, atoms=None):
    '''
    Function that takes the value 'p' of degree of freedom 'dof' of atom 'at'
    and returns the negative of the total potential energy of full system of
    atoms. Used by the numerical derivative method.
    '''
    if atoms==None:
       return 0
    sve = (atoms[at].position)[dof]
    (atoms[at].position)[dof] = p
    pot = atoms.get_potential_energy()
    (atoms[at].position)[dof] = sve
    return -pot

################################################################################
def perform_numerical_derivative_check(vc, atoms, heading, dashwidth):
    '''
    Perform a numerical derivative check for the ASE atoms object in 'atoms'.
    '''
    # compute analytical forces (negative gradient of cohesive energy)
    energy = atoms.get_potential_energy()
    forces = atoms.get_forces()

    # Loop over atoms and compute numerical derivative check
    sg = [MaxStepGenerator(base_step=1e-4, num_steps=14,
                           use_exact_steps=True, step_ratio=1.6, offset=0),
          MaxStepGenerator(base_step=1e-3, num_steps=14,
                           use_exact_steps=True, step_ratio=1.6, offset=0),
          MaxStepGenerator(base_step=1e-2, num_steps=14,
                           use_exact_steps=True, step_ratio=1.6, offset=0)]
    nsg = 4 # number of step generators to try
    Dnegpot = [nd.Derivative(negpot, full_output=True)]
    for i in range(nsg-1):
        Dnegpot.append(nd.Derivative(negpot, step=sg[i], full_output=True))
    forces_num    = np.zeros(shape=(len(atoms),3), dtype=float, order='C')
    forces_uncert = np.zeros(shape=(len(atoms),3), dtype=float, order='C')
    forces_failed = np.zeros(shape=(len(atoms),3), dtype=int, order='C')
    for at in range(0,len(atoms)):
        for dof in range(0,3):
            p = (atoms[at].position)[dof]
            errmin_sg = 1e30
            failed_to_get_deriv = True
            for i in range(nsg):
                try:
                    val, info = Dnegpot[i](p,at=at,dof=dof,atoms=atoms)
                    if abs(val-forces[at,dof])<errmin_sg :
                        errmin_sg = abs(val-forces[at,dof])
                        val_best = val
                        info_best = info
                        failed_to_get_deriv = False
                except:
                    # Failed to compute derivative, so skip this value
                    (atoms[at].position)[dof] = p  # Restore value which may have been
                                                   # left changed when exception was
                                                   # generated.
            if failed_to_get_deriv :
                # Failed all attempts for this at/dof, assume this is
                # because the potential is being evaluated outside its
                # legal range and skip. (TODO: A more careful check
                # on the reason for the error would be good.)
                forces_failed[at,dof] = 1
            else:
                val = val_best
                info = info_best
                forces_num[at,dof] = val
                forces_uncert[at,dof] = info.error_estimate

    # Identify outliers using a box plot construction with fences
    # (See http://www.itl.nist.gov/div898/handbook/prc/section1/prc16.htm)
    # We'll take all results above the upper outer fence to be outliers.
    #
    # create flattened version of forces_uncert without failures
    numterms = 3*len(atoms) - np.sum(forces_failed)
    forces_uncert_without_failures = np.zeros(numterms, dtype=float)
    i = -1
    for at in range(0,len(atoms)):
        for dof in range(0,3):
            if not forces_failed[at,dof]:
                i += 1
                forces_uncert_without_failures[i] = forces_uncert[at,dof]
    uncert_lower_quartile  = np.percentile(forces_uncert_without_failures, 25)
    uncert_upper_quartile  = np.percentile(forces_uncert_without_failures, 75)
    uncert_interquartile_range = uncert_upper_quartile - uncert_lower_quartile
    uncert_upper_fence = uncert_upper_quartile + 3*uncert_interquartile_range

    # Initialize for printing
    frmt_head = \
                "{0:>6}  {1:>4} {2:>3} {3:>25} {4:>25} {5:>15} {6:>15}"
    frmt_line = \
        {True:  "{0: 6d}  {1:4} {2: 3d} {3: 25.15e} {4: 25.15e} {5: 15.5e} {6: 15.5e} {7:1}", \
         False: "             {2: 3d} {3: 25.15e} {4: 25.15e} {5: 15.5e} {6: 15.5e} {7:1}"}
    frmt_fail = "{0: 6d}  {1:4} {2: 3d} {3:>25} {4:>25} {5:>15} {6:>15} {7:1}"
    vc.rwrite('Comparison of analytical forces obtained from the model, '
              'the force computed as a numerical derivative')
    vc.rwrite('of the energy, the difference between them, and the uncertainty '
              'in the numerical estimate of the force.')
    vc.rwrite('The computed equilibrium lattice constant for the crystal (a0) '
              'is given in the heading.')
    vc.rwrite('')
    vc.rwrite(heading)
    vc.rwrite('-'*dashwidth)
    args = ('Part', 'Spec', 'Dir', 'Force_model', 'Force_numer', \
            '|Force diff|', 'uncertainty')
    vc.rwrite(frmt_head.format(*args))
    vc.rwrite('-'*dashwidth)

    # Identify max error and print numerical derivative results
    eps_prec = np.finfo(float).eps
    errmax = 0.
    at_least_one_result_discarded = False
    at_least_one_failure = False
    for at in range(0,len(atoms)):
        for dof in range(0,3):
            if forces_failed[at,dof]:
                # skip at/dof for which force was not computed
                args = (at+1, atoms[at].symbol, dof+1, "", "", "", "", "F")
                vc.rwrite(frmt_fail.format(*args))
                at_least_one_failure = True
                continue
            forcediff = abs(forces[at,dof]-forces_num[at,dof])
            den = max(abs(forces_num[at,dof]), eps_prec)
            if forces_uncert[at,dof] < uncert_upper_fence:
                # Result is not an outlier. Include it in determining max error
                lowacc_mark = ' '
                if forcediff/den > errmax:
                    errmax = forcediff/den
                    at_max  = at
                    dof_max = dof
            else:
                lowacc_mark = '*'
                at_least_one_result_discarded = True
            # Print results line
            args = (at+1, atoms[at].symbol, dof+1, \
                    forces[at,dof], forces_num[at,dof], \
                    forcediff, forces_uncert[at,dof], lowacc_mark)
            vc.rwrite(frmt_line[dof==0].format(*args))
            if (dof==2):
               vc.rwrite('-'*dashwidth)
    if at_least_one_result_discarded:
        vc.rwrite('* Starred lines are suspected outliers and are not '
                  'included when determining the error.')
        vc.rwrite('  A calculation is considered an outlier if it has an '
                  'uncertainty that lies at an abnormal')
        vc.rwrite('  distance from the other uncertainties in this set of '
                  'calculations.  Outliers are determined')
        vc.rwrite('  using the box plot construction with fences. '
                  'An outlier could indicate a problem with the')
        vc.rwrite('  the numerical differentiation or problems with the '
                  'potential energy, such as dicontinuities.')
    if at_least_one_failure:
        if at_least_one_result_discarded:
            vc.rwrite('')
        vc.rwrite('WARNING: Numerical derivative could not be computed for '
                  'at least one atom/dof.')
        vc.rwrite('         Failed calculations indicated with an "F". This '
                  'can be due to attempts')
        vc.rwrite('         to evaluate the model outside its legal range '
                  'during a numerical derivative.')
        vc.rwrite('         calculation. These lines are ignored when '
                  'computing the error.')

   # Print summary
    vc.rwrite('')
    vc.rwrite( \
        'Maximum error obtained for particle = {0:d}, direction = {1:d}:'. \
        format(at_max+1, dof_max+1))
    vc.rwrite('')
    vc.rwrite('              |F_model - F_numer|')
    vc.rwrite('    error = ----------------------- = {0:.5e}'. \
              format(errmax))
    vc.rwrite('              max{|F_numer|, eps}')
    vc.rwrite('')
    vc.rwrite('')

    return errmax

################################################################################
def cubic_cell_energy(alat, atoms, ncells_per_side):
    '''
    Calculate the energy of the passed 'atoms' structure containing a
    cubic structure with 'ncells_per_side'. Scale to lattice constant
    'alat' (passed as a nd array of length 1) and return the energy.
    '''
    acell = alat[0]*ncells_per_side
    atoms.set_cell([acell, acell, acell],scale_atoms=True)
    return atoms.get_potential_energy()

################################################################################
def do_vc(model, vc):
    '''
    Perform Numerical Derivative Check VC
    '''
    # Get supported species
    species = KIM_get_supported_species_list(model)
    species = kimvc.remove_species_not_supported_by_ASE(species)
    species.sort()

    # Basic cell parameters
    lattice_constant = 2.5
    max_lattice_constant = 10.0
    pert_amp = 0.1*lattice_constant
    ncells_per_side = 2
    seed = 13
    random.seed(seed)

    # Print Vc info
    dashwidth = 101
    vc.rwrite('')
    vc.rwrite('-'*dashwidth)
    vc.rwrite('Results for KIM Model      : %s' % model.strip())
    vc.rwrite('Supported species          : %s' % ' '.join(species))
    vc.rwrite('')
    vc.rwrite('random seed                = %d'    % seed)
    vc.rwrite('perturbation amplitude     = %0.3f' % pert_amp)
    vc.rwrite('number unit cells per side = %d'    % ncells_per_side)
    vc.rwrite('-'*dashwidth)
    vc.rwrite('')

    # Initialize variables
    errmaxmax = 0.
    got_atleast_one = False
    alat_ave = 0.

    # Perform numerical derivative check for monotatomic systems
    for spec in species:
        # find equilibrium lattice constant, so that the numerical derivatives
        # of all potentials are evaluated in a similar portion of their
        # potential energy surface, making comparisons between potentials
        # more menaningful.
        calc = KIM(model)
        alat = lattice_constant
        not_done = True
        got_equil = False
        while not_done:
            atoms = FaceCenteredCubic(size=(1, 1, 1),
                    latticeconstant=alat, symbol=spec, pbc=True)
            atoms.set_calculator(calc)
            try:
                res = scipy.optimize.minimize(cubic_cell_energy, alat,
                      args=(atoms,1), method='Nelder-Mead', tol=1e-6)
                alat = res.x[0]
                not_done = False
                got_equil = True
            except:
                # failed for some reason (assume it's because of KIM error)
                alat += 0.5
                if alat > max_lattice_constant:
                    not_done = False
        # Create a non-periodic crystal with the equilibrium lattice constant
        if not got_equil:
           vc.rwrite('FAILED to get equilibrium constant for species =  %s; skipping' % spec)
           continue
        alat_ave += alat
        atoms = FaceCenteredCubic(
                size=(ncells_per_side, ncells_per_side, ncells_per_side),
                latticeconstant=alat, symbol=spec, pbc=False)
        atoms.set_calculator(calc)
        # Rescale and perturb crystal for numerical derivative calculation
        kimvc.rescale_to_get_nonzero_forces(atoms, 0.01)
        kimvc.randomize_positions(atoms, pert_amp)
        kimvc.perturb_until_all_forces_sizeable(atoms, pert_amp)
        aux_file = 'config-'+spec+'.xyz'
        vc.vc_files.append(aux_file)
        vc.write_aux_ase_atoms(aux_file, atoms, 'xyz')
        # Perform numerical derivatice test
        a0string = "a0 = {}".format(alat)
        heading = 'MONOATOMIC STRUCTURE -- Species = ' + spec \
                + ' | ' + a0string \
                + ' | (Configuration in file "' + aux_file + '")'
        try:
            errmax = perform_numerical_derivative_check(vc, atoms, \
                                                        heading, dashwidth)
            errmaxmax = max(errmax, errmaxmax)
            got_atleast_one = True
        except:
            pass
    alat_ave = alat_ave/len(species)

    # Perform numerical derivative check for mixed system
    if len(species)>1:
        while True:
            atoms = FaceCenteredCubic(
                    size=(ncells_per_side, ncells_per_side, ncells_per_side),
                    latticeconstant=alat_ave, symbol="H", pbc=False)
            if len(atoms) < len(species):
                ncells_per_side += 1
            else:
                break
        kimvc.randomize_species(atoms, species)
        calc = KIM(model)
        atoms.set_calculator(calc)
        # Find equilibrium lattice constant
        res = scipy.optimize.minimize(cubic_cell_energy, alat_ave,
            args=(atoms,ncells_per_side), method='Nelder-Mead', tol=1e-6)
        alat = res.x[0]
        # Change periodicity to false
        atoms.set_pbc([False, False, False])
        # Rescale and perturb crystal for numerical derivative calculation
        kimvc.rescale_to_get_nonzero_forces(atoms, 0.01)
        kimvc.randomize_positions(atoms, pert_amp)
        kimvc.perturb_until_all_forces_sizeable(atoms, pert_amp)
        aux_file = 'config-'+''.join(species)+'.xyz'
        vc.vc_files.append(aux_file)
        vc.write_aux_ase_atoms(aux_file, atoms, 'xyz')
        # Perform numerical derivatice test
        a0string = "a0 = {}".format(alat)
        heading = 'MIXED STRUCTURE -- Species = ' + ' '.join(species) \
                + ' | ' + a0string \
                + ' | (Configuration in file "' + aux_file + '")'
        try:
            errmax = perform_numerical_derivative_check(vc, atoms, \
                                                        heading, dashwidth)
            errmaxmax = max(errmax, errmaxmax)
            got_atleast_one = True
        except:
            pass

    if got_atleast_one:

        # Compute grade
        vc.rwrite('='*dashwidth)
        vc.rwrite('Grade is based on maximum error across all systems '
                  '(for a model supporting multiple species):')
        vc.rwrite('')
        vc.rwrite('Maximum error     = {:15.5e}'.format(errmaxmax))
        vc.rwrite('Machine precision = {:15.5e}'.format(np.finfo(float).eps))
        vc.rwrite('')
        vc_grade, vc_comment \
            = kimvc.vc_letter_grade_machine_precision(errmaxmax)

        vc.rwrite('Note: Lower grades can reflect errors in the analytical '
                  'forces, a lack of smoothness in the potential, or '
                  'loss of precision in the numerical derivative computation.')
        vc.rwrite('')

        return vc_grade, vc_comment

    else:
        raise RuntimeError('Failed to compute a single numerical derivative.')
        return None, None

################################################################################
#
#   MAIN PROGRAM
#
###############################################################################
if __name__ == '__main__':

    vcargs = {"vc_name"        : "vc-forces-numerical-derivative",
              "vc_author"      : __author__,
              "vc_description" : kimvc.vc_stripall(__doc__),
              "vc_category"    : "consistency",
              "vc_grade_basis" : "graded",
              "vc_files"       : [],
              "vc_debug"       : False}  # Set to True to get exception traceback info

    # Get the model extended KIM ID:
    model = input("Model Extended KIM ID = ")

    # Execute VC
    kimvc.setup_and_run_vc(do_vc, model, **vcargs)