#Nathan Fox
#Homework 26
#I give permission for this work to be posted online

#Read procedures from class/last homework
read(`C26.txt`):
read(`C27.txt`):

#Import packages
with(`Optimization`):

#Help procedure
Help:=proc():
    print(` FindRectangle(F, x, y, precision, pos:=true) `):
    print(` EstimateArea(F, x, y, N, pos:=true) `):
    print(` EstimateAreaMarkov(F, x, y, N, delta, pos:=true) `):
end:

##AUXILIARY PROCEDURES##

#FindRectangle(F, x, y, precision, pos:=true): Find a rectangle that
#F(x,y)^2<1 is contained in pretty tightly.  If pos is true, then it
#must also have x>=0 and y>=0 (and these will be the min x and the
#min y, even if they are not tightest)
#The fit tightness is specified by precision
#returns minx, maxx, miny, maxy
FindRectangle:=proc(F, x, y, precision, pos:=true)
local minx, maxx, miny, maxy, findthing:
    findthing:=proc(ismax, vr)
    local thing, last, cur, sol:
        #First, we need to find a point in the region.
        if pos then
            cur:=Maximize(0, {F^2<=1, x>=0, y>=0}):
        else
            cur:=Maximize(0, {F^2<=1}):
        fi:
        if vr = x then
            thing:=op(cur[2][1])[2]:
        else
            thing:=op(cur[2][2])[2]:
        fi:
        
        #This worked, but it was slow for regions subtending
        #small angles from the origin, so, we eliminated it
        #
        ##First, we need to find a point in the region.
        ##To do this, we will randomly intersect it with lines
        ##through the origin until one works
        #ra:=rand(0..1):
        #while true do
        #    slp:=evalf(tan(raR()*Pi/2)):
        #    if not pos then
        #        slp:=slp*(-1)^ra():
        #        sol:=fsolve({F^2=1, y=slp*x}, {x,y}):
        #    else
        #        sol:=fsolve({F^2=1, y=slp*x}, {x=0..infinity, y=0..infinity}):
        #    fi:
        #    if type(sol, set) then
        #        #Sometimes, fsolve is dumb
        #        if abs(subs(sol, F^2)-1) > precision then
        #            next:
        #        fi:
        #        if vr = x then
        #            thing:=subs(sol, x):
        #        else
        #            thing:=subs(sol, y):
        #        fi:
        #        break:
        #    fi:
        #
        #od:
        
        sol:=solve(subs({vr=thing}, F)^2<=1):
        while sol <> NULL and (not(pos) or evalf(op(sol)[2]) >=0) do
            last:=thing:
            thing:=2*abs(thing)*ismax+ismax:
            sol:=solve(subs({vr=thing}, F)^2<=1):
            #Deal with stupid corner cases
            if type(sol, numeric) then
                sol:=[sol, sol]:
            fi:
        od:
        while abs(thing-last) > precision do
            cur:=(thing+last)/2:
            sol:=solve(subs({vr=cur}, F)^2<=1):
            #Deal with stupid corner cases
            if type(sol, numeric) then
                sol:=[sol, sol]:
            fi:
            if sol = NULL or (pos and evalf(op(sol)[2]) < 0) then
                thing:=cur:
            else
                last:=cur:
            fi:
        od:
        return thing
    end:
    
    if pos then
        minx:=0:
    else
        minx:=findthing(-1, x):
    fi:
    maxx:=findthing(1, x):
    if pos then
        miny:=0:
    else
        miny:=findthing(-1, y):
    fi:
    maxy:=findthing(1, y):
    
    return minx, maxx, miny, maxy:
end:


##PROBLEM 1##

#EstimateArea(F, x, y, N): inputs a polynomial F in x,y with
#positive coefficients, and a large positive integer N, and
#estimates the area of the region
#{(x,y); F(x,y)^2<1, x>=0, y>=0}
#by picking N random points in a rectangle that contains the above region, and counts those that are inside it.
#
#NOTE: If you want to drop the x>=0, y>=0 condition, add a fifth
#argument false (pos)
EstimateArea:=proc(F, x, y, N, pos:=true)
local minx, maxx, miny, maxy, precision, i, ct, rx, ry:
    precision:=10^(-3):
    minx, maxx, miny, maxy:=FindRectangle(F, x, y, precision, pos):
    
    ct:=0:
    for i from 1 to N do
        rx:=raR()*(maxx-minx)+minx:
        ry:=raR()*(maxy-miny)+miny:
        if subs({x=rx, y=ry}, F)^2 < 1 then
            ct:=ct+1:
        fi:
    od:
    return (ct/N)*(maxx-minx)*(maxy-miny):
end:

#EstimateArea(x^2+y^2,x,y,10000); returned 0.787912117516531
#EstimateArea(x^3+2*x*y+y^3,x,y,10000); returned 0.578969213630107

##PROBLEM 2##

#EstimateAreaMarkov(F, x, y, N, delta): inputs a polynomial F in x,
#with positive coefficients, and a large positive integer N, and a
#delta less than 1, and estimates the area of
#{(x,y); F(x,y)^2<1, x>=0, y>=0}
#by picking N random points in a rectangle that contains the
#region F(x,y).
#
#NOTE: If you want to drop the x>=0, y>=0 condition, add a fifth
#argument false (pos)
EstimateAreaMarkov:=proc(F, x, y, N, delta, pos:=true)
local minx, maxx, miny, maxy, precision, i, ct, cx, cy, cxn, cyn:
    precision:=10^(-3):
    minx, maxx, miny, maxy:=FindRectangle(F, x, y, precision, pos):
    
    ct:=0:
    #Start at a random point
    cx:=raR()*(maxx-minx)+minx:
    cy:=raR()*(maxy-miny)+miny:
    for i from 1 to N do
        cxn:=cx+raR()*2*delta-delta:
        cyn:=cy+raR()*2*delta-delta:
        if cxn < minx or cxn > maxx or cyn < miny or cyn > maxy then
            #do nothing
        else
            cx:=cxn:
            cy:=cyn:
        fi:
        if subs({x=cx, y=cy}, F)^2 < 1 then
            ct:=ct+1:
        fi:
    od:
    return (ct/N)*(maxx-minx)*(maxy-miny):    
end:

#EstimateAreaMarkov(x^2+y^2, x, y, 10000, .2); returned 0.810335229487526
#EstimateAreaMarkov(x^3+2*x*y+y^3, x, y, 10000, .3); returned 0.560166965890867