Python NLP.plot方法代码示例

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
def milpTransfer(originProb):
    newProb = NLP(originProb.f, originProb.x0)
    originProb.fill(newProb)
    newProb.discreteVars = originProb.discreteVars
    def err(s): # to prevent text output
        raise OpenOptException(s)
    newProb.err = err
    for fn in ['df', 'd2f', 'c', 'dc', 'h', 'dh']:
        if hasattr(originProb, fn) and getattr(originProb.userProvided, fn) or originProb.isFDmodel:
            setattr(newProb, fn, getattr(originProb, fn))
    
    newProb.plot = 0
    newProb.iprint = -1
    newProb.nlpSolver = originProb.nlpSolver 
    return newProb

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
    def run(self, plot=True):
        """
        Solves the optimization problem.
        """        
        # Initial try
        p0 = self.get_p0()
        
        #Lower bounds and Upper bounds (HARDCODED FOR QUADTANK)
        lbound = N.array([0.0001]*len(p0))
        if self.gridsize == 1:
            ubound = [10.0]*(self.gridsize*self.nbr_us)
        else:
            ubound = [10.0]*(self.gridsize*self.nbr_us) + [0.20,0.20,0.20,0.20,N.inf]*((self.gridsize-1))

        
        #UPPER BOUND FOR VDP
        #ubound = [0.75]*(self.gridsize*self.nbr_us)+[N.inf]*((self.gridsize-1)*self.nbr_ys)
        
        if self.verbosity >= Multiple_Shooting.NORMAL:
            print 'Initial parameter vector: '
            print p0
            print 'Lower bound:', len(lbound)
            print 'Upper bound:', len(ubound)

        # Get OpenOPT handler
        p_solve = NLP(self.f,p0,lb = lbound, ub=ubound,maxFunEvals = self.maxFeval, maxIter = self.maxIter, ftol=self.ftol, maxTime=self.maxTime)
        
        #If multiple shooting is preformed or single shooting
        if self.gridsize > 1:
            p_solve.h  = self.h
        
        if plot:
            p_solve.plot = 1

        self.opt = p_solve.solve(self.optMethod)        
        
        return self.opt

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
    def fit_node(self,index):
        qnode=self.qlist[index]
        print qnode.q
        th=qnode.th_condensed['a3']
        counts=qnode.th_condensed['counts']
        counts_err=qnode.th_condensed['counts_err']
        print qnode.th_condensed['counts'].std()
        print qnode.th_condensed['counts'].mean()
        maxval=qnode.th_condensed['counts'].max()
        minval=qnode.th_condensed['counts'].min()
        diff=qnode.th_condensed['counts'].max()-qnode.th_condensed['counts'].min()\
            -qnode.th_condensed['counts'].mean()
        sig=qnode.th_condensed['counts'].std()

        if diff-2*sig>0:
            #the difference between the high and low point and
            #the mean is greater than 3 sigma so we have a signal
            p0=findpeak(th,counts,1)
            print 'p0',p0
            #Area center width Bak
            center=p0[0]
            width=p0[1]
            sigma=width/2/N.sqrt(2*N.log(2))
            Imax=maxval-minval
            area=Imax*(N.sqrt(2*pi)*sigma)
            print 'Imax',Imax
            pin=[area,center,width,0]





            if 1:
                p = NLP(chisq, pin, maxIter = 1e3, maxFunEvals = 1e5)
                #p.lb=lowerm
                #p.ub=upperm
                p.args.f=(th,counts,counts_err)
                p.plot = 0
                p.iprint = 1
                p.contol = 1e-5#3 # required constraints tolerance, default for NLP is 1e-6

    # for ALGENCAN solver gradtol is the only one stop criterium connected to openopt
    # (except maxfun, maxiter)
    # Note that in ALGENCAN gradtol means norm of projected gradient of  the Augmented Lagrangian
    # so it should be something like 1e-3...1e-5
                p.gradtol = 1e-5#5 # gradient stop criterium (default for NLP is 1e-6)
        #print 'maxiter', p.maxiter
        #print 'maxfun', p.maxfun
                p.maxIter=50
    #    p.maxfun=100

        #p.df_iter = 50
                p.maxTime = 4000
        #r=p.solve('scipy_cobyla')
            #r=p.solve('scipy_lbfgsb')
                #r = p.solve('algencan')
                print 'ralg'
                r = p.solve('ralg')
                print 'done'
                pfit=r.xf
                print 'pfit openopt',pfit
                print 'r dict', r.__dict__

            if 0:
                print 'curvefit'
                print sys.executable
                pfit,popt=curve_fit(gauss2, th, counts, p0=pfit, sigma=counts_err)
                print 'p,popt', pfit,popt
                perror=N.sqrt(N.diag(popt))
                print 'perror',perror
                chisqr=chisq(pfit,th,counts,counts_err)
                dof=len(th)-len(pfit)
                print 'chisq',chisqr
            if 0:
                oparam=scipy.odr.Model(gauss)
                mydatao=scipy.odr.RealData(th,counts,sx=None,sy=counts_err)
                myodr = scipy.odr.ODR(mydatao, oparam, beta0=pfit)
                myoutput=myodr.run()
                myoutput.pprint()
                pfit=myoutput.beta
            if 1:
                print 'mpfit'
                p0=pfit
                parbase={'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}
                parinfo=[]
                for i in range(len(p0)):
                    parinfo.append(copy.deepcopy(parbase))
                for i in range(len(p0)):
                    parinfo[i]['value']=p0[i]
                fa = {'x':th, 'y':counts, 'err':counts_err}
                #parinfo[1]['fixed']=1
                #parinfo[2]['fixed']=1
                m = mpfit(myfunct_res, p0, parinfo=parinfo,functkw=fa)
                if (m.status <= 0):
                    print 'error message = ', m.errmsg
                params=m.params
                pfit=params
                perror=m.perror
                #chisqr=(myfunct_res(m.params, x=th, y=counts, err=counts_err)[1]**2).sum()
                chisqr=chisq(pfit,th,counts,counts_err)
#.........这里部分代码省略.........

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
        # 2nd arg - start point
        #p = NLP(f, x0, c=c,  gtol=gtol, contol=contol, iprint = 50, maxIter = 10000, maxFunEvals = 1e7, name = 'NLP_1')
    x0 = rd.uniform(-0.8, 0.8, num_dof*(2*N+1))
    # x0[2*N] = 0
    # x0[2*(2*N)] = 0
    # x0[3*(2*N)] = 0
    #x0 = rd.uniform(-0.1, 0.1, 2*N+1).tolist() + rd.uniform(-0.01, 0.01, 2*N+1).tolist() + rd.uniform(-.001, 0.001, 2*N+1).tolist()

    #x0 = range(0, num_dof*(2*N+1))
    #p = NLP(f, x0, c=c, h=h, gtol=gtol, contol=contol, iprint = 1, maxIter = 1000, maxFunEvals = 1e7, name = 'NLP_1')
    p = NLP(f, x0, c=c, h=h, gtol=gtol, contol=contol, iprint = 1, maxIter = 700, maxFunEvals = 1e7, name = 'NLP_1')

    #p = NLP(f, x0,  gtol=gtol, contol=contol, iprint = 50, maxIter = 10000, maxFunEvals = 1e7, name = 'NLP_1')

        #optional: graphic output, requires pylab (matplotlib)
    p.plot = True

    solver = 'ralg'
    #solver = 'scipy_cobyla'
        #solver = 'algencan'
        #solver = 'ipopt'
        #solver = 'scipy_slsqp'

        # solve the problem

    r = p.solve(solver, plot=0) # string argument is solver name


        # r.xf and r.ff are optim point and optim objFun value
        # r.ff should be something like 132.05
        #print r.xf, r.ff

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
        #lowerm[0:3]=[-1,-1,-1]
        upperm=N.ones(len(p0))
    if 1:
        p = NLP(Entropy, p0, maxIter = 1e3, maxFunEvals = 1e5)

    if 0:
        p = NLP(chisq, p0, maxIter = 1e3, maxFunEvals = 1e5)


    if 0:
        p = NLP(max_wrap, p0, maxIter = 1e3, maxFunEvals = 1e5)
    if 0:
        p.lb=lowerm
        p.ub=upperm
        p.args.f=(h,k,l,fq,fqerr,x,z,cosmat_list,coslist,flist)
        p.plot = 0
        p.iprint = 1
        p.contol = 1e-5#3 # required constraints tolerance, default for NLP is 1e-6

    # for ALGENCAN solver gradtol is the only one stop criterium connected to openopt
    # (except maxfun, maxiter)
    # Note that in ALGENCAN gradtol means norm of projected gradient of  the Augmented Lagrangian
    # so it should be something like 1e-3...1e-5
        p.gradtol = 1e-5#5 # gradient stop criterium (default for NLP is 1e-6)
        #print 'maxiter', p.maxiter
        #print 'maxfun', p.maxfun
        p.maxIter=50
    #    p.maxfun=100

        #p.df_iter = 50
        p.maxTime = 4000

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
    p.maxiter=10
#    p.maxfun=100
    
# see also: help(NLP) -> maxTime, maxCPUTime, ftol and xtol
# that are connected to / used in lincher and some other solvers

    # optional: check of user-supplied derivatives
    #p.checkdf()
    #p.checkdc()
    #p.checkdh()

# last but not least:
# please don't forget,
# Python indexing starts from ZERO!!

    p.plot = 0
    p.iprint = 1
    #p.df_iter = 50
    p.maxTime = 4000
    print 'solving'
    r = p.solve('algencan')
    print 'done'
    pout=r.xf
    print 'solution:', pout

    print len(pout)
    
    P_up,P_down=transform_p(pout,Mx,Mz,M)
    P=P_up-P_down
    
    

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
solvers = [
    'ralg',
    'lincher',
    'gsubg',
    'scipy_slsqp',
    'scipy_cobyla',
    'interalg',
    'auglag',
    'ptn',
    'mma'
    ]


pro = NLP(f=objLegrand.cost,x0=theta,df=objLegrand.sensitivity,lb=lb,ub=ub,A=A,b=b)
pro.plot = True
rRalg = pro.solve('ralg')

pro = NLP(f=objLegrand.cost,x0=theta3,df=objLegrand.sensitivity,lb=lb,ub=ub,A=A,b=b)
pro.plot = True
rRalg = pro.solve('ralg')



pro = NLP(f=objLegrand.cost,x0=theta,df=objLegrand.sensitivity,lb=lb,ub=ub,A=A,b=b)
pro.plot = True
rSLSQP = pro.solve('scipy_slsqp')

# we try to refine the solution

pro = NLP(f=objLegrand.cost,x0=rPswarm.xf,df=objLegrand.sensitivity,lb=lb,ub=ub,A=A,b=b)

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
def fitpeak(x,y,yerr):
    maxval=x.max()
    minval=x.min()
    diff=y.max()-y.min()-y.mean()
    sig=y.std()
    print 'diff',diff,'std',sig
    if diff-1*sig>0:
        #the difference between the high and low point and
        #the mean is greater than 3 sigma so we have a signal
        p0=findpeak(x,y,2)
        print 'p0',p0
        #Area center width Bak area2 center2 width2
        center1=p0[0]
        width1=p0[1]
        center2=p0[2]
        width2=p0[3]
        sigma=width/2/N.sqrt(2*N.log(2))
        ymax=maxval-minval
        area=ymax*(N.sqrt(2*pi)*sigma)
        print 'ymax',ymax
        pin=[area,center1,width1,0,area,center2,width2]





        if 1:
            p = NLP(chisq, pin, maxIter = 1e3, maxFunEvals = 1e5)
            #p.lb=lowerm
            #p.ub=upperm
            p.args.f=(x,y,yerr)
            p.plot = 0
            p.iprint = 1
            p.contol = 1e-5#3 # required constraints tolerance, default for NLP is 1e-6

# for ALGENCAN solver gradtol is the only one stop criterium connected to openopt
# (except maxfun, maxiter)
# Note that in ALGENCAN gradtol means norm of projected gradient of  the Augmented Lagrangian
# so it should be something like 1e-3...1e-5
            p.gradtol = 1e-5#5 # gradient stop criterium (default for NLP is 1e-6)
    #print 'maxiter', p.maxiter
    #print 'maxfun', p.maxfun
            p.maxIter=50
#    p.maxfun=100

    #p.df_iter = 50
            p.maxTime = 4000
    #r=p.solve('scipy_cobyla')
        #r=p.solve('scipy_lbfgsb')
            #r = p.solve('algencan')
            print 'ralg'
            r = p.solve('ralg')
            print 'done'
            pfit=r.xf
            print 'pfit openopt',pfit
            print 'r dict', r.__dict__
        if 1: 
            print 'mpfit'
            p0=pfit
            parbase={'value':0., 'fixed':0, 'limited':[0,0], 'limits':[0.,0.]}
            parinfo=[]
            for i in range(len(p0)):
                parinfo.append(copy.deepcopy(parbase))
            for i in range(len(p0)): 
                parinfo[i]['value']=p0[i]
            fa = {'x':x, 'y':y, 'err':yerr}
            #parinfo[1]['fixed']=1
            #parinfo[2]['fixed']=1
            m = mpfit(myfunct_res, p0, parinfo=parinfo,functkw=fa)
            if (m.status <= 0): 
                print 'error message = ', m.errmsg
            params=m.params
            pfit=params
            perror=m.perror
            #chisqr=(myfunct_res(m.params, x=th, y=counts, err=counts_err)[1]**2).sum()
            chisqr=chisq(pfit,x,y,yerr)
            dof=m.dof
            #Icalc=gauss(pfit,th)
            #print 'mpfit chisqr', chisqr
        ycalc=gauss(pfit,x)

        if 1:
            width_x=N.linspace(p0[0]-p0[1],p0[0]+p0[1],100)
            width_y=N.ones(width_x.shape)*(maxval-minval)/2
            pos_y=N.linspace(minval,maxval,100)
            pos_x=N.ones(pos_y.shape)*p0[0]
            if 0:
                
                pylab.errorbar(th,counts,counts_err,marker='s',linestyle='None',mfc='black',mec='black',ecolor='black')
                pylab.plot(width_x,width_y)
                pylab.plot(pos_x,pos_y)
                pylab.plot(x,ycalc)
                pylab.show()

    else:
        #fix center
        #fix width
        print 'no peak'
        #Area center width Bak
        area=0
#.........这里部分代码省略.........

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
solvers = [
    'ralg',
    'lincher',
    'gsubg',
    'scipy_slsqp',
    'scipy_cobyla',
    'interalg',
    'auglag',
    'ptn',
    'mma'
    ]


pro = NLP(f=objLegrand.cost,x0=theta,df=objLegrand.adjointInterpolate1,lb=lb,ub=ub,A=A,b=b)
pro.plot = True
rRalg = pro.solve('ralg',xtol=1e-20)

pro = NLP(f=objLegrand.cost,x0=theta,df=objLegrand.sensitivity,lb=lb,ub=ub,A=A,b=b)
pro.plot = True
rRalg = pro.solve('lincher')

pro = NLP(f=objLegrand.cost,x0=theta,df=objLegrand.sensitivity,lb=lb,ub=ub,A=A,b=b)
pro.plot = True
rSLSQP = pro.solve('scipy_slsqp')

# we try to refine the solution

pro = NLP(f=objLegrand.cost,x0=rPswarm.xf,df=objLegrand.sensitivity,lb=lb,ub=ub,A=A,b=b)
pro.plot = True
rSLSQP = pro.solve('ralg')

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]
p.ub[3*N/4:] = inf

#p.ub[4] = 4
#p.lb[5], p.ub[5] = 8, 15

#A = zeros((K, N))
#b = zeros(K)
#for i in xrange(K):
#    A[i] = 1+cos(i+arange(N))
#    b[i] = sin(i)
#p.A = A
#p.b = b

#p.Aeq = zeros(p.n)
#p.Aeq[100:102] = 1
#p.beq = 11

p.contol = 1e-6
p.plot = 1
p.maxFunEvals = 1e7
p.name = 'nlp_4'
p.debug=1
solver = 'ralg'
solver = 'gsubg'
#solver = 'scipy_cobyla'
solver = 'algencan'
#solver = 'ipopt'

r = p.solve(solver, xlabel = 'time', fTol = 10, debug=0, maxIter = 5500, plot=0, maxTime=1000, ftol = 1e-8, xtol = 1e-6, iprint=1, showLS=0, showFeas=0, show_hs=0)

# 需要导入模块: from openopt import NLP [as 别名]
# 或者: from openopt.NLP import plot [as 别名]

#.........这里部分代码省略.........
        def eqZ(i): 
            def equation(x):
                px, pz = x[spx:spz], x[spz:spf]
                return px[i] - pz[i]
            return equation

        def eqpz(j, J):      
            def equation(x):
                F, Z = Table.unflatten(index=parameter.factors, columns=parameter.industries, sam=x[sF:sZ]), x[sZ:spx]
                return Z[j] - self.calibration.b[J] * np.prod([F[J][H] ** self.calibration.beta[J][H] for H in parameter.factors])
            return equation

        def eqpf(H):
            def equation(x):
                F = Table.unflatten(index=parameter.factors, columns=parameter.industries, sam=x[sF:sZ])
                return sum(F[J][H] for J in parameter.industries) - self.calibration.FF[H]
            return equation

        def eqF(h, H, j, J):      
            def equation(x):
                F, Z, pz, pf = Table.unflatten(index=parameter.factors, columns=parameter.industries, sam=x[sF:sZ]), x[sZ:spx], x[spz:spf], x[spf:epf]
                pf = np.array([float(x[spf:epf]), 1])
                return F[J][H] - self.calibration.beta[J][H] * pz[j] * Z[j] / pf[h]
            return equation
        self.calibration = copy(calibration)        
        self.calibration.X0 = self.calibration.X0['HH']
        j = i = len(parameter.industries)
        h = len(parameter.factors)
        sX = 0
        sF = sX + i
        sZ = sF + i * h
        spx = sZ + j
        spz = spx + i
        spf = spz + j
        epf = spf + h - 1
        epf_numerair = epf + 1
        lb = np.array([(np.float64(0.001))] * epf)
        self.x = x = np.empty(epf, dtype='f64')
        x[sX:sF] = self.calibration.X0.data
        x[sF:sZ] = self.calibration.F0.as_matrix().flatten()
        x[sZ:spx] = self.calibration.Z0.data
        x[spx:spz] = [1] * i
        x[spz:spf] = [1] * j
        x[spf:epf] = [1] * (epf - spf)       
        self.t = x[:]
        if debug:
            self.x = x = np.array([21.1] * epf)
        print x
        xnames = [] * (epf_numerair)
        xnames[sX:sF] = self.calibration.X0.names
        xnames[sF:sZ] = [i+h+' ' for i in parameter.industries for h in parameter.factors]
        xnames[sZ:spx] = self.calibration.Z0.names
        xnames[spx:spz] = parameter.industries
        xnames[spz:spf] = parameter.industries
        xnames[spf:epf] = parameter.factors
        xnames[epf] = parameter.factors[-1]

        xtypes = [] * epf_numerair
        xtypes[sX:sF] = ['X0'] * len(self.calibration.X0)
        xtypes[sF:sZ] = ['F'] * (len(parameter.industries) + len(parameter.factors))
        xtypes[sZ:spx] = ['F0'] * len(self.calibration.F0)
        xtypes[spx:spz] = ['pb'] * len(parameter.industries)
        xtypes[spz:spf] = ['pz'] * len(parameter.industries)
        xtypes[spf:epf] = ['pf'] * len(parameter.factors)
        xtypes[epf] = ['pf']

        self.xnametypes = ['%s %s' % (xnames[i], xtypes[i]) for i in range(epf - 1)]

        constraints = []
        for i, I in enumerate(parameter.industries):    
            industry, Industry = unlink2(i, I)
            constraints.append(eqX(industry, Industry)) 
            constraints.append(eqpx(industry))
            constraints.append(eqZ(industry))
            constraints.append(eqpz(industry, Industry))
                    
            

        for F in parameter.factors:
            Factor = unlink(F)
            constraints.append(eqpf(Factor))

        for f, F in enumerate(parameter.factors):
            for i, I in enumerate(parameter.industries):
                industry, Industry = unlink2(i, I)
                factor, Factor = unlink2(f, F)
                constraints.append(eqF(factor, Factor, industry, Industry))


        self.UU = UU = lambda x: - np.prod([x[i] ** self.calibration.alpha[i] for i in range(len(parameter.industries))])



        p = NLP(UU, x, h=constraints, lb=lb, iprint = 50, maxIter = 10000, maxFunEvals = 1e7, name = 'NLP_1')        
        p.plot = debug
        self.r = p.solve('ralg', plot=0)

        if debug: 
            for i, constraint in enumerate(constraints):
                print(i, '%02f' % constraint(self.r.xf))