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

def btb(cls, pred_tr, pred_te, labels_tr, labels_te, data_tr, data_te, verbose=0):

        train_preds = np.clip(pred_tr, -0.99, 8.99)
        test_preds = np.clip(pred_te, -0.99, 8.99)
        num_classes = 8

        # train offsets
        offsets = np.ones(num_classes) * -0.5
        offset_train_preds = np.vstack((train_preds, train_preds, labels_tr))
        for j in range(num_classes):
            def train_offset(x):
                return -apply_offset(offset_train_preds, x, j)

            offsets[j] = fmin_powell(train_offset, offsets[j])

        # apply offsets to test
        data = np.vstack((test_preds, test_preds))
        for j in range(num_classes):
            data[1, data[0].astype(int) == j] = data[0, data[0].astype(int) == j] + offsets[j]

        pred_class = np.round(np.clip(data[1], 1, 8)).astype(int)
        return pred_class

def btb2(cls, pred_tr, pred_te, labels_tr, labels_te, data_tr, data_te, verbose=0):
        train_preds = np.clip(pred_tr, -0.99, 8.99)
        test_preds = np.clip(pred_te, -0.99, 8.99)
        num_classes = 8

        # train offsets
        offsets = np.array([0.1, -1, -2, -1, -0.8, 0.02, 0.8, 1])
        # offsets = np.ones(num_classes) * -0.5
        offset_train_preds = np.vstack((train_preds, train_preds, labels_tr))
        for j in range(num_classes):
            def train_offset(x):
                return -apply_offset(offset_train_preds, x, j)

            offsets[j] = fmin_powell(train_offset, offsets[j])

        # apply offsets to test
        data = np.vstack((test_preds, test_preds))
        for j in range(num_classes):
            data[1, data[0].astype(int) == j] = data[0, data[0].astype(int) == j] + offsets[j]

        pred_class = np.round(np.clip(data[1], 1, 8)).astype(int)
        return pred_class

def testMaxent(self):
        k = 2
        data = csb.io.load(self.data_fn)
        model = MaxentModel(k)
        model.sample_weights()
        posterior = MaxentPosterior(model, data[:100000] / 180. * numpy.pi)

        model.get() * 1.

        x0 = posterior.model.get().flatten()
        target = lambda w:-posterior(w, n=50)
        x = fmin_powell(target, x0, disp=False)

        self.assertTrue(x != None)
        self.assertTrue(len(x) == k * k * 4)

        posterior.model.set(x)
        posterior.model.normalize(True)

        xx = numpy.linspace(0 , 2 * numpy.pi, 500)
        fx = posterior.model.log_prob(xx, xx)

        self.assertAlmostEqual(posterior.model.log_z(integration='simpson'),
                               posterior.model.log_z(integration='trapezoidal'),
                               places=2)

        self.assertTrue(fx != None)
        z = numpy.exp(log_sum_exp(numpy.ravel(fx))) 
        self.assertAlmostEqual(z * xx[1] ** 2, 1., places=1)

def powell(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_powell(f, x0, callback=store, ftol=1e-12)
    return all_x_i, all_y_i, all_f_i

def bm_sigsq_optim(tree,traits,rate=1):
    # tree_utils2.assign_sigsq(tree,[rate])
    start = [rate]
    opt = optimize.fmin_powell(calc_like_sigsq, start, args = (tree,traits),full_output=False,disp=True)
    return [tree.get_newick_repr(True),opt]

def opt_magshift(lcs, rawdispersionmethod, verbose=True, trace=False):
    """
    I optimize the magnitude shifts between lightcurves. This is usually a first thing to do.
    First curve does not get shifted.

    You might think that this can be done by just using the median mag of a curve, but here we use the dispersion method,
    to get it right even in case of curves with non-overlapping regions.
    """

    couplelist = [couple for couple in [[lc1, lc2] for lc1 in lcs for lc2 in lcs] if couple[0] != couple[1]]

    def setshifts(p): # p are the absolute magshifts.
        for (l, shift) in zip(lcs[1:], p):
            l.magshift = shift
            # We don't use the relative shiftmag() !
        if trace:
            pycs.gen.util.trace(lcs)

    inipars = np.array([l.magshift for l in lcs[1:]])
    if verbose:
        print "Starting magshift optimization ..."
        print "Initial pars : ", inipars

    def errorfct(p):
        setshifts(p)
        d2values = np.array([rawdispersionmethod(*couple)["d2"] for couple in couplelist])
        return np.mean(d2values)

    minout = spopt.fmin_powell(errorfct, inipars, full_output=1, disp=verbose)
    popt = minout[0]
    if popt.shape == ():
        popt = np.array([popt])

    if verbose:
        print "Optimal pars : ", popt
    setshifts(popt)
    # We do not return anything

def optcflat(self, verbose = False):
        """
        Optimizes only the "border coeffs" so to get zero slope at the extrema
        Run optc() first ...
        This has to be done with an iterative optimizer
        """

        full = self.getc(m=1)
        inip = self.getc(m=1)[[0, 1, -2, -1]] # 4 coeffs

        def setp(p):
            full[[0, 1, -2, -1]] = p
            self.setcflat(full)

        if verbose:
            print "Starting flat coeff optimization ..."
            print "Initial pars : ", inip

        def errorfct(p):
            setp(p)
            return self.r2(nostab=False) # To get the same as optc would return !

        minout = spopt.fmin_powell(errorfct, inip, full_output=1, disp=verbose)
        popt = minout[0]
        if popt.shape == ():
            popt = np.array([popt])

        if verbose:
            print "Optimal pars : ", popt
        setp(popt)
        return self.r2(nostab=False) # We include the stab points, like optc does.
        # This last line also updates self.lastr2 ...

def optcmltv(lcs, spline, verbose=True):
    """
    I will optimize the coefficients of the spline so to minimize the mltv.
    I do not use the microlensing of the lcs at all !

    Simple powell optimization, slow. A pity.

    Add BOK and time shifts in there and it might be bingo !

    Would be more efficient if we add knots on the fly
    """


    inic = spline.getc(m=2)

    def setc(c):
        spline.setc(c, m=2)

    def errorfct(c):
        setc(c)
        (tv, dist) = mltv(lcs, spline, weight=False)
        print "put weight"
        return tv + 0.1*spline.tv()

    minout = spopt.fmin_powell(errorfct, inic, full_output=1, disp=verbose)
    copt = minout[0]

    # We find a common shift to all coeffs so that the level matches

    meanc = np.mean(spline.getc(m=2))
    meanmag = np.mean(np.concatenate([l.getmags(noml = True) for l in lcs]))

    setc(copt)
    spline.c += meanmag - meanc

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 btb_rank(cls, pred_tr, pred_te, labels_tr, labels_te, data_tr, data_te, verbose=0):

        pred_tr = pred_tr.argsort().argsort() / float(len(pred_tr))
        pred_te = pred_te.argsort().argsort() / float(len(pred_te))

        init = np.zeros(7)
        for i in range(7):
            init[i] = np.where(labels_tr <= (i + 1), 1.0, 0.0).sum() / float(len(labels_tr))
        print init

        init2 = np.array(init)

        for j in range(7):

            def lossfunc(x):
                ary = np.array(init2)
                ary[j] = x
                _loss = -1.0 * model_util.evaluate_class(labels_tr, cls.to_class(pred_tr, ary))
                if verbose == 1:
                    print _loss
                return _loss

            init2[j] = fmin_powell(lossfunc, init2[j])

        pred_clss = cls.to_class(pred_te, init2)
        return pred_clss

def opt_fluxshift(lcs, sourcespline, verbose = True, trace=False):
    """
    Optimizes the flux shift and the magshift of the lcs (not the first one)
    to get the best fit to the "sourcespline". Does not touch the microlensing, nor the spline.
    So this is a building block to be used iteratively with the other optimizers.
    Especially of the sourcespline, as we fit here even on regions not well constrained by the spline !

    The spline should typically well fit to the first curve.
    """


    for l in lcs[1:]: # We don't touch the first one.

        if verbose:
            print "Fluxshift optimization on %s ..." % (l)

        minfs = l.getminfluxshift()
        maxfs = 100000.0
        inip = l.fluxshift


        inip = (0, 0)
        # p = (fluxshift, magshift)
        def setp(p):
            fs = p[0]*1000.0
            ms = p[1]*0.1

            if fs < minfs:
                l.setfluxshift(minfs, consmag = False)
            else:
                l.setfluxshift(fs, consmag = False)

            l.magshift = ms


        def errorfct(p):    
            setp(p)
            return pycs.gen.spl.r2([l], sourcespline)

        minout = spopt.fmin_powell(errorfct, inip, full_output=1, xtol=0.001, disp=verbose)         

        popt = minout[0]    
        setp(popt)
        if verbose:
            print "Done with %s ..." % (l)
    # We do not return anything