import numpy as np
import math

# Matrix functions

def matrix(xs, ker):
    '''Given vectors xs = (x_1, ..., x_n) and a symmetric kernel function ker
    that satisfies ker(x, x) = 1, return the matrix (ker(x_j, x_k)).'''
    n = len(xs)
    g = np.zeros((n, n))
    for i in range(n):
        g[i,i] = 1
        for j in range(i):
            g[i,j] = ker(xs[i], xs[j])
            g[j,i] = g[i,j]

    return g

def is_semipos(g):
    '''Return True if the symmetric matrix g is semipositive.'''
    return all((v >= 0 for v in np.linalg.eigvalsh(g)))

# Kernel functions

def symmetric(xs, ys):
    '''"Kernel" that yields a random symmetric matrix with 1s on the diagonal
    and whose entries satisfy 0 \leq b_{jk} \leq 1.'''
    return np.random.rand()

# If 0 \leq xs \leq 1 and ys = (1, ..., 1) then K(xs, ys) = avg(xs).
# If xs \leq ys, then K(xs, ys) = avg(xs) / avg(ys)
def K(xs, ys):
    '''(xs, ys) -> (sum min(x, y)) / (sum max(x, y))'''
    a = 0
    b = 0
    for x, y in zip(xs, ys):
        a += min(x, y)
        b += max(x, y)
    return a / b

# Linear-algebraic upper bound on the kernel K.
def Kbound(xs, ys):
    '''(x, y) -> 1 - |x-y|/|x+y|'''
    a = 0
    b = 0
    for x, y in zip(xs, ys):
        a += (x - y)**2
        b += (x + y)**2

    return (math.sqrt(b) - math.sqrt(a)) / math.sqrt(b)

# Upper bound on Kbound
def InnerProd(xs, ys):
    '''(x, y) -> 2 \sqrt{<x, y>} / |x + y|'''
    a = 0
    b = 0
    for x, y in zip(xs, ys):
        a += x * y
        b += (x + y)**2

    return 2 * math.sqrt(a) / math.sqrt(b)

if __name__ == '__main__':
    n = 4
    weight = 100

    #ker = symmetric
    ker = K
    #ker = Kbound
    #ker = InnerProd

    pos = 0
    neg = 0
    i = 1
    while True:
        xs = weight * np.random.rand(n,n)
        g = matrix(xs, ker)

        if not is_semipos(g):
            neg += 1
        else:
            pos += 1

        if i % 500 == 0:
            print("""
====
Tried {} matrices
- {} positive,
- {} negative
- P(positive) = {:.06f}
""".strip().format(i, pos, neg, pos / i))

        i += 1