positive polynomial optimization + taylor model sketch idea

euler.py

import cvxpy as cvx
import numpy as np
import sos
import sympy as sy
import matplotlib.pyplot as plt
#raised chebyehve
t = sy.symbols("t")

N = 5
# seems like even N becomes infeasible.
terms = [sy.chebyshevt(n,t) + 1 for n in range(N)] # raised chebyshev functions are positive on interval [-1,1]
print(terms)

'''
for i in range(1,4):
    ts = np.linspace(-1,1,100)
    #print(ts)
    #rint(sy.lambdify(t,terms[i], 'numpy')(ts))
    plt.plot( ts , sy.lambdify(t,terms[i])(ts))
plt.show()
'''

vdict = {}

l,d = sos.polyvar(terms) # lower bound on solution
vdict.update(d)

w,d = sos.polyvar(terms, nonneg=True) # width of tube. Width is always positive (nonneg)
vdict.update(d)

u = l + w # upper curve is higher than lower by width

def picard(t,f):
    return sy.integrate(f, [t,-1,t]) + np.exp(-1) # picard integration on [-1,1] interval with initial cond x(-1)=1/e

ui = picard(t,u)
li = picard(t,l)

c = []

def split(y , N): # split a polynomial into lower an upper parts.
    yp = sy.poly(y, gens=t)
    lower = sum([ c*t**p for (p,), c in  yp.terms() if p < N]) 
    #upper = sum([ c*x**p for (p,), c in  yp.terms() if p > N])
    upper = y - lower
    return lower,upper

terms = [sy.chebyshevt(n,t) + 1 for n in range(N+1)]

# ui <= u
lowerui, upperui = split(ui, N) # need to truncate  highest power of u using interval method
print(lowerui)
print(upperui)
du = upperui.subs(t,1)  #Is this where the even dependence of N comes from?
#c += [ du >= sos.cvxify(upperui.subs(t,1), vdict), du >= sos.cvxify(upperui.subs(t,1)]  # , upperui.subs(t,-1))
print(du)
lam1,d = sos.polyvar(terms,nonneg=True) # positive polynomial
vdict.update(d)
# This makes the iterated interval inside the original interval
c += sos.poly_eq(  lowerui + du + lam1 , u , vdict) # write polynomial inequalities in slack equality form


# l <= li 
# 
lam2, d = sos.polyvar(terms,nonneg=True)
vdict.update(d)
c += sos.poly_eq(  l + lam2 , li , vdict) # makes new lower bound higher than original lower bound

obj = cvx.Minimize( sos.cvxify(w.subs(t ,0.9), vdict) ) # randomly picked reasonable objective. Try minimax?
#obj = cvx.Maximize( sos.cvxify(l.subs(t ,1), vdict) )
print(c)
prob = cvx.Problem(obj, c)

res = prob.solve(verbose=True) #solver=cvx.CBC
print(res)
lower = sy.lambdify(t, sos.poly_value(l , vdict))
upper = sy.lambdify(t, sos.poly_value(u , vdict))

#plt.plot(ts, upper(ts) - np.exp(ts) ) # plot differences
#plt.plot(ts, lower(ts) - np.exp(ts)  )
ts = np.linspace(-1,1,100)
plt.plot(ts, upper(ts) , label= "upper")
plt.plot(ts, lower(ts) , label= "lower")
plt.plot(ts, np.exp(ts) , label= "exact")
#plt.plot(ts,np.exp(ts) - lower(ts) )
plt.legend()
plt.show()

''' 
if I need to add in 
interval rounding to get closure
is there a point to this? Is it actually simpler in any sense?

Collecting up chebyshev compoentns and chebysehv splitting would perform 
lanczos economization. That'd be coo

What about a bvp

Get iterative formulation.

And what are the requirements
1. We need an iterative contractive operator
2. We need to confirm all functions that forall t, l <= f <= u 
map to in between li and ui. This part might be challenging
3. Get the interval contracting and small.

x <= a
y = Lx
Lx <= La ? Yes, if positive semi definite. Otherwise we need to split it.
No. Nice try. not component wise inequality.

Secondary question: finitely confirming a differential operator is positive semi definite
forall x, xLx >= 0 ?
Similar to the above. Make regions in space.





Value function learning is contractive.
hmm.


Piecewise lyapunov functions



Being able to use an LP makes it WAY faster, WAY more stable, and opens up sweet MIPpurtunities.


'''