MLP - what did author have in mind? I am currently going through Machine learning and pattern recognition and I have following dilemma:
A algorithm for MLP is presented in that way

and corresponding python code to equation 4.9 (which is of interest to me) is written that way
deltah = self.hidden*self.beta*(1.0-self.hidden)*(np.dot(deltao,np.transpose(self.weights2)))

where self.beta is a constant (not shown in algorithm), self.hidden is hidden layer neurons' output values, deltao is error of output layer and self.weights2 are weights of output layer neurons.
My question is, why author is multiplying deltao matrix by weights2 matrix and why weights2 matrix is tranposed? Why (and how) is it equivalent to 4.9 equation?
I know that it is required get matrices to right dimensions to make multiplication possible but why does it make sense from algorithm point of view? Why can't i somehow multiply weights2 by deltao just like it is stated in equation 4.9?
Full python code is here:
# Code from Chapter 4 of Machine Learning: An Algorithmic Perspective (2nd Edition)
# by Stephen Marsland (http://stephenmonika.net)

# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.

# Stephen Marsland, 2008, 2014

import numpy as np

class mlp:
""" A Multi-Layer Perceptron"""

def __init__(self,inputs,targets,nhidden,beta=1,momentum=0.9,outtype='logistic'):
    """ Constructor """
    # Set up network size
    self.nin = np.shape(inputs)[1]
    self.nout = np.shape(targets)[1]
    self.ndata = np.shape(inputs)[0]
    self.nhidden = nhidden

    self.beta = beta
    self.momentum = momentum
    self.outtype = outtype

    # Initialise network
    self.weights1 = (np.random.rand(self.nin+1,self.nhidden)-0.5)*2/np.sqrt(self.nin)
    self.weights2 = (np.random.rand(self.nhidden+1,self.nout)-0.5)*2/np.sqrt(self.nhidden)

def earlystopping(self,inputs,targets,valid,validtargets,eta,niterations=100):

    valid = np.concatenate((valid,-np.ones((np.shape(valid)[0],1))),axis=1)

    old_val_error1 = 100002
    old_val_error2 = 100001
    new_val_error = 100000

    count = 0
    while (((old_val_error1 - new_val_error) > 0.001) or ((old_val_error2 - old_val_error1)>0.001)):
        count+=1
        print count
        self.mlptrain(inputs,targets,eta,niterations)
        old_val_error2 = old_val_error1
        old_val_error1 = new_val_error
        validout = self.mlpfwd(valid)
        new_val_error = 0.5*np.sum((validtargets-validout)**2)

    print "Stopped", new_val_error,old_val_error1, old_val_error2
    return new_val_error

def mlptrain(self,inputs,targets,eta,niterations):
    """ Train the thing """    
    # Add the inputs that match the bias node
    inputs = np.concatenate((inputs,-np.ones((self.ndata,1))),axis=1)
    change = range(self.ndata)

    updatew1 = np.zeros((np.shape(self.weights1)))
    updatew2 = np.zeros((np.shape(self.weights2)))

    for n in range(niterations):

        self.outputs = self.mlpfwd(inputs)

        error = 0.5*np.sum((self.outputs-targets)**2)
        if (np.mod(n,100)==0):
            print "Iteration: ",n, " Error: ",error    

        # Different types of output neurons
        if self.outtype == 'linear':
            deltao = (self.outputs-targets)/self.ndata
        elif self.outtype == 'logistic':
            deltao = self.beta*(self.outputs-targets)*self.outputs*(1.0-self.outputs)
        elif self.outtype == 'softmax':
            deltao = (self.outputs-targets)*(self.outputs*(-self.outputs)+self.outputs)/self.ndata 
        else:
            print "error"

        deltah = self.hidden*self.beta*(1.0-self.hidden)*(np.dot(deltao,np.transpose(self.weights2)))

        updatew1 = eta*(np.dot(np.transpose(inputs),deltah[:,:-1])) + self.momentum*updatew1
        updatew2 = eta*(np.dot(np.transpose(self.hidden),deltao)) + self.momentum*updatew2
        self.weights1 -= updatew1
        self.weights2 -= updatew2

        # Randomise order of inputs (not necessary for matrix-based calculation)
        #np.random.shuffle(change)
        #inputs = inputs[change,:]
        #targets = targets[change,:]

def mlpfwd(self,inputs):
    """ Run the network forward """

    self.hidden = np.dot(inputs,self.weights1);
    self.hidden = 1.0/(1.0+np.exp(-self.beta*self.hidden))
    self.hidden = np.concatenate((self.hidden,-np.ones((np.shape(inputs)[0],1))),axis=1)

    outputs = np.dot(self.hidden,self.weights2);

    # Different types of output neurons
    if self.outtype == 'linear':
        return outputs
    elif self.outtype == 'logistic':
        return 1.0/(1.0+np.exp(-self.beta*outputs))
    elif self.outtype == 'softmax':
        normalisers = np.sum(np.exp(outputs),axis=1)*np.ones((1,np.shape(outputs)[0]))
        return np.transpose(np.transpose(np.exp(outputs))/normalisers)
    else:
        print "error"

def confmat(self,inputs,targets):
    """Confusion matrix"""

    # Add the inputs that match the bias node
    inputs = np.concatenate((inputs,-np.ones((np.shape(inputs)[0],1))),axis=1)
    outputs = self.mlpfwd(inputs)

    nclasses = np.shape(targets)[1]

    if nclasses==1:
        nclasses = 2
        outputs = np.where(outputs>0.5,1,0)
    else:
        # 1-of-N encoding
        outputs = np.argmax(outputs,1)
        targets = np.argmax(targets,1)

    cm = np.zeros((nclasses,nclasses))
    for i in range(nclasses):
        for j in range(nclasses):
            cm[i,j] = np.sum(np.where(outputs==i,1,0)*np.where(targets==j,1,0))

    print "Confusion matrix is:"
    print cm
    print "Percentage Correct: ",np.trace(cm)/np.sum(cm)*100

 A: In formula $4.9$:
 $$\delta_h(\zeta) = \bar{g}(h_\zeta)\sum_{k=1}^N w_\zeta \delta_o(k),\text{ where }\bar{g}=g(h_\zeta)(1-g(h_\zeta)) ,$$ 
we have a $\zeta$-th neuron in a hidden layer which is connected to $N$ output neurons with weights $w_\zeta$. $\delta_o(k)$ is the error on the output neuron $k$, $w_\zeta$ is a vector of weights and it has length $N$.
So $$\sum_{k=1}^N w_\zeta \delta_o(k)= w_{\zeta1}*\delta_o(1)+w_{\zeta2}*\delta_o(2)...+w_{\zeta N}*\delta_o(N) $$
Each $\delta_o(k)$ is a vector of some length $O$, depending on the shape of the output.
$$\sum_{k=1}^N w_\zeta \delta_o(k)=w_{\zeta1}  \begin{pmatrix}\delta_{o1}(1)\\\delta_{o2}(1)\\ ...\\\delta_{oO}(1) \end{pmatrix}+w_{\zeta2}  \begin{pmatrix}\delta_{o1}(2)\\\delta_{o2}(2)\\ ...\\\delta_{oO}(2) \end{pmatrix} +...+w_{\zeta N}  \begin{pmatrix}\delta_{o1}(N)\\\delta_{o2}(N)\\ ...\\\delta_{oO}(N) \end{pmatrix}  $$
This is a vector matrix multiplication (note the tranposition) :
$$
\begin{pmatrix}w_{\zeta1}&w_{\zeta2}&... w_{\zeta N}  \end{pmatrix} \begin{pmatrix}\delta_{o1}(1)&\delta_{o2}(1)&...&\delta_{oO}(1)\\
\delta_{o1}(2)&\delta_{o2}(2)&...&\delta_{oO}(2)\\
...\\
\delta_{o1}(N)&\delta_{o2}(N)&...&\delta_{oO}(N)\\\end{pmatrix} $$
There are $M$ neurons in the hidden layer (plus a bias neuron)
$$
   \begin{pmatrix}
    w_{11}  & w_{12} & ... &w_{1N}  \\
    w_{21} & w_{22} & ... &w_{2N}  \\
    ...  \\
    w_{M+1,1} & w_{M+1,2} &... & w_{M+1,N} \\
    \end{pmatrix} \cdot \begin{pmatrix}\delta_{o1}(1)&\delta_{o2}(1)&...&\delta_{oO}(1)\\
\delta_{o1}(2)&\delta_{o2}(2)&...&\delta_{oO}(2)\\
...\\
\delta_{o1}(N)&\delta_{o2}(N)&...&\delta_{oO}(N)\\\end{pmatrix} $$
Thus we have $(M+1 \times N) \cdot (N\times O)=(M+1)\times O$ resulting matrix.
Or we could represent it as a matrix vector multiplication:
$$
\begin{pmatrix}\delta_{o1}(1)&\delta_{o1}(2)&...&\delta_{o1}(N)\\
\delta_{o2}(1)&\delta_{o2}(2)&...&\delta_{o2}(N)\\
...\\
\delta_{oO}(1)&\delta_{oO}(2)&...&\delta_{oO}(N)\\\end{pmatrix}
\begin{pmatrix}w_{\zeta1}&w_{\zeta2}&... w_{\zeta N}  \end{pmatrix}^T= $$
$$
=\begin{pmatrix}\delta_{o1}(1)&\delta_{o2}(2)&...&\delta_{o1}(N)\\
\delta_{o2}(1)&\delta_{o2}(2)&...&\delta_{o2}(N)\\
...\\
\delta_{oO}(1)&\delta_{oO}(2)&...&\delta_{oO}(N)\\\end{pmatrix}    \begin{pmatrix}
    w_{11}  & w_{12} & ... &w_{1N}  \\
    w_{21} & w_{22} & ... &w_{2N}  \\
    ...  \\
    w_{M+1,1} & w_{M+1,2} &... & w_{M+1,N} \\
    \end{pmatrix}^T
$$
getting $(O \times N) \cdot (N \times M+1) = O \times (M+1)$ matrix.
