#OK to post homework
#Lucy Martinez, 02-10-2024, Assignment 7

with(numtheory):
with(NumberTheory):

# Problem 2:  p. 59, lines 7 and 8 from the bottom
# (that start with (a) and (b)) has two typos. Find them!
#
# Solution:
# The typos are in the channel noise box. The book has the following table:
# [Message, Codeword, Channel Noise, Received Vec, Decoded Word, Received Msg]
# [  01,      0101,     0101,          0001,         0101,         01]
# [  01,      0101,     0101,          0100,         0000,         00]
#
# But it should be:
#    [Message, Codeword, Channel Noise, Received Vec, Decoded Word, Received Msg]
# (a) [  01,      0101,     0100,          0001,         0101,         01]
# (b) [  01,      0101,     0001,          0100,         0000,         00]
# 
# To verify this, we need to check that
# codeword + channel noise = received vector (mod 2)
# Hence, using the updated channel noise:
# 0101 + 0100 = 0001, which agrees with the received vector
# 0101 + 0001 = 0100, which agrees with the received vector



# Problem 3, Part 1: Using Maple's rand(1..n) command write a procedure LC(N) that
# inputs a positive integer N, larger than 1, and outputs 1 with probability 1/N
# and 0 otherwise. For example, LC(10) should output 1 about one tenth of the times.
# [Hint: use rand(1..N)() and output 1 it is 1, and 0 otherwise]
LC:=proc(N) local ra:
ra:=rand(1..N)():
if ra=1 then
 return 1:
else
 return 0:
fi:
end:

# Problem 3, Part 2: Test LC(N) out by doing add(x[LC(10)],i=1..10000);
# and make sure that it is close to 1000*x[1]+9000*x[0]

# Answer: Running this outputs (as desired)
# 985*x[1] + 9015*x[0]
# 960*x[1] + 9040*x[0]
# 1054*x[1] + 8946*x[0]
# 1046*x[1] + 8954*x[0]





# Problem 4: Write a procedure TestCode(q,n,C,N,K) that inputs a code q-ari code,C,
# on n letters (a subset of GF(q)^n), a positive integer N and
# a large positive integer K and K times does the following:
# (1) Constructs, once and for all the table T:=DecodeT(q,n,C)
# (2) Picks a random member v, of C
#  for each component of v, with probabiliy 1/N, changes it (you can add 1 mod q),
#  getting a (possibly) distorted vector v1. 
#   See whether T[v1]<>v. If it is add it to the counter
# (3) output evalf(co/K)
TestCode:=proc(q,n,C,N,K) local co,i,T,ra,j,pr,v,v1:
co:=0:
for i from 1 to K do
 ra:=rand(1..nops(C))():
 T:=DecodeT(q,n,C):
 v:=C[ra]:
 v1:=v:
 for j from 1 to nops(v) do
  pr:=LC(N):
  if pr=1 then
   v1:=subsop(j=v1[j]+ 1 mod q,v1):
  fi:
 od:
 if T[v1]<>v then
   co:=co+1:
 fi:
od:
evalf(co/K):
end:




# Problem 5: With the above program, run it with q=2, n=7, K=10000, C=GRC(2,7,3)
# and N=2,3,4,5,6,7,8,9,10,15,20,30
# What were the failure rates for each choice of N above?

# Answer:

# N=2, TestCode(2,7,GRC(2,7,3),2,10000); outputs 0.9375000000
# N=3, TestCode(2,7,GRC(2,7,3),3,10000); outputs 0.7390000000
# N=4, TestCode(2,7,GRC(2,7,3),4,10000); outputs 0.5549000000
# N=5, TestCode(2,7,GRC(2,7,3),5,10000); outputs 0.4216000000
# N=6, TestCode(2,7,GRC(2,7,3),6,10000); outputs 0.3347000000
# N=7, TestCode(2,7,GRC(2,7,3),7,10000); outputs 0.2694000000
# N=8, TestCode(2,7,GRC(2,7,3),8,10000); outputs 0.2154000000
# N=9, TestCode(2,7,GRC(2,7,3),9,10000); outputs 0.1764000000
# N=10, TestCode(2,7,GRC(2,7,3),10,10000); outputs 0.1526000000
# N=15, TestCode(2,7,GRC(2,7,3),15,10000); outputs 0.07140000000
# N=20, TestCode(2,7,GRC(2,7,3),20,10000); outputs 0.04500000000
# N=30, TestCode(2,7,GRC(2,7,3),30,10000); outputs 0.02080000000



# Problem 6: Finish procedure SA(q,n,M) that we just started.
# Test it with q=2,n=4, M=[[1,0,1,1],[0,1,0,1]]
# and see whether you get the same as the standard array for the small code
# in Example 6.6 on p. 59 of Hill's book (make sure that your top row agrees with his).
SA:=proc(q,n,M) local SL,C,C1,A,i,a1,j:
C:=LtoC(q,M):
C1:=C:
C:=C minus {[0$n]}:
SL:=[[[0$n],op(C)]]:

# removes first row of SL
A:=Fqn(q,n) minus {op(SL[1])}:

# continue in this way until all the cosets are listed and
#  every vector in Fqn(q,n) appears exactly once
for i from 1 while A<>{} do
 a1:=minW(A):
 SL:=[op(SL),[seq(a1+C1[j] mod q,j=1..nops(C1))]]:
 A:=A minus {op(SL[i+1])}:
od:
SL:
end:

# SA(2,4,[[1,0,1,1],[0,1,0,1]]) returns
# [[[0, 0, 0, 0], [0, 1, 0, 1], [1, 0, 1, 1], [1, 1, 1, 0]], [[1, 0, 0, 0], [1, 1, 0, 1], [0, 0, 1, 1], [0, 1, 1, 0]], [[0, 1, 0, 0], [0, 0, 0, 1], [1, 1, 1, 1], [1, 0, 1, 0]], [[0, 0, 1, 0], [0, 1, 1, 1], [1, 0, 0, 1], [1, 1, 0, 0]]]
# The second and third column are switched but this array agrees with the one
#  in the book (pg 59 in book)



#minW(A): inputs a list of vectors A and finds a vector of smallest weight
minW:=proc(A) local n,cha,i:
n:=nops(A[1]):
cha:=A[1]:
for i from 2 to nops(A) do
 if HD(A[i],[0$n])<=HD(cha,[0$n]) then
  cha:=A[i]:
 fi:
od:
cha:
end:





####################################### From Previous Classes Procedures:
# C7.txt, Feb. 8, 2024
Help7:=proc():
 print(`NN(C,v),DecodeT(q,n,C),GLC1(q,M,d),GLC(q,n,d),SA(q,n,M)`):
end:


# NN=Nearest Neighbor
# NN(C,v): inputs a code C (subset of Fq(q,n) where n:=nops(C))
#  and a vector v, and finds the set of members of C closest to v
NN:=proc(C,v) local i,rec,cha:
# extracts first element of C, set of champions
cha:={C[1]}:

# the recurrent hamming distance which will be updated
rec:=HD(v,C[1]):

for i from 2 to nops(C) do
 if HD(v,C[i])<rec then
  cha:={C[i]}:
  rec:=HD(v,C[i]):
 elif HD(v,C[i])=rec then
  cha:=cha union {C[i]}:
 fi:
od:
cha:
end:



# DecodeT(q,n,C): Constructs once and for all a table that sends every member
#  of Fqn(q,n) to one of its nearest neighbors
DecodeT:=proc(q,n,C) local T,S,v:
S:=Fqn(q,n):

for v in S do
 T[v]:=NN(C,v)[1]:
od:
op(T):
end:


# GLC1(q,M,d): tries to add a new member to the current basis M
# out of the other vectors that still has minimum weight d
GLC1:=proc(q,M,d) local n,A,c,M1:
n:=nops(M[1]):
A:=Fqn(q,n) minus LtoC(q,M):
for c in A do
 M1:=[op(M),c]:
 if MinW(q,M1)=d then
  return(M1):
 fi:
od:
FAIL:
end:

# Generate good linear codes
# GLC(q,n,d): inputs q,n,d (d for min distance), greedily and randomly tries to
# to get as large a linear code given by a basis as possible
GLC:=proc(q,n,d) local M,M1:
M:=[[1$d,0$(n-d)]]:
M1:=GLC1(q,M,d):
while M1<>FAIL do
 M:=M1:
 M1:=GLC1(q,M,d):
od:
M:
end:


# SA(q,n,M): inputs a basis M of a linear [n,nops(M),d]-code, and ouputs
# Slepian's Standard Array as a matrix of vectors containing all the vectors
# in GF(q)^n (aka Fqn(q,n) ) such that the first row is an ordering of the 
# members of the actual code ( LtoC(q,M) ) with [0$n] as the first entry and the 
# the first columns are the coset representatives 
#SA:=proc(q,n,M) local SL,C,A:
#C:=LtoC(q,M):
#C:=C minus {[0$n]}:
#SL:=[[0$n],op(C)]:

# removes first row of SL
#A:=Fqn(q,n) minus {op(SL[1])}:
# look at page 58 in the book (page 35 in pdf)

# write a function that finds the minimum weight, call it minW(A) where A
#  is the list of vectors and finds a vector of smallest weight
#  (min number of zero entries) 
#   Then, we will need to use minW(A) to choose the next coset representative a1
#   The next row (of SL) is a1+SL[1] 
#   keep updating A until you run out of vectors in Fqn(q,n)
#end:

##################################### From previous classes
#Feb. 5, 2024, C6.txt#C5.txt, Feb. 1, 2024
Help6:=proc(): print(`LtoC(q,M), MinW(q,M)`):end:


Help5:=proc(): print(`Nei(q,c), SP(q,c,t), GRC(q,n,d), GRCgs(q,n,d) , MinD(C), CV(S,n)`):
print(`BDtoC(BD,n)`):
end:

#Old code
#Jan. 29, 2024 C4.txt
Help4:=proc(): print(`Fqn(q,n), HD(u,v), RV(q,n) , RC(q,n,d,K), SPB(q,n,t), BDfano(), BDex212() `):end:


#Alphabet {0,1,...q-1}, Fqn(q,n): {0,1,...,q-1}^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:

#Def. (n,M,d) q-ary code is a subset of Fqn(q,n) with M elements with
#minimal Hamming Distance d between any two members
#It can detect up to d-1 errors
#
#If d=2*t+1 correct t errors.
#
# Hamming distance between u and v: counts the number of places in which
#  u and v differs
#HD(u,v): The Hamming distance between two words (of the same length)
HD:=proc(u,v) local i,co:
co:=0:
for i from 1 to nops(u) do
  if u[i]<>v[i] then
      co:=co+1:
  fi:
od:
co:
end:

#SPB(q,n,d): 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:


#RV(q,n): A random word of length n in {0,1,..,q-1}
RV:=proc(q,n) local ra,i:
ra:=rand(0..q-1):
[seq( ra(), i=1..n)]:
end:

#RC(q,n,d,K): inputs q,n,d, and K and keeps picking K+1 (thanks to Nuray) random vectors
#whenver the new vector is not distance <=d-1 from all the previous one 
RC:=proc(q,n,d,K) local C,v,i,c:
C:={RV(q,n)}:
for i from 1 to K do
 v:=RV(q,n):
  if min(seq(HD(v,c), c in C))>=d then
   C:=C union {v}:
  fi:
od:
C:
end:

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:

#end of old code

#Nei(q,c): 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:

#SP(q,c,t): the set of all vectors in Fqn(q,n) whose distance is <=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:

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:


#GRCgs(q,n,d): George Spahn's version
GRCgs:=proc(q,n,d) local S,A,v:
print(`Warning: use at your own risk`):
A:=Fqn(q,n):
S:={}:
while A<>{} do:
 v:=A[rand(1..nops(A))()]:
 S:=S union {v}:
 A:=A minus SP(q,v,d-1):
od:
S:
end:

#MinD(C): The minimal (Hamming) distance of the code C
MinD:=proc(C) local i,j:
min( seq(seq(HD(C[i],C[j]),j=i+1..nops(C)), i=1..nops(C))):
 end:

#CV(S,n): the characteristic vector of the subset S of {1,...,n}
CV:=proc(S,n) local v,i:
v:=[]:
for i from 1 to n do
 if member(i,S) then
  v:=[op(v),1]:
 else
  v:=[op(v),0]:
 fi:
od:
v:
end:

BDtoC:=proc(BD,n) local s, C:
C:={seq(CV(s,n),s in BD)}:
C:=C union subs({0=1,1=0},C):
C union {[0$n],[1$n]}:
end:

##end of old stuff
#LtoC(q,M): inputs a list of basis vectors for our linear code over GF(q)
#outputs all the codewords (the actual subset of GF(q)^n with q^k elements
LtoC:=proc(q,M) local n,k,C,c,i,M1:
option remember:
k:=nops(M):
n:=nops(M[1]):
if k=1 then
 RETURN({seq(i*M[1] mod q,i=0..q-1) }):
fi:
M1:=M[1..k-1]:
C:=LtoC(q,M1):
{seq(seq(c+i*M[k] mod q,i=0..q-1),c in C)}:
end:

#MinW(q,M): The minimal weight of the Linear code generated by M over GF(q)
MinW:=proc(q,M) local n,C,c:
n:=nops(M[1]):
C:=LtoC(q,M):

min( seq(HD(c,[0$n]), c in C minus {[0$n]} )):

end: