3834.sagews

ttr2 = (0, 0); ttr3 = (1, 0); ttr1 = (0.5, 0.866025403784439); ttr = (ttr1, ttr2, ttr3)

defaultcubatur = [(0.736712498969535, 0.237932366475534, 0.025355134554931, 0.104166666666667), (0.333333333333333, 0.333333333333333, 0.333333333333333, 0.0625)]

def S0 (P, (tr1, tr2, tr3), f=lambda x,y:x , cubatur=defaultcubatur):
    answer = 0;
    for (alpha1, alpha2, alpha3, w) in cubatur:
       answer = answer + (
    f(P(alpha1*tr1[0]+alpha2*tr2[0]+alpha3*tr3[0], alpha1*tr1[1]+alpha2*tr2[1]+alpha3*tr3[1]), (alpha1, alpha2, alpha3)) +
    f(P(alpha1*tr1[0]+alpha3*tr2[0]+alpha2*tr3[0], alpha1*tr1[1]+alpha3*tr2[1]+alpha2*tr3[1]), (alpha1, alpha3, alpha2)) + 
    f(P(alpha2*tr1[0]+alpha1*tr2[0]+alpha3*tr3[0], alpha2*tr1[1]+alpha1*tr2[1]+alpha3*tr3[1]), (alpha2, alpha1, alpha3)) +
    f(P(alpha3*tr1[0]+alpha1*tr2[0]+alpha2*tr3[0], alpha3*tr1[1]+alpha1*tr2[1]+alpha2*tr3[1]), (alpha3, alpha1, alpha2)) +
    f(P(alpha2*tr1[0]+alpha3*tr2[0]+alpha1*tr3[0], alpha2*tr1[1]+alpha3*tr2[1]+alpha1*tr3[1]), (alpha2, alpha3, alpha1)) +
    f(P(alpha3*tr1[0]+alpha2*tr2[0]+alpha1*tr3[0], alpha3*tr1[1]+alpha2*tr2[1]+alpha1*tr3[1]), (alpha3, alpha2, alpha1)))*w
    return answer

def isvalidpoint ((i,j), size):
    return (i>=0)&(i<size)&(j>=0)&(j<=i)
    
def weight ((i0,j0), (i1,j1), size): #return 1 if true and 2 if false
    if (i0,j0) == (i1,j1):
        return len(neartr((i0,j0),size))*2
    if (j0==0)&(j1==0)|((i0-j0)==0)&((i1-j1)==0)|(i0==size-1)&(i1==size-1):
        return 1
    return 2

        
def isvalidtr (tr, size):
    return reduce (lambda b, point: b & isvalidpoint (point, size), tr, True)

def neartr ((i,j), size):
    return filter (lambda tr: isvalidtr(tr, size), [((i,j), (i,j-1), (i-1,j-1)), ((i,j), (i-1,j-1), (i-1,j)), ((i,j), (i-1,j), (i,j+1)), ((i,j), (i,j+1), (i+1,j+1)), ((i,j), (i+1, j+1), (i+1,j)), ((i,j), (i+1,j), (i,j-1))])

def weights ((i,j),size):
    return map(lambda (i0,j0): (weight((i,j),(i0,j0), size)*1/12, pointtonom((i0,j0))),
        filter (lambda point: isvalidpoint(point, size), [(i,j), (i,j-1), (i-1,j-1), (i-1,j), (i,j+1), (i+1, j+1), (i+1,j)]))


def pointtonom ((i,j)):
    return i*(i+1)/2+j

def pointtocord ((i,j), size):
    return (ttr1[0]+i*(ttr2[0]-ttr1[0])/(size-1)+ j*(ttr3[0]-ttr2[0])/(size-1),
            ttr1[1]+i*(ttr2[1]-ttr1[1])/(size-1)+ j*(ttr3[1]-ttr2[1])/(size-1))
            
def trtocord (tr, size):
    return map (lambda point:pointtocord(point,size), tr)

def multbyg (f, point, size):
    return reduce(lambda x,y: x+y, 
                  map(lambda tr: S0(f, tr, lambda x,(alpha1,alpha2,alpha3):x*alpha1), 
                      map(lambda tr:trtocord(tr,size),
                          neartr(point,size))))

def rightside (f, size):
    result=[]
    for i0 in [0..(size-1)]:
        for j0 in [0..i0]:
            result.append((multbyg (f, (i0,j0), size), i0, j0))
    return result

def zero(size, a=0):
    result = []
    for i0 in [0..(size-1)]:
        for j0 in [0..i0]:
            result.append((a, i0, j0))
    return result

def multbyA (x, size):
#    print 'start'
    result = []
    for i0 in [0..(size-1)]:
        for j0 in [0..i0]:
            result.append((reduce (lambda tmp, (a,n): tmp+a*x[n][0], weights((i0,j0),size), 0), i0, j0))
#            result.append((weights((i0,j0),size), i0, j0))
#    print 'end'
    return result

def l2prod (l1, l2):
    return reduce (lambda x, (y1,y2): x+y1[0]*y2[0], zip(l1,l2), 0)
    
def l2norm (l1):
    return numerical_approx(sqrt(reduce(lambda x, y: x+y[0]*y[0], l1, 0)))

def zipwith (l1,l2,f):
    return map(f, zip(l1,l2))

def mult (l1, alpha):
    return map(lambda (x,i,j): (alpha*x,i,j),l1)
    
def minus (l1, l2):
    return zipwith (l1,l2, lambda ((x1,i1,j1),(x2,i2,j2)):((x1-x2),i1,j1))

def solveLanc (size,F,eps):
    c0 = zero (size)
    r0 = minus(F,multbyA(c0,size))
    xi1 = l2norm(r0)
    beta = xi1
    if (beta < eps):
        return c0
    v1 = mult(r0,1/beta)
    lambda1 = 0
    beta1 = 0 
    w = multbyA(v1,size)
    alpha1 = l2prod(w,v1)
    eta1 = alpha1-lambda1*beta1
    p1 = mult(v1,1/eta1)
    c1 = minus(c0, mult(p1, -xi1))
    if (l2norm(minus(multbyA(c1, size), F))/beta < eps):
        return c1
    w = minus (w, mult(v1,alpha1))
    while True:
        print 'c1'
        print c1
        beta2 = l2norm(w)
        v2 = mult(w,1/beta2)
        w = minus(multbyA (v2,size), mult(v1, beta2))
        alpha2 = l2prod(w,v2)
        lambda2 = beta2/eta1
        xi2 = -lambda2*xi1
        eta2 = alpha2 - lambda2*beta2
        p2 = mult(minus(v2, mult(p1,beta2)), 1/eta2)
        c2 = minus (c1, mult(p2, -xi2))
        if (l2norm(minus(multbyA(c2, size), F))/beta < eps):
            return c2
        w = minus (w, mult(v2, alpha2))
        eta1 = eta2
        p1 = p2
        c1 = c2
        xi1=xi2

def solveCG (size, F, eps):
    c0 = zero(size)
    r0 = minus(F, multbyA(c0, size))
    p0 = r0
    beta = l2norm (r0)
    if (beta < eps):
        return c0
    while (l2norm(r0)/beta > eps):
#        print 'c0'
#        print c0
        alpha0 = l2prod(r0,r0)/l2prod(multbyA(p0, size), p0)
        c1 = minus (c0, mult(p0, -alpha0))
        r1 = minus (r0, mult(multbyA(p0, size),alpha0))
        beta0 = l2prod(r1,r1)/l2prod(r0,r0)
        p1 = minus(r1, mult(p0, -beta0))
        c0 = c1
        p0 = p1
        r0 = r1
    return c0

def errortr (f, tr, resh, size):
    return S0(f, (pointtocord(tr[0], size), pointtocord(tr[1], size), pointtocord(tr[2], size)), lambda x,(alpha1,alpha2,alpha3):(x - alpha1*resh[pointtonom(tr[0])][0] - alpha2*resh[pointtonom(tr[1])][0] - alpha3*resh[pointtonom(tr[2])][0])^2)

def split (size):
    split = []
    for i in [1..size-1]:
        for j in [0..i-1]:
            split.append (((i,j),(i-1,j),(i,j+1)))
    for i in [1..size-2]:
        for j in [0..i-1]:
            split.append (((i,j),(i,j+1),(i+1,j)))
    return split

def error (f, answ, size):
    return reduce (lambda a,b:a+b,
                   map(lambda tr: errortr(f,tr, answ, size), split(size)))

def test (f, size, eps):
    return error (f, solveCG(size, rightside (f, size), eps), size)

def tests (f, eps):
    return (test(f, 9, eps), test (f, 17, eps), test (f, 33, eps), test (f, 65, eps))

tests (lambda x,y:1, 1e-8)

(7.90603477079775e-25, 1.95721323540227e-14, 3.41497026894274e-14, 4.69870562654060e-13)

tests (lambda x,y:x+y, 1e-8)

(1.78399205377903e-14, 6.35929066185160e-15, 7.30622044520375e-14, 5.01785128201991e-13)

tests (lambda x,y: x*x, 1e-8)

(0.000332724780072078, 0.0000884425941907039, 0.0000227573485156689, 5.76948928539569e-6)