import numpy as np

class IsingLattice:

    E = 0.0
    E2 = 0.0
    M = 0.0
    M2 = 0.0

    n_cycles = 0

    def __init__(self, n_rows, n_cols):
        self.n_rows = n_rows
        self.n_cols = n_cols
        self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

    def energy(self):
        "Return the total energy of the current lattice configuration."
        energy = 0.0
        for r in range(0,self.n_rows):
            for c in range(0,self.n_cols):
                energy = energy - (self.lattice[r][c]*self.lattice[r][c-1]) - (self.lattice[r][c]*self.lattice[r-1][c])
        return energy

    def magnetisation(self):
        "Return the total magnetisation of the current lattice configuration."
        magnetisation = np.sum(self.lattice)
        return magnetisation

    def montecarlostep(self, T):
        # complete this function so that it performs a single Monte Carlo step
        energy = self.energy()
        magnetisation = self.magnetisation()
        #the following two lines will select the coordinates of the random spin for you
        random_i = np.random.choice(range(0, self.n_rows))
        random_j = np.random.choice(range(0, self.n_cols))
        self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
        deltaE=self.energy() - energy
        
        #the following line will choose a random number in the rang e[0,1) for you
        random_number = np.random.random()
        if deltaE>0:
            if random_number > np.exp(-deltaE/T):
                 self.lattice[random_i][random_j]=-self.lattice[random_i][random_j]
            
        self.E = self.E + self.energy()
        self.E2 = self.E**2 + self.energy()**2
        self.M = self.M + self.magnetisation()
        self.M2 = self.M**2 + self.magnetisation()**2
        self.n_cycles = self.n_cycles + 1
        return (self.energy(),self.magnetisation())

    def statistics(self):
        # complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
        if self.n_cycles > 0:
            averageE=self.E/self.n_cycles
            averageE2=self.E2/self.n_cycles
            averageM=self.M/self.n_cycles
            averageM2=self.M2/self.n_cycles
        else:
            averageE=self.E
            averageE2=self.E2
            averageM=self.M
            averageM2=self.M2
        return (averageE,averageE2,averageM,averageM2)
