mrqcof

gistfile1.txt

    MACHEPS = np.finfo(np.float64).eps

    def estimate_mrqcof(self, xk, *args):
        '''
        Estimates the gradient and hessian matrix for the fit
        (differential_evolution can't do this).
        '''
        temp_parameters = np.copy(self.parameters)

        derivmatrix = np.zeros((len(xk), self.numpoints), np.float64)
        ei = np.zeros((len(xk),), np.float64)
        epsilon = (pow(MACHEPS, 1. / 3)
                   * np.fmax(np.fabs(xk), 0.1))

        alpha = np.zeros((len(xk), len(xk)), np.float64)
        beta = np.zeros(len(xk), np.float64)

        #this is wasteful of function evaluations. Requires 2 evaluations for
        #LHS and RHS
        for k in range(len(xk)):
            temp_parameters[self.fitted_parameters] = xk
            d = epsilon[k]

            rcost = self.energy(xk[k] + d, *args)
            lcost = self.energy(xk[k] - d, *args)
            beta[k] = (rcost - lcost) / 2. / d

            temp_parameters[self.fitted_parameters[k]] = xk[k] + d
            f2 = self.model(temp_parameters, *args)
            temp_parameters[self.fitted_parameters[k]] = xk[k] - d
            f1 = self.model(temp_parameters, *args)

            derivmatrix[k, :] = (f2 - f1) / 2. / d

        for i in range(len(xk)):
            for j in range(i + 1):
                val = np.sum(derivmatrix[i] * derivmatrix[j] / self.edata**2)
                alpha[i, j] = val;
                alpha[j, i] = val;

        return alpha, beta