How to use the forcebalance.finite_difference.fdwrap function in forcebalance

def test_fd_stencils(self):
        """Check finite difference stencils return approximately correct results"""
        func = lambda x: x[0]**2
        fd_stencils = [function for function in dir(forcebalance.finite_difference) if re.match('^f..?d.p$',function)]

        self.logger.debug("Comparing fd stencils against some simple functions\n")
        for func in self.functions:
            for p in range(1):
                for x in range(10):
                    input = [0,0,0]
                    input[p]=x
                    f=forcebalance.finite_difference.fdwrap(func[0], input, p)
                    for stencil in fd_stencils:
                        fd = eval("forcebalance.finite_difference.%s" % stencil)
                        result = fd(f,.0001)
                        if re.match('^f..d.p$', stencil):
                            self.assertAlmostEqual(result[0], func[1](input,p), places=3)
                            self.assertAlmostEqual(result[1], func[2](input,p), places=3)
                        else:
                            self.assertAlmostEqual(result, func[1](input,p), places=3)

self.FF.make(mvals_)
            moments = self.engine.multipole_moments(polarizability='polarizability' in self.ref_moments, optimize=self.optimize_geometry)
            # Unpack from dictionary.
            return self.unpack_moments(moments)

        self.FF.make(mvals)
        ref_momvals = self.unpack_moments(self.ref_moments)
        calc_moments = self.engine.multipole_moments(polarizability='polarizability' in self.ref_moments, optimize=self.optimize_geometry)
        calc_momvals = self.unpack_moments(calc_moments)

        D = calc_momvals - ref_momvals
        dV = np.zeros((self.FF.np,len(calc_momvals)))

        if AGrad or AHess:
            for p in self.pgrad:
                dV[p,:], _ = f12d3p(fdwrap(get_momvals, mvals, p), h = self.h, f0 = calc_momvals)
                
        Answer['X'] = np.dot(D,D)
        for p in self.pgrad:
            Answer['G'][p] = 2*np.dot(D, dV[p,:])
            for q in self.pgrad:
                Answer['H'][p,q] = 2*np.dot(dV[p,:], dV[q,:])

        if not in_fd():
            self.calc_moments = calc_moments
            self.objective = Answer['X']

        return Answer

   @param[in] xyzs List of OpenMM positions
   @param[in] settings OpenMM settings for creating the System
   @param[in] boxes Periodic box vectors
   @return G First derivative of the energies in a N_param x N_coord array
   @return Hd Second derivative of the energies (i.e. diagonal Hessian elements) in a N_param x N_coord array

   """

   G        = np.zeros((FF.np,len(xyzs)))
   Hd       = np.zeros((FF.np,len(xyzs)))
   if not AGrad:
      return G, Hd
   E0       = energy_driver(mvals, pdb, FF, xyzs, settings, boxes)
   for i in range(FF.np):
      # Not doing the three-point finite difference anymore.
      G[i,:] = f1d2p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,boxes=boxes),h,f0=E0)
   return G, Hd

Derivative of the energy in a FF.np x length array.
    
    """
    G = np.zeros((FF.np, length))
    
    if not AGrad:
        return G
    def energy_driver(mvals_):
        FF.make(mvals_)
        return engine.energy()

    ED0 = energy_driver(mvals)
        
    for i in pgrad:
        logger.info("%i %s\r" % (i, (FF.plist[i] + " "*30)))
        EDG, _   = f12d3p(fdwrap(energy_driver, mvals, i), h, f0=ED0)

        G[i,:]   = EDG[:]
    return G

logger.info("%s\n" % (FF.plist[i] + " "*30))
        ED1      = fdwrap(energy_driver,mvals,i)(h)
        E1       = ED1[:,0]
        D1       = ED1[:,1:]
        b = np.exp(-(E1-E0)/kT)
        b /= np.sum(b)
        if 'h_' in property_kwargs:
            property_kwargs['h_'] = H0.copy() + (E1-E0)
        if 'd_' in property_kwargs:
            property_kwargs['d_'] = D1.copy()
        S = -1*np.dot(b,np.log(b))
        InfoContent = np.exp(S)
        if InfoContent / len(E0) < 0.1:
            logger.warn("Warning: Effective number of snapshots: % .1f (out of %i)\n" % (InfoContent, len(E0)))
        P1 = property_driver(b=b,**property_kwargs)
        EDM1      = fdwrap(energy_driver,mvals,i)(-h)
        EM1       = EDM1[:,0]
        DM1       = EDM1[:,1:]
        b = np.exp(-(EM1-E0)/kT)
        b /= np.sum(b)
        if 'h_' in property_kwargs:
            property_kwargs['h_'] = H0.copy() + (EM1-E0)
        if 'd_' in property_kwargs:
            property_kwargs['d_'] = DM1.copy()
        S = -1*np.dot(b,np.log(b))
        InfoContent = np.exp(S)
        if InfoContent / len(E0) < 0.1:
            logger.warn("Warning: Effective number of snapshots: % .1f (out of %i)\n" % (InfoContent, len(E0)))
        PM1 = property_driver(b=b,**property_kwargs)
        G[i] = (P1-PM1)/(2*h)
    if 'h_' in property_kwargs:
        property_kwargs['h_'] = H0.copy()

def FDCheckG(self):
        """ Finite-difference checker for the objective function gradient.

        For each element in the gradient, use a five-point finite difference
        stencil to compute a finite-difference derivative, and compare it to
        the analytic result.

        """

        Adata        = self.Objective.Full(self.mvals0,Order=1)['G']
        Fdata        = np.zeros(self.FF.np,dtype=float)
        printcool("Checking first derivatives by finite difference!\n%-8s%-20s%13s%13s%13s%13s" \
                  % ("Index", "Parameter ID","Analytic","Numerical","Difference","Fractional"),bold=1,color=5)
        for i in range(self.FF.np):
            Fdata[i] = f1d7p(fdwrap(self.Objective.Full,self.mvals0,i,'X',Order=0),self.h)
            Denom = max(abs(Adata[i]),abs(Fdata[i]))
            Denom = Denom > 1e-8 and Denom or 1e-8
            D = Adata[i] - Fdata[i]
            Q = (Adata[i] - Fdata[i])/Denom
            cD = abs(D) > 0.5 and "\x1b[1;91m" or (abs(D) > 1e-2 and "\x1b[91m" or (abs(D) > 1e-5 and "\x1b[93m" or "\x1b[92m"))
            cQ = abs(Q) > 0.5 and "\x1b[1;91m" or (abs(Q) > 1e-2 and "\x1b[91m" or (abs(Q) > 1e-5 and "\x1b[93m" or "\x1b[92m"))
            logger.info("\r    %-8i%-20s% 13.4e% 13.4e%s% 13.4e%s% 13.4e\x1b[0m\n" \
                  % (i, self.FF.plist[i][:20], Adata[i], Fdata[i], cD, D, cQ, Q))

logger.info("\rNow working on" + str(d) + 50*' ' + '\r')
            odir = os.path.join(os.getcwd(),d)
            #if os.path.exists(odir):
            #    shutil.rmtree(odir)
            if not os.path.exists(odir): os.makedirs(odir)
            apath = os.path.join(odir, "current")
            submit_psi(apath, d, mvals)
            for p in range(self.FF.np):
                def subjob(mvals_,h):
                    apath = os.path.join(odir, str(p), str(h))
                    submit_psi(apath, d, mvals_)
                    #logger.info("Will set up a job for %s, parameter %i\n" % (d, p))
                    return 0.0
                if self.callderivs[d][p]:
                    if AHess:
                        f12d3p(fdwrap(subjob, mvals, p, h=self.h), h = self.h, f0 = 0.0)
                    elif AGrad:
                        if self.bidirect:
                            f12d3p(fdwrap(subjob, mvals, p, h=self.h), h = self.h, f0 = 0.0)
                        else:
                            f1d2p(fdwrap(subjob, mvals, p, h=self.h), h = self.h, f0 = 0.0)

def func1(arg):
                mvals = list(mvals0)
                mvals[pidxj] += arg
                return f1d5p(fdwrap(self.Objective.Full,mvals,pidxi,'X',Order=0),self.h)
            return func1

    @param[in] settings OpenMM settings for creating the System
    @param[in] boxes Periodic box vectors
    @return G First derivative of the energies in a N_param x N_coord array

    """

    G        = np.zeros((FF.np,len(xyzs)))
    if not AGrad:
        return G
    E0       = energy_driver(mvals, pdb, FF, xyzs, settings, simulation, boxes)
    CheckFDPts = False
    for i in range(FF.np):
        G[i,:], _ = f12d3p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h,f0=E0)
        if CheckFDPts:
            # Check whether the number of finite difference points is sufficient.  Forward difference still gives residual error of a few percent.
            G1 = f1d7p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h)
            dG = G1 - G[i,:]
            dGfrac = (G1 - G[i,:]) / G[i,:]
            print "Parameter %3i 7-pt vs. central derivative : RMS, Max error (fractional) = % .4e % .4e (% .4e % .4e)" % (i, np.sqrt(np.mean(dG**2)), max(np.abs(dG)), np.sqrt(np.mean(dGfrac**2)), max(np.abs(dGfrac)))
    return G

def func1(arg):
                mvals = list(mvals0)
                mvals[pidxj] += arg
                return f1d5p(fdwrap(self.Objective.Full,mvals,pidxi,'X',Order=0),self.h)
            return func1