HW4 solution, for HW4 Due Monday September 19, 2016.
---------------------------------------------------

SOLUTION TO (a)

### The idea is to test the inequality at MANY randomly chosen
###   t,y1,y2  combinations. It is true that the function 
###   could fail to be convex in some local domains, so it is 
###   a good idea to generate the random y1 and y2 points close 
###   together sometimes.

### Most of the complexity would be to code this to apply generally 
###   for vector arguments, so to keep the solution simple, I code 
###   ONLY for ffcn taking scalar arguments. I intended a slightly
###   more ambitious function to be called "ConvexTst", which is 
###   why the error message in this function refers to that name.

ConvexTst.Scalar = function(ffcn, Dom, npt1=1000, npt2=1000, 
             diam = .01, seed=NULL, tol=1.e-10, Min=-100, Max=100) {
    ## ffcn must be a function of vector argument x with dimension 
    ##   (if specified) equal to number of rows of matrix Dom
       K = nrow(Dom)
       if(K>1) stop("Vector argument of ffcn. Use ConvexTst instead")
       if(!is.null(seed)) .Random.seed = seed
       Dmin = max(Dom[1,1],Min)
       Dmax = min(Dom[1,2],Max)
    ## Restrict domain to (-Min,Max) in case of infinite endpoints
       y1=runif(npt1, Dmin, Dmax)
       y2=runif(npt1, Dmin, Dmax)
       t = runif(npt1)
       Convex=T
       ### Assume ffcn does NOT vectorize
       for(i in 1:npt1) 
           Convex = Convex & ( ffcn(y1[i]+t[i]*(y2[i]-y1[i])) 
              <= t[i]*ffcn(y2[i])+(1-t[i])*ffcn(y1[i]) + tol)
    ### NB. In this inequality, I allow a small extra tolerance
    ### "tol" for computational error. Taking tol=0 is OK too.
    ###
    ### If 1st test passed, do 2nd test with y1, y2 closer than diam
       if(Convex) {
          y1=runif(npt2, Dmin+diam, Dmax-diam)
          y2=runif(npt2, y1-diam,y1+diam)
          t= runif(npt2)
          for(i in 1:npt2)
              Convex = Convex & ( ffcn(y1[i]+t[i]*(y2[i]-y1[i])) 
              <= t[i]*ffcn(y2[i])+(1-t[i])*ffcn(y1[i]) + tol)
       }
       Convex   }

> ConvexTst.Scalar( function(x) x^2, array(c(-2,-2,3,3),c(2,2)) )
Error in ConvexTst.Scalar(function(x) x^2, array(c(-2, -2, 3, 3), c(2,  : 
  Vector argument of ffcn. Use ConvexTst instead

> ConvexTst.Scalar( function(x) x^2, array(c(-2,3),c(1,2)) )
[1] TRUE
> ConvexTst.Scalar( function(x) x^2, array(c(-2,3),c(1,2)), 
            npt1=1e6,npt2=1e6 )
[1] TRUE                         ### took about 15secs on my laptop

> ConvexTst.Scalar( function(x) x^2-.5*x^3, array(c(-2,1),c(1,2)) )
[1] FALSE
> ConvexTst.Scalar( function(x) x^2-.3*x^3, array(c(-2,1),c(1,2)) )
[1] TRUE

### Both of these answers are correct, by 2nd deriv test
### Now do test on requested function  h

> h = function(x, probvec, cvec)  
    c(outer(x,cvec, function(x,y) abs(x-y)) %*% probvec)

> curve(h(x, c(.2,.5,.3), c(3,1,-1)), -3,7)
 ### You can see this curve is convex but not twice differentiable
 ###    since it has corners

> ConvexTst.Scalar(function(x) h(x, c(.2,.5,.3), c(3,1,-1)), 
+             array(c(-3,7),c(1,2)) )
[1] TRUE
> ConvexTst.Scalar(function(x) h(x, c(.2,.5,.3), c(3,1,-1)), 
+             array(c(-3,7),c(1,2)), npt1=1e5, npt2=1e5 )
[1] TRUE                     ### this one took around 10 seconds


### Now for speedup: 1st test not needed at all, just do 2nd test
###    using the Hessian (2nd deriv) at y1: do not need t,y2

> ConvexTst.S2 = function(ffcn.Hess, Dom, npt2=1000, 
             seed=NULL, tol=1.e-10, Min=-100, Max=100) {
    ## ffcn must be a function of vector argument x with dimension 
    ##   (if specified) equal to number of rows of matrix Dom
       K = nrow(Dom)
       if(K>1) stop("Vector argument of ffcn. Use ConvexTst instead")
       if(!is.null(seed)) .Random.seed = seed
       Dmin = max(Dom[1,1],Min)
       Dmax = min(Dom[1,2],Max)
    ## Restrict domain to (-Min,Max) in case of infinite endpoints
       Convex=T
       ### Assume ffcn.Hess does NOT vectorize
       y1=runif(npt2, Dmin, Dmax)
       for(i in 1:npt2)
           Convex = Convex & ( ffcn.Hess(y1[i]) >= -tol )
       Convex   }


> ConvexTst.S2( function(x) 2-3*x, array(c(-2,1),c(1,2)), npt2=1e6 )
[1] FALSE             ### about 4 secs



SOLUTION TO (b)

> f1 = function(x,a)   3*x+a*x^2
  f2 = function(x,b)   (x-b)^4

> h = function(x) 3*(x-5)^4 + 2* ( (x-5)^4 )^2
> xv = runif(100)

> sum(abs(f1(f2(xv,5),2) - h(xv) ))
[1] 0

> Compos1 = function(f1,f2) {
          function(x,a,b)  f1( f2(x,b), a)
  }

> Compos1(f1,f2)(xv[1],2,5)
[1] 574008.9
> f1( f2(xv[1],5), 2)
[1] 574008.9

### Complete solution might try a few additional a,b pairs and 
###  at least one more example  f1,f2

> sum(abs(f1(f2(xv,2),1) - Compos1(f1,f2)(xv,1,2) ))
[1] 0
> g1 = function(x,a) x+a; g2=function(x,b) b*x
> sum(abs(g1(g2(xv,2),1) - Compos1(g1,g2)(xv,1,2) ))
[1] 0
> sum(abs(g1(g2(xv,8),10) - Compos1(g1,g2)(xv,10,8) ))
[1] 0