import numpy as np

def rowaddmult(A, i, j, c):
    A[j] += A[i] * c
    return A

def backsub(U,v):
    Copy = U.copy()
    rows,cols = Copy.shape
    ret = np.ones((rows))
    for i in range(rows-1,-1,-1):
        x = (v[i] - U[i,i+1:rows]@ret[i+1:rows]) / U[i,i]
        ret[i] = x
    return ret

def rowred(A):
    Copy = A.copy()
    rows,cols = Copy.shape
    for a in range(rows):
        for b in range(rows):
            if (b > a and Copy[a,a] != 0.0):
                Copy = rowaddmult(Copy, a, b, (-1)*Copy[b,a]/Copy[a,a])
                #print(Copy)
    return Copy

def solve(Aug):
    Ud = rowred(Aug)
    rows, cols = Ud.shape
    return backsub(Ud[:,0:rows], Ud[:,rows])



A = np.zeros((200, 3))
A[0:100,0] = (np.random.rand(100) * 5) + 5
A[0:100,1] = (np.random.rand(100) * 5) + 5
A[100:,0] = (np.random.rand(100) * 4) + 8
A[100:,1] = (np.random.rand(100) * 4) + 4
A[:,2] = np.ones(200)

v = np.ones(200)
v[100:] = -1
#v = np.zeros(200)

M = np.ones((3,4))
M[:, 0:3] = A.T@A
M[:, 3] = A.T@v
x = solve(M)

def determinefalses(A, x):
    #"Is this point in the cluster on the top left?'
    #The probability of getting a false positive (point is on top left when correct is point is on bottom right.);
    #aka
    #The probability of getting a false negative (that is, a 'no' answer when the correct answer is 'yes'.).
    false_pos = []
    false_neg = []

    for i in range(0, 100):
        fxy = (x[0] * A[i,0]) + (x[1] * A[i,1]) + x[2]
        if fxy < 0:#double check this!
            false_pos.append((A[i,0], A[i,1]))

    for i in range(100, 200):
        fxy = (x[0] * A[i,0]) + (x[1] * A[i,1]) + x[2]
        if fxy > 0:#double check this!
            false_neg.append((A[i,0], A[i,1]))

    return [len(false_pos) / 100, len(false_neg) / 100]

falses = determinefalses(A, x)

print(f'Probability of False Positive: {falses[0]}\nProbability of False Negative: {falses[1]}')



%matplotlib inline
import matplotlib.pyplot as plt
x1 = np.linspace(5,12,1000) # 1000 linearly spaced numbers - read up on the np.linspace() command!
y1 = -xhat[0]*x1/xhat[1] - xhat[2]/xhat[1]
plt.plot(x1, y1, 'b')

plt.scatter(A[0:100,0],A[0:100,1],c='b')
plt.scatter(A[100:,0],A[100:,1],c='r')