# Test points 
xx=c(0,0.05,0.236,0.382,0.5,0.618,0.786,0.886,1)

x=xx[c(1,3,4,6,9)]
y=c(1,1.25,2.5,5,10)

plot(x,y)

# If p(x)=a+bx+x^2
# Let A=[x_i^(j-1)] and u=(a,b,...,f,g). We need to minimize |p(x_i)-y_i|.
# Let us constraining by min r, with r>=Au-y and r>=y-Au, p(x)>=1, p(x)<=10.
u0=c(5,0,0,6)  #p(x)=1 is an initial guess.
m=length(x)
n=length(u0)-1

A=matrix(0,m,n); Id=0*(1:(4*m)); Id10=0*(1:m)
for (i in 1:m){
for (j in 1:n){A[i,j]=x[i]^(j-1)}
Id[i]=1
Id10[i]=10
Id[m+i]=1
}
A

Id
Id10

B=matrix(0,4*m,n+1)
B[1:m,1:n]=-A
B[(m+1):(2*m),1:n]=A
B[(2*m+1):(3*m),1:n]=-A
B[(3*m+1):(4*m),1:n]=A
B[,n+1]=Id
B

yy=Id
yy[1:m]=-y        # restrição >=-y
yy[(m+1):(2*m)]=y # restrição >=y

yy[(2*m+1):(3*m)]=-Id10 # restrição <=10
yy[(3*m+1):(4*m)]=Id[1:m]  # # restrição >=1
## Solves linear and quadratic programming problems
## but needs a feasible starting value
#
# from example(solve.QP) in 'quadprog'
# 
fQP <- function(u) {u[length(u)]}
Amat       <- B
bvec       <- yy

gQP <- function(u) {v=0*(1:length(u))
    v[length(u)]=1
    v
}
####
result=constrOptim(u0, fQP,gQP, ui = B, ci = bvec)
result
ur=result$par

poli<-function(t){
a=ur[1]
b=ur[2]
c=ur[3]
return(a+b*t+c*t^2)
}

curve(poli,0,1)
points(x,y)

print("Start X")
print(x)
print("Start Y")
print(y)
print("Approximated z=f(x)")
print(poli(x))
print("p(x)=a+bx+cx^2"); print("a, b, c"); print(ur[1:3])