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

#Read packages
with(`combinat`):

#Read procedures from class
read(`C3.txt`):
read(`C3a.txt`):

#Help procedure
Help:=proc():
    print(` MS(N, k, M, eps) , PostV(L, M) `):
end:

##PROBLEM 2##

#MS(N, k, M, eps): inputs a positive integer N, a positive integer
#k (smaller than N, in practice much smaller), a large integer M,
#and a small number eps, and outputs the number of times the
#estimated probability was not in the interval [1/2-eps,1/2+eps]
MS:=proc(N, k, M, eps) local count, i, cu:
    count:=0:
    for i from 1 to M do
        cu:=S(N, k):
        if abs(cu-1/2) > eps then
            count:=count + 1:
        fi:
    od:
    return count:
end:

#MS(10000, 100, 1000, 3/100) outputted 467
#MS(10000, 100, 1000, 1/100) outputted 770
#MS(10000, 500, 1000, 1/100) outputted 635

##PROBLEM 3##

#PostV(L, M): checks that L and M are good inputs, and returns
#FAIL if it is not, and prints out an explanation why it is not.
#Rememeber that L is a list, M is a list of lists each of its
#elements has the same length, L can only have probabilities,
#they should add up to 1, etc. etc.

PostV:=proc(L, M) local i, j, err, lg, checkProbList:
    #Error message encapsulator
    err:=proc(msg)
        print(cat(`Error: `, msg)):
        return FAIL:
    end:
    #Probability list checker encapsulator
    checkProbList:=proc(lst, nme, tolerance) local tot, i:
        tot:=0:
        for i from 1 to nops(lst) do
            if not(type(lst[i], numeric)) then
                return err(cat(nme, `[`, i, `] is not a number`)):
            elif lst[i] < 0 then
                return err(cat(nme, `[`, i, `] is less than 0`)):
            elif lst[i] > 1 then
                return err(cat(nme, `[`, i, `] is greater than 1`)):
            fi:
            tot:=tot + lst[i]:
        od:
        #Allow floating-point probabilities
        if abs(tot - 1) > tolerance then
            return err(cat(`Entries in `, nme, ` do not add to 1`)):
        fi:
        return true:
    end:
    #Check validity of L
    if not(type(L, list)) then
        return err(`L is not a list`):
    fi:
    if checkProbList(L, "L", 1e-7) = FAIL then
        return FAIL:
    fi:
    #Check validity of M
    if not(type(M, list)) then
        return err(`M is not a list`):
    fi:
    if nops(M) <> nops(L) then
        return err(`L and M are not the same size`):
    fi:
    lg:=nops(M[1]):
    for i from 1 to nops(M) do
        if nops(M[i]) <> lg then
            return err(`Not all entries of M are the same length`):
        fi:
        if checkProbList(M[i], cat("M[", i, "]"), 1e-7) = FAIL then
            return FAIL:
        fi:
    od:
    #Check that no division by zero will occur
    for j from 1 to nops(M[1]) do
        if add(L[i]*M[i][j], i=1..nops(L)) = 0 then
            return err(cat(`Index `, j, ` results in division by zero`)):
        fi:
    od:
    return [seq([seq(L[i]*M[i][j]/add(L[i]*M[i][j], i=1..nops(L)), i=1..nops(L))], j=1..nops(M[1]))]:
end: