Python scipy.optimize 模块，fmin_bfgs() 实例源码

def funktion(R):
    print('size of Matrix: '+str(R.shape[1])+' rows, '+str(R.shape[2])+' colums ')
    for l in range(R.shape[1]):
        print(str(l)+' of '+str(R.shape[1]))
        for n in range(R.shape[2]):
            Q=R[:,l,n]
            if isnan(sum(Q))==True:
                x=array([0.25,0.25,0.25])
                z=array([ NaN,  NaN,  NaN])
            else:
                x=array([0.25,0.25,0.25])
                z = opti.fmin_bfgs(my_cost, x, args = (Q,r),full_output=False,disp=False,retall=False,gtol=0.01)
                a_neu[:,l,n]=z
                f_orig[l,n]=dot((dot(r,x)-Q).transpose(),(dot(r,x)-Q))
                f_neu[l,n]=dot((dot(r,z)-Q).transpose(),(dot(r,z)-Q))
                residuals[l,n]=sum(z)
    return a_neu


####---snow1:

def bfgs(x0, f, f_prime, hessian=None):
    all_x_i = [x0[0]]
    all_y_i = [x0[1]]
    all_f_i = [f(x0)]
    def store(X):
        x, y = X
        all_x_i.append(x)
        all_y_i.append(y)
        all_f_i.append(f(X))
    optimize.fmin_bfgs(f, x0, f_prime, callback=store, gtol=1e-12)
    return all_x_i, all_y_i, all_f_i

def opt_ts(rslcs, method="weights", verbose=True):
    """
    I optimize the timeshifts between the rslcs to minimize the wtv between them.
    Note that even if the wtvdiff is only about two curves, we cannot split this into optimizing
    AB AC AD in a row, as this would never calculate BC, and BC is not contained into AB + AC.
    """
    rslcsc = [rs.copy() for rs in rslcs]

    # No need for reverse combis, as wtvdiff is symmetric.
    #couplelist = [couple for couple in [[rs1, rs2] for rs1 in rslcsc for rs2 in rslcsc] if couple[0] != couple[1]]

    indexes = np.arange(len(rslcsc))
    indlist = [c for c in [[i1, i2] for i1 in indexes for i2 in indexes] if c[1] > c[0]]
    couplelist = [[rslcsc[i1], rslcsc[i2]] for (i1, i2) in indlist]

    inishifts = [rs.timeshift for rs in rslcsc[1:]]

    def errorfct(timeshifts):
        if timeshifts.shape == ():
            timeshifts = np.array([timeshifts])
        for (rs, timeshift) in zip(rslcsc[1:], timeshifts):
            rs.timeshift = timeshift

        tvs = np.array([wtvdiff(rs1, rs2, method=method) for (rs1, rs2) in couplelist])
        ret = np.sum(tvs)
        print timeshifts, ret
        return ret

    if verbose:
        print "Starting time shift optimization ..."
        print "Initial pars : ", inishifts

    minout = spopt.fmin_powell(errorfct, inishifts, xtol=0.001, full_output=1, disp=verbose)
    #minout = spopt.fmin_bfgs(errorfct, inishifts, maxiter=None, full_output=1, disp=verbose, retall=0, callback=None)

    popt = minout[0]

    r2 = errorfct(popt) # This sets popt, and the optimal ML and source.

    if verbose:
        print "Optimal pars : ", popt

def f(R):
    for i in range(R.shape[1]):
        print(str(i)+' of '+str(R.shape[1]))
        Q=R[:,i]
        if isnan(sum(Q))==True:
            x=array([0.25,0.25,0.25])
            z=array([ NaN,  NaN, NaN])
        else:
            x=array([0.25,0.25,0.25])
            z = opti.fmin_bfgs(my_cost, x, args = (Q, r),full_output=False,disp=False,retall=False,gtol=0.01)
        Z_opti[:,i]=z
        f_orig[i]=dot((dot(r,x)-Q).transpose(),(dot(r,x)-Q))
        f_neu[i]=dot((dot(r,z)-Q).transpose(),(dot(r,z)-Q))
        residuals[i]=sum(z)
    return Z_opti

def f(R):
    for i in range(R.shape[1]):
        print(str(i)+' of '+str(R.shape[1]))
        Q=R[:,i]
        if isnan(sum(Q))==True:
            x=array([0.25,0.25,0.25,0.25])
            z=array([ NaN,  NaN, NaN, NaN])
        else:
            x=array([0.25,0.25,0.25,0.25])
            z = opti.fmin_bfgs(my_cost, x, args = (Q, r),full_output=False,disp=False,retall=False,gtol=0.01)
        Z_opti[:,i]=z
        f_orig[i]=dot((dot(r,x)-Q).transpose(),(dot(r,x)-Q))
        f_neu[i]=dot((dot(r,z)-Q).transpose(),(dot(r,z)-Q))
        residuals[i]=sum(z)
    return Z_opti

def test_spatial_median_2d():
    X = np.array([0., 0., 1., 1., 0., 1.]).reshape(3, 2)
    _, median = _spatial_median(X, max_iter=100, tol=1.e-6)

    def cost_func(y):
        dists = np.array([norm(x - y) for x in X])
        return np.sum(dists)

    # Check if median is solution of the Fermat-Weber location problem
    fermat_weber = fmin_bfgs(cost_func, median, disp=False)
    assert_array_almost_equal(median, fermat_weber)
    # Check when maximum iteration is exceeded a warning is emitted
    assert_warns(ConvergenceWarning, _spatial_median, X, max_iter=30, tol=0.)

def __call__(self, net, input, target):
        from scipy.optimize import fmin_bfgs
        if 'disp' not in self.kwargs:
            self.kwargs['disp'] = 0
        self.kwargs['maxiter'] = self.epochs

        x = fmin_bfgs(self.fcn, self.x.copy(), fprime=self.grad, callback=self.step,
                      **self.kwargs)
        self.x[:] = x

def optimiseparams(svdata, initialguessparams, trueparams, filterrunner, includejumps, randomstate):
    def objectivefunction(transformedparameterndarray):
        mockinfinity = 10000.
        print(transformedparameterndarray)
        parameterndarray = inversetransformparameterndarray(transformedparameterndarray, includejumps)
        params = sv.Params(*parameterndarray)
        print(params)
        if params.persistence >= 0.99 or params.persistence <= -0.99:
            print('Parameter out of bound: persistence = %f\n' % params.persistence)
            return mockinfinity
        if params.cor <= -0.99 or params.cor >= 0.99:
            print('Parameter out of bound: cor = %f\n' % params.cor)
            return mockinfinity
        if params.voloflogvar <= 0.01:
            print('Parameter out of bound: voloflogvar = %f\n' % params.vologlogvar)
            return mockinfinity
        if params.jumpintensity < 0. or params.jumpintensity > 1.:
            print('Parameter out of bound: jumpintensity = %f\n' % params.jumpintensity)
            return mockinfinity
        if params.jumpvol < 0.:
            print('Parameter out of bound: jumpvol = %f\n' % params.jumpvol)
            return mockinfinity
        loglikelihood = filterrunner(svdata, params, randomstate).stochfilter.loglikelihood
        loglikelihood = np.asscalar(loglikelihood)
        print('Loglikelihood: %f\n' % loglikelihood)
        return -loglikelihood

    print('True parameters:')
    print(trueparams)
    trueparamsloglikelihood = -objectivefunction(transformparameterndarray(np.array(trueparams), includejumps))
    print('True params loglikelihood: %f\n' % trueparamsloglikelihood)

    print('Initial guess parameters:')
    print(initialguessparams)

    print('Running the optimisation routine (BFGS)...')
    res = opt.fmin_bfgs(
            objectivefunction,
            x0=transformparameterndarray(np.array(initialguessparams), includejumps),
            epsilon=0.1111,
            disp=True
            )

    res = sv.Params(*inversetransformparameterndarray(res, includejumps))

    print('Result:')
    print(res)

    return res