chapman.py 文件源码

def chapman_fit(alt,
                ne,
                x0=[1e6, 300, 50],
                bounds=[(1, None),
                        (150, 500),
                        (30, 80)],
                verbose=False,
                **kwds):
    """
    """
    # optimization setup
    z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
    chapman = chapman_sym(z, Nm, Hm, H_O)
    dNm = SYM.diff(chapman, Nm)
    dHm = SYM.diff(chapman, Hm)
    dH_O = SYM.diff(chapman, H_O)
    chapman_f = SYM.lambdify((z, Nm, Hm, H_O),
                             chapman,
                             modules='numexpr')
    dNm_f = SYM.lambdify((z, Nm, Hm, H_O),
                         dNm,
                         modules='numexpr')
    dHm_f = SYM.lambdify((z, Nm, Hm, H_O),
                         dHm,
                         modules='numexpr')
    dH_O_f = SYM.lambdify((z, Nm, Hm, H_O),
                          dH_O,
                          modules='numexpr')
    # define cost function
    y = NP.asarray(ne)
    def J(x):
        Nm, Hm, H_O = x
        if verbose:
            print('-' * 80)
            print(x)
        y_hat = NP.array([chapman_f(z, Nm, Hm, H_O) for z in alt])
        diff = y - y_hat
        J1 = NP.array([dNm_f(z, Nm, Hm, H_O) for z in alt])
        J2 = NP.array([dHm_f(z, Nm, Hm, H_O) for z in alt])
        J3 = NP.array([dH_O_f(z, Nm, Hm, H_O) for z in alt])
        return (NP.dot(diff, diff),
                NP.array([-2 * NP.sum(diff * J1),
                          -2 * NP.sum(diff * J2),
                          -2 * NP.sum(diff * J3)]))
    # minimize cost function
    x_star, f, d = scipy.optimize.fmin_l_bfgs_b(J,
                                                x0,
                                                bounds=bounds,
                                                **kwds)
    assert d['warnflag'] == 0
    return x_star