# OK to post homework
# Robert Dougherty-Bliss, 2020-03-22, Assignment 15

randomize():
with(Statistics):

# 1.

# Just playing around with C15.txt.
read "C15.txt":

# ExactAv can do symbolic stuff as well. That's pretty cool.
simplify(ExactAv(2, 2, x, y, M, E, M));
simplify(ExactAv(2, 2, x, y, E * M^2, E, M));

# I'm a little confused. Where's y? It doesn't show up here. Isn't that weird?

# (There was a bug in C15.txt, but I fixed it. Essentially, I think that it
# always assumed that y = 1. I was also suspicious that the expected value of M
# was always zero. That happens when y = 1, but not always. There should now be
# y's in the output.)

# Okay, now the homework.

for k from 1 to 10 do
    NickAv(4, 4, 1.3, 1.0, 1000, 10000, E, E, M);
od;

(* Output from the above:

                                              10.57600000

                                              10.37920000

                                              11.59080000

                                              11.41560000

                                              10.14600000

                                              10.18560000

                                              10.47040000

                                              10.96240000

                                              11.21280000

                                              10.41040000

    Seems pretty close!

*)

# 2.

# (This takes a lot longer.)

for k from 1 to 10 do
    NickAv(40, 40, 1.3, 1.0, 1000, 10000, E, E, M);
od;

(* Output:

                                              854.7888000

                                              931.6104000

                                              904.2460000

                                              853.7996000

                                              857.2900000

                                              839.6788000

                                              878.7384000

                                              871.0188000

                                              929.8976000

                                              870.1108000

    Not so close!

*)

# 3.

# This is just NickAv, but with a "starting matrix" parameter.
NickAvGiven:=proc(A1, M,N,x,y,K1,K2,G,En,Ma) local A,i,su:
A := A1:

for i from 1 to K1 do
 A:=OSM(A,M,N,x,y):
od:

su:=0:

for i from 1 to K2 do
 A:=OSM(A,M,N,x,y):
 su:=su+subs({En=Ene(A),Ma=Mag(A)},G):
od:

evalf(su/K2):
end:

NickAvONE := proc(M,N,x,y,K1,K2,G,En,Ma) local A, i, j:
    A := [seq([seq(1, j=1..N)], i=1..M)]:
    NickAvGiven(A, M, N, x, y, K1, K2, G, En, Ma);
end proc:

for k from 1 to 10 do
    NickAvONE(4, 4, 1.3, 1.0, 1000, 10000, E, E, M);
od;

(* Output:

                                              10.76800000

                                              11.34440000

                                              10.89880000

                                              10.47200000

                                              9.966800000

                                              11.67960000

                                              11.10320000

                                              11.26560000

                                              11.53400000

                                              10.94040000
*)


NickAvALT := proc(M,N,x,y,K1,K2,G,En,Ma) local A, i, j:
    A := [seq([seq((-1)^(i + j), j=1..N)], i=1..M)]:
    NickAvGiven(A, M, N, x, y, K1, K2, G, En, Ma);
end proc:

for k from 1 to 10 do
    NickAvALT(4, 4, 1.3, 1.0, 1000, 10000, E, E, M);
od;

(* Output:

                                              11.55840000

                                              10.05040000

                                              11.09120000

                                              11.16080000

                                              11.26960000

                                              10.39000000

                                              11.24440000

                                              10.14600000

                                              10.69200000

                                              10.37720000

*)

# 4.

with(ArrayTools):

# Shrinks a 3^k by 3^k {-1, 1}-matrix into a 3^(k - 1) by 3^(k - 1) matrix
# using the "rule of the majority" for each block of size 3.
Shrink := proc(A, k) local i, j, x, y, L, s:
    L := Array([]):

    for i from 1 to 3^(k - 1) do
        for j from 1 to 3^(k - 1) do
            x := 3 * i - 1:
            y := 3 * j - 1:

            s := A[x - 1, y - 1] + A[x - 1, y] + A[x - 1, y + 1] +
                 A[x, y - 1] + A[x, y] + A[x, y + 1] +
                 A[x + 1, y + 1] + A[x + 1, y] + A[x + 1, y + 1]:

            if s > 0 then
                Append(L, 1):
            else
                Append(L, -1):
            fi:
        od:
    od:

    Reshape(L, [3^(k - 1), 3^(k - 1)]);
end proc:

interface(rtablesize=infinity):

A := Nick(81, 81, 0.9, 1, 10);

for i from 1 to 4 do
    A := Shrink(A, 4 - i + 1);
od;