# Please do not post homework
# Aurora Hiveley, 02/01/24, Assignment 5

Help:= proc(): print(`Fqn(q,n), Nei(q,c), SP(q,c,t), GRC(q,n,d), SPB(q,c,t), card_array(q,nLB,nUB,d_list), spb_array(q,nLB,nUB,d_list), parabd(BD,v), BDex212(), BDfano(), plotkinbound(n,d), bound_compare() `): end:


### hill textbook chapter 2 exercises

## 2.1: construct binary (n,M,d) codes where possible
# (6,2,6) = { (0,0,0,0,0,0), (1,1,1,1,1,1) }

# (3,8,1) = { (0,0,0), (0,0,1), (0,1,0), (1,0,0), (0,1,1), (1,0,1), (1,1,0), (1,1,1) }

# (4,8,2) = { (0,0,0,0), (0,0,1,1), (0,1,0,1), (1,1,0,0), (1,0,1,0), (1,0,0,1), (1,1,1,1) }

# (5,3,4) = cannot exist. distance 4 in words of length n=5
# pick your first word. WLOG, (0,0,0,0,0)
# pick your second word. it must differ from the first by at least 4 bits, so only one letter can be the same
# WLOG, pick (1,1,1,1,0)
# now notice that we cannot pick another word which differs from both (00000) and (11110) in 4 positions
# for the first letter, WLOG we pick 0. then that is THE letter which must be in common with the first word
# for the second letter, we MUST differ from word 1, so we make the second letter a 1, which agrees with word 2
# now note regardless of what we pick for letter 3, it is impossible to still differ by distance at least 4 
# from both words 1 and 2

# (8,30,3) = cannot exist
# using a helpful hint from natasha: words must be at least distance 3 apart, so each has a sphere of radius two around itself
# spheres may overlap between words so long as the *centers* of spheres do not fall in any other sphere
# thus the most any sphere could overlap with another is by 1, and this is equivalent to formalizing *non-overlapping* spheres of radius 1
# there are q^n = 2^8 = 256 total words in F(q,n) = F(2,8)
# each has 8 words inside of a sphere of radius 8 around it, so there are 9 elements per sphere
# for our code to have 30 words, we need 30*9 = 270 elements in F(2,8)
# since we only have 256, the pigeonhole principle says that this code cannot exist

## 2.2: if there exists a code (n,M,d), then there also exists a code (n-1,M',d)
## where M' <= M/2
# credit to kaylee for this idea!
# consider the binary tree producing all word of length k at level k of the tree. 
# consider level n, where M of the binary strings are selected for our code.
# note that the level above this consists of all words of length n-1, 
# and there are half as many of these words as there are of length n.
# worst case scenario, if all words of M are in pairs with a common root in level n-1,
# then we can select only one word from each pair in M to form M', meaning M' has half as many words as M.
# however, if words of M are not all paired up, then we may be able to keep more than M/2 
# since there are more than M/2 roots in level n-1 which are eligible.
# from this conclusion, we see that at best, the words of M' in level n-1 each produce two words for M in level n.
# then A2(n,d) <= 2 A2(n-1,d) since this was an upper bound for how many words we could select in level n.

## 2.3: prove that Aq(3,2)=q^2
# consider the set of all words of length 3 with distance at least two between them
# delete the last letter from each word, and note that the remaining "sub-words" consisting of the first two letters
# must all be distinct, otherwise the original words would not have differed by at least 2 letters to begin with.
# then there are q*q = q^2 distinct sub-words of length two, and the last letter is effectively forced by the distance of 2 between words.
# thus Aq(3,2) = q^2

## 2.4: let En be the subset of binary strings of length n with even weight.
## then En is the code obtained by adding a parity check to binary strings of length n-1.
## additionally, En is an (n, 2^{n-1}, 2) code
# the set of even weight strings of length n can be obtained from strings of length n-1 by
# adding a 0 in position n if the string of length n-1 had even weight
# or adding a 1 in position n if the string of length n-1 had odd weight.
# then En is a code consisting of words of length n, 
# of which there are 2^{n-1} since we build En from all strings of length n-1
# since each word of length n-1 is unique, the base words must differ from each other
# in at least one position. then, to build the even weight string of length n, 
# we add another distinct bit, so each word in En differs in at least 2 positions from all other words



### GRC(q,n,d) activities

## copied from C5

# constructs all words from alphabet {0,1,...q-1} of length n
Fqn:=proc(q,n) local S,a,v:
option remember:
if n=0 then
 RETURN({[]}):
fi:

S:=Fqn(q,n-1):

{seq(seq([op(v),a],a=0..q-1), v in S)}:

end:

# constructs set of all the neighbors of the vector c in Fqn(q,n)
Nei:=proc(q,c) local n,i,a: 
n:=nops(c):
{seq(seq([op(1..i-1,c),a,op(i+1..n,c)], a=0..q-1) , i=1..n)}:
end:

## sphere construction, set of all vectors in Fqn(q,n) w dist <= t from c
SP:=proc(q,c,t) local S,s,i:
S:={c}:
for i from 1 to t do
 S:=S union {seq(op(Nei(q,s)),s in S)}:
od:
S:
end:

## greedy random code, generates set of words in code s.t. dist btwn any two words at least d
GRC:=proc(q,n,d) local S,A,v:
A:=Fqn(q,n):
S:={}:
while A<>{} do:
 v:=A[1]:
 S:=S union {v}:
 A:=A minus SP(q,v,d-1):
od:
S:
end:

## sphere packing bound: the best you can hope for (as far as the size of C) for q-ary (n,2*t+1) code
SPB:=proc(q,n,t) local i:
trunc(q^n/add(binomial(n,i)*(q-1)^i,i=0..t)):
end:


## builds an array of code cardinalities for d in d_list and LB <= n <= UB
## generalized proc from code written for hw4

card_array:= proc(q,nLB,nUB,d_list) local UB, LB, A, i,n:
# normalize range on n for indexing purposes in array
UB:= nUB - (nLB - 1):
LB:= 1:

A:= array(1..UB, 1..nops(d_list), []): 	# initialize empty table

for i from 1 to nops(d_list) do
	for n from nLB to nUB do
		A[n - (nLB-1),i] := nops(GRC(q,n,d_list[i]));
	od:
od:

op(A);

end:

## builds an array of sphere packing bounds for d in d_list and LB <= n <= UB

spb_array:= proc(q,nLB,nUB,d_list) local UB, LB, A, i,n,t:
# normalize range on n for indexing purposes in array
UB:= nUB - (nLB - 1):
LB:= 1:

A:= array(1..UB, 1..nops(d_list), []): 	# initialize empty array

for i from 1 to nops(d_list) do
	for n from nLB to nUB do
        t := (d_list[i]-1)/2:    # d = 2t + 1
		A[n - (nLB-1),i] := SPB(q,n,t);
	od:
od:

op(A);

end:


## problem 3 calls
# q = 2; d = 3,5,7; 5 <= n <= 20
ds := [3,5,7]:
# card_array(2,5,20,ds);
# spb_array(2,5,20,ds);

## problem 4 calls
# q = 3; d = 3,5,7; 5 <= n <= 10
# card_array(3,5,10,ds);
# spb_array(3,5,10,ds);

## problem 5 calls
# q = 5; d = 3,5,7; 5 <= n <= 7
# card_array(5,5,7,ds);
# spb_array(5,5,7,ds);


## obviously the SPB >= M since SPB is the best that we can do. 
## the distance between M and SPB increases as n increases, 
## which makes sense since in the selection process there are more opportunities
## to pick a sub-optimal word for the code.
## the same is true for bigger values of q since there are more words
## of a fixed length to pick from (|F(q,n)| = q^n).



### input potential block design, output [b,v,r,k,lambda] or fail if not a block design
# BD is a set of sets, v is the number of sets (potential blocks)

with(combinat):

parabd := proc(BD, v) local b,r,k,lambda, i, elements,ele,r1, pairs,p,l1 :
# check each condition of definition one by one

# each block contains exactly k points
k := nops(BD[1]):
for i from 2 to v do
    if not k = nops(BD[i]) then
        return FAIL 
    fi:
od:


# each point lies in exactly r blocks
# initialize r
r:= 0:
elements := [BD[1][1]]:
for i from 1 to v do
    if member(BD[1][1], BD[i]) then
        r += 1:
    fi:
od:

# build elements set
elements := BD[1] :
for i from 2 to v do
    elements := elements union BD[i] :
od:

# check each element
for ele in elements do 
    r1 := 0:
    for i from 1 to v do
        if member(ele,BD[i]) then 
            r1 += 1 :
        fi:
    od:

    if not r = r1 then 
        return FAIL
    fi:
od:


# each pair of points occurs together in exactly lambda blocks
pairs := choose(elements,2):

# initialize lambda
lambda := 0:
for i from 1 to v do 
    if (member(pairs[1][1],BD[i]) and member(pairs[1][2],BD[i])) then
        lambda += 1:
    fi:
od:

# check each pair
for p in pairs do 
    l1 := 0:
    for i from 1 to v do
        if member(p[1],BD[i]) and member(p[2],BD[i]) then 
            l1 += 1:
        fi:
    od:

    if not lambda = l1 then 
        return FAIL
    fi:
od:


[BD,v,r,k,lambda];
end:


## copied from C5 for test calls
BDfano:=proc():
{{1,2,4},{2,3,5},{3,4,6},{4,5,7},{5,6,1},{6,7,2},{7,1,3}}:
end:

BDex212:=proc():
{{1,3,4,5,9},
{2,4,5,6,10},
{3,5,6,7,11},
{1,4,6,7,8},
{2,5,7,8,9},
{3,6,8,9,10},
{4,7,9,10,11},
{1,5,8,10,11},
{1,2,6,9,11},
{1,2,3,7,10},
{2,3,4,8,11}
}
end:




### plotkinbound activity
## calculates bound given in exercise 2.21 of hill
plotkinbound := proc(n,d) :
    2*floor(d/(2*d - n));
end:

## compares plotkin results to sphere packing bound
## outputs array representation for d from 2 to 20 and n from 2 to 2d of which bigger
bound_compare := proc() local n,d,s,p,C:
# initialize array
C:= array(1..19, 1..38): # rows are n's, columns are d's

    for d from 2 to 20 do
        for n from 2 to (2*d - 1) do    # must have n < 2d to avoid divide by zero
            s:= SPB(2,n,(d-1)/2):   # d = 2t + 1
            p:= plotkinbound(n,d):

            if s < p then 
                C[d-1, n-1]:= 'p':
            elif p < s then  
                C[d-1, n-1]:= 's':
            else
                C[d-1, n-1]:= 'e':
            fi:
        od:

        # make the array look nicer
        for n from 2*d to 39 do
            C[d-1, n-1]:= 0 ;
        od:

    od:
op(C);
end:

## in general, the sphere bound is tighter 
## (a few = cases, one instance of plotkin being tighter)