\chapter{Mean Field Annealing}

As a prelude to mean field annealing (MFA), a summary of the
continuous Hopefield model is given below.

\section{The Hopfield Net - A summary of the continuous model}

\begin{itemize}
\item The network is built by connecting a large number of neurons to
each other to form a dense, parallel array. In general, any $i^{th}$
neuron in the net is described by two variables: its current state
(input) $u_i$,  and its output $v_i$.

The output is related to the input by a monotonically increasing
transfer function: $v_i = g_{\lambda}(u_i)$.

\item The net's dynamics are governed by the first order differential
equation:
\be{FODE.1}
{{\bf du} \over {\bf dt}} = 
- \eta {\bf u} + {\bf W v} + {\bf I^b}\;.
\ee

\item $g_{\lambda}(u_i) = {1 \over {1 + e^{-2\lambda u_i}}}\;$, for
large $\lambda$, this function approximates a step function,
forcing $v_i$ to $[0,1]\;\Longrightarrow\;$ (the output of the net
lies on the unit hypercube).

\item $W_{ij}$ is the connection strength between neuron $i$ and neuron
$j$. 

\item $I_i^b$ is an input bias to neuron $i$.

\item Neurons may be updated individually at random (asynchronously)
or all together (synchronously).
\end{itemize}

\subsection{Operation of the Network}

\vspace{0.1cm}
\begin{figure}[tbp]
\begin{center}
\begin{picture}(320,325)
\psfig{figure=hopfield_net.eps,height=320pt,width=325pt}
\end{picture}
\end{center}
\caption{A schematic view of the Hopfield network.}
\label{Hop-model}
\end{figure}

The network (see Figure \ref{Hop-model}) operates as follows:

\begin{itemize}
\item The network's dynamics are governed by a set of nonlinear
coupled differential equations of the form:

\be{FODE.2}
{{du_i} \over {dt}} = - u_i + \sum_{j=1}^{n} W_{ij} v_j + I_{i}^{b}\;.
\ee

or the general form of it:

\be{FODE.3}
{d {\bf u} \over {dt}} = - \eta {\bf u} + {\bf W v} + {\bf I^b}\;.
\ee

where {\bf u}, {\bf W}, {\bf v}, and {\bf I} are vectors.

\item The $I_{i}^b$ corresponds to an input bias to neuron $i$.
\item The $\eta$ term introduces an element of decay into the
network's dynamics, which can be used to ensure that ${\bf u}$
variables remain bounded in magnitude.
\end{itemize}

\subsection{Updating of the States of the Neurons}

To make sure that stability of the of the net is guaranteed,
$u_i$ is updated continuously by evaluating ${d \over dt}u_i$ as a
function of the total input to the $i^{th}$ neuron.

\be{State.1}
{\bf \dot{u}} = 
-\eta{\bf u}+{\bf W v}+{\bf I^{b}}\;,
\ee

where 

\be{State.2}
v = g(u)\;.
\ee

Equations \ref{State.1} and \ref{State.2} govern the dynamics of the
network.

\subsection{Stability of the Net}

Given that {\bf W} is symmetric, Hopfield proved stability for the
network by showing that with the above dynamic equations
(\ref{State.1} and \ref{State.2}) the network output {\bf v} evolves
so as to minimize the {\it Liapunov function}:

\be{Liapunov}
E^{liap}(v) = 
\left( -{1 \over 2} \sum_{ij} W_{ij} v_i v_j -\sum_i I_{i}^b v_i \right) 
+ \eta{\displaystyle \sum_i}\, {\displaystyle \int_{0}^{v_i}}\, g^{-1}(V)\,dV\;.
\ee

\begin{itemize}
\item $v = (v_1 \ldots v_n)$ is denoting the net state, therefore we
have an energy function defined throughout $v \in [0,1]^n$ and
dynamics which guarantee convergence to minima.

\item For a suitable $W$, the minima of the first term lie at the
hypercube vertices. The second term is minimized at the hypercube
center but has negligible impact at low $\eta$.
\end{itemize}

\subsection{Updating the States of the Net}
Updating of the states (or the output of the neurons) is carried
out either asynchronously or synchronously. From a simulation point of
view, we found that synchronous updating a bit more efficient, but
from a theoretical point of view, it does not make much of a 
difference on the operation of the network as to which updating mode
is used.

\subsection{Mapping the Problem onto the Network}
A representation of the problem is needed in which the feasible
solutions lie at hypercube vertices. Each hypercube vertex is regarded
as a configuration of the problem, with an associated energy (or cost)
function. Next, the partition function is then formed by summing the
Boltzmann factors of the admissible configurations. A quadratic energy
function (to be mapped onto the first term of \ref{Liapunov}) for
which the minima correspond to solutions and the depth of each minimum
reflects the solution quality.

Hopfield and Tank used the TSP problem as an example, which they have
mapped onto the net using a permutation matrix representation of the TSP
tours.

\subsection{Shortcomings of the Hopfield Model}
The net's partition function sums over a vast number of
configurations (e.g. for the TSP, all of the ${2^{N}}^2$
configurations are admissible). These configurations, again, for the
TSP, are nothing like valid tours, and even though these offenders
have small Boltzmann weights their large number inevitably affects the
thermal average quantities. 


\section{MFA: The Potts Approach for Timetabling}

\subsection{Introduction and Motivation}
One of the potential drawbacks of using simulated annealing for hard 
optimization problems is that finding a good solution can often take an 
unacceptably long time. 
Mean-field annealing (MFA) attempts to avoid this problem by using a
deterministic approximation to simulated annealing, by attempting to average
over the statistics of the annealing process. 
The result is improved execution speed at the expense of solution quality. 
Although not strictly a continuous descent technique, MFA is closely related to
the Hopfield neural network, on which we would like to state few more
points. We clearly know that the original Hopfield algorithm had
difficulty constraining the network into valid solutions. Therefore,
one possible way forward is to reduce the computational burden being
placed on the net. For example, in the case of the TSP, this can be
done by constraining each city to be ``on'' only once, by enforcing
$\sum_a V_{ia} = 1\; \forall i$, ($V_{ia} = 1/0$, which means that the
city $i$ is (is not) the $a^{th}$ city visited in the 
tour), rather than relying on an energy penalty term to try to
explicitly do this, as is the case in Hopfield model. 

This idea was first utilized and analyzed for the TSP 
case by, among others, Peterson and S$\ddot{o}$derberg
(ref: Peterson-89). It greatly enhances the degree of convergence to 
valid solutions. Later, Peterson and colleagues somewhat extended the
technique and applied it to a case of high school timetabling
(ref: Peterson-91).
In (ref: Peterson-89) the normalization (definition \ref{norm-dfn})
occurs in the context of a {\it mean field annealing} approach. This
generally yields the same final network equations as the ``neuronal
circuit'' approach of Hopfield, but the exposition is a little bit
clearer, as it is laid out  in a statistical mechanics
framework\footnote{This equivalence was also noted in
Hopfiled and Tanks's 1985 paper (ref: Hop-85). Martin Simmen brought 
this point to our attention.}. The idea is to regard each hypercube
vertex as a configuration, with an associated energy (or cost) ${\it
E}$. The partition function is then formed by summing the Boltzmann
factors of the  admissible configurations. For example, in the
Hopfield algorithm, for the TSP, all of the $2^{N^2}$ configurations
are admissible, whereas if each city is restricted to being visited
only once, then only $N^N$ vertices are admissible. After taking a
mean field approximation, saddlepoint equations are derived (see next
sections), the solutions of which pick out the dominant states of the
network at the current temperature $T$.

This statistical mechanics characterization gives a somewhat clearer
picture of why the normalized model (which corresponds
to a type of Potts model (ref: Wu-82) in physics, in that each
city represented by a spin which can be in only one of $N$ states)
ought to perform better than Hopfield model. As indicated in
(ref: Peterson-89), the partition function of Hopfield algorithm sums
over a vast number of configurations which are nothing like valid
tours in the TSP, and even though these offenders have small
Boltzmann weights their large number inevitably affects the thermal
average quantities. The other advantage that draws us to this
derivation is that, by introducing a temperature, the concept of
annealing, i.e. reducing $T$ during a run, can be justifiably
employed. 

Th Mean-field annealing algorithm has been successfully applied to, as
mentioned above, high school class scheduling. For scheduling, it is
advantageous to use a  Potts neural encoding to specify discrete
neural variables (or neurons) for the problem. This is defined in its
simplest form as a mapping of events onto space-time slots, for
example an event $i$, in this case a professor-class pair ($p,q$), is
mapped onto a space-time slot $a$, in this case a classroom-timeslot
pair ($x,t$).  Now, the Potts neurons $S_{ia}$ are defined to be $1$
if event $i$ takes place in space-time slot $a$, and $0$ otherwise.
In this way, the constraints involved can be embedded in the neural
net in  terms of the weights $w_{i,j}$ of the network, which encode a 
Potts normalization condition such as $\sum_{a} S_{ia} = 1$.

For a full derivation of the mean-field
annealing algorithm from its roots in statistical physics, see 
Krogh and Hertz (ref: Hertz-91).

\subsection{Brief Overview}

The basic idea is that it is possible to approximate the actual cost 
or energy function $E$, which is a function of discrete neural
variables $S_{ia}$, by an effective energy function $E^{'}$ that can
be represented in terms of continuous variables $U_{ia}$ and 
$V_{ia}$. These are known as mean field variables, since $V_{ia}$ is
an approximation to the average value of $S_{ia}$  at a given
temperature $T$.

This approach effectively smooths out the energy function and makes it 
easier to find the minimum value, which is obtained by solving 
the saddle point equations 
${\partial E^{'} \over \partial V_{ia}} = 0$ and  
${\partial E^{'} \over \partial U_{ia}} = 0$, 
which generate a set of self-consistent mean field theory (MFT) equations 
in terms of the mean field variables $U$ and $V$:

\be{spe-1}
U_{ia} = - {1 \over T} {\partial E \over \partial V_{ia}}
\ee

\be{spe-2}
V_{ia} = {e^{U_{ia}} \over \sum_b e^{U_{ib}}}. 
\ee

The MFA algorithm involves solving equations \ref{spe-1} and
\ref{spe-2} at a series of progressively lower temperatures $T$: this
process is known as temperature annealing. The critical temperature
$T_c$, which sets the scale of $T$, is estimated by expanding equation
\ref{spe-2} around the trivial fixed-point (ref: Peterson-89,
Peterson-91) $V^{(0)}_{ia}
= {1 \over N_a}$, where $N_a$ is the number of possible states of each
of the network neurons. For example, for the events defined by
professor-class pairs $(p,q)$ mapped onto classroom-timeslots $(x,t)$,
we have $N_p N_q$ neurons, each of which has $N_x N_t$ possible
states, in which case $V^{(0)}_{pq;xt} = {1 \over N_x N_t}$. 

Equations \ref{spe-1} and \ref{spe-2} can be solved iteratively using 
either synchronous or serial updating. The iterative dynamics to
evolve the mean field variables toward a self-consistent solution 
is explained in detail by Peterson at al. (ref: Peterson-91). 
The solutions correspond to stable states of the Hopfield
model (ref: Hopfield-85).

Observe from equation \ref{spe-2} that any solution to the MF
equations respects a continuous version of the Potts condition 

\be{potts-cond}
\sum_a V_{ia} = 1  \;\;\;\forall\;i .
\ee

\subsection{The Mean-Field Annealing Algorithm}

The generic MFA algorithm appears in Figure~\ref{MFA-algorithm}. 

At high temperatures $T$, the mean-field solutions will be states 
near the fixed-point symmetrical maximum entropy state $V_{ia} =
1/N_a$. At low temperatures, finding a mean-field solution will be
equivalent to using the Hopfield model, which is highly sensitive to
the initial conditions and known to be ineffective for hard
problems. MFA improves over the Hopfield model by using
annealing to slowly decrease the temperature in order to sidestep
these problems.
 
These characteristics are similar to those of simulated annealing,
which is no surprise since both it and the mean-field method compute
thermal averages over Gibbs distributions of discrete states, the
former stochastically and the latter through a deterministic
approximation. It is therefore natural to couple the mean-field method
with the concept of annealing from high to low temperatures. 

In addition to the structure of the energy function, there are three
major interdependent issues which arise in completely specifying 
a mean-field annealing algorithm for a timetabling problem:

\begin{itemize}
\item The values of the coefficients of terms in the energy function.
\item The types of dynamics used to find solutions of the MFT
 equations at each $T$.
\item The annealing schedule details, i.e. the initial temperature
 $T(0)$, the rules for deciding when to reduce $T$ and
 by how much, and the termination criteria.
\end{itemize}

Peterson {\it et al.} (ref: Peterson-91) introduced a quantity 
called {\it saturation}, $\Sigma$, defined as
\be{saturation}
\Sigma = {1 \over N_i} \sum_{ia} V^{2}_{ia} \;\; ,
\ee

where $N_{i}$ is the number of events (in this case the number of 
professor-class pairs).
This characterizes the degree of clustering of events in time and/or 
space, $\Sigma_{min} = {1 \over N_a}$ corresponds to high temperature,
whereas $\Sigma_{max} = 1$ means that all the $V_{ia}$ have converged
to 0 or 1 values, indicating that each event has been assigned to a
space-time slot.


\begin{figure}[tbp]
% \vspace{0.25cm}
% \hrule
\medskip
\begin{center}
\begin{enumerate}
\item Choose a problem and encode the constraints into weights $\{w_{ij}\}$.
\item Find the approximate phase transition temperature by linearizing
  equation (\ref{spe-2}).
\item Add a self-coupling $\beta$-term if necessary. In a neural net, 
  this corresponds to a feedback connection from a neuron to itself.
\item Initialize the neurons $V_{ia}$ to high temperature values
  ${1 \over N_a}$ plus a small random term such as ${rand[-1,1] \times 0.001}$; 
  and set $T(0) = T_c$.
\item Until ($\Sigma \geq 0.99$) do:
\begin{itemize}
  \item At each $T(n)$, update all $U_{ia}$ and $V_{ia}$ by iterating to
     a solution of the mean field equations.
  \item $T(n+1) = \alpha T(n)$, we chose $\alpha = 0.9$
\end{itemize}
\item The discrete values $S_{ia}$ that specify the schedule are obtained 
   by rounding the mean field values $V_{ia}$ to the nearest integer 
   (0 or 1). 
\item Perform {\it greedy heuristics} if needed to account for possible 
   imbalances or rule violations.
\end{enumerate}
\end{center}
% \smallskip
% \hrule
% \smallskip 
\caption{The Generic Mean-Field  Annealing Algorithm}
\label{MFA-algorithm}
% \vspace{0.25cm}
\end{figure}

\subsection{Observations on the Algorithm}

The first step of Figure~\ref{MFA-algorithm} is to map the constraints
of the problem into the neural net connection weights.
In our implementation, at each $T_n$ the MFA algorithm
(Figure~\ref{MFA-algorithm})  performs one update per neural variable 
(defined as one sweep) 
with sequential updating using equations \ref{spe-1} and \ref{spe-2}. 
After reaching a saturation value close to 1 (we chose $\Sigma = 0.99$)
we check whether the
obtained solutions are valid, i.e. $E_{hard} = 0$. If this is not the
case the network is re-initialized and is allowed to
resettle. We repeat this procedure a number of times until the best 
solution is found. 
A similar procedure was carried out on high school
scheduling by Peterson {\it et al.}

The MFA implementation was a little more complicated 
than the implementation of simulated annealing and the expert system, 
since it had many more parameters to handle, and it was 
often more difficult to find optimal values for these parameters.  
For example, one complication is the computation of the critical
temperature $T_c$, which involved an iterative procedure of a
linearized dynamic system.  On the other hand, we observed that the 
convergence time was indeed much less than any of the convergence
times of the simulated annealing using the three annealing schedules
studied.

\subsection{Definitions and Constraints}
Here we define all of the terms and concepts used throughout this
chapter for the implementation of class scheduling on the Potts
model. 

Let ($N_p$) to be a number of professors that give ($N_q$) classes in 
($N_x$) rooms at ($N_t$) time slots. Event-wise, $(p,q)$ event
-- denoting a professor $p$ giving a class $q$ -- takes place in the
space-time slot $(x,t)$ -- denoting room $x$ at time $t$.

\bd{dfn-0}
\ \ Let $S_{pq;xt}$ denote the Potts (or multi-state) neuron, where
\be{form-1}
S_{pq;xt} =
\twoif {1} {\text{if event $(p,q)$ takes place in $(x,t)$ ,}}
       {0} {\text{otherwise.}}
\ee
\ed

\bd{const-dfn.1}
\ \ We have dealt with two types of constraints:
\begin{enumerate}
\item Those involving professors, rooms, subjects, and time
periods.
\item Those dealing with students scheduling and preferences.
\end{enumerate}
In addition, some of these constraints are usually classified to be
{\bf hard}, others are {\bf medium} and the rest are {\bf soft}.
\ed

\bd{const-dfn.2}
Here, we define a number of {\bf hard constraints} of the first
type (see definition \ref{const-dfn.1}) .
\begin{itemize}
\item An event $(p,q)$ should occupy precisely one space-time slot
$(x,t)$.
\end{itemize}
This constraint can be embedded in the neural net in terms of the
Potts normalization condition:
\be{form-2}
\sum_{x,t} S_{pq,xt} = 1
\ee
\begin{itemize}
\item Different events $(p_1,q_1)$ and $(p_2,q_2)$ should not occupy
the same space-time slot $(x,t)$.
\end{itemize}
\be{form-3}
E_{XT} =
{1 \over 2}\sum_{x,t}\sum_{p_1,q_1}\sum_{p_2,q_2}S_{p_1q_1;xt}S_{p_2q_2;xt} =  
{1 \over 2}\sum_{x,t}{\left(\sum_{pq}S_{pq;xt}\right)}^2
\ee
\begin{itemize}
\item A teacher $p$ should have at most one class at a time.
\end{itemize}
\be{form-4}
E_{PT} =
{1 \over 2}\sum_{p,t}\sum_{q_1,x_1}\sum_{q_2,x_2}S_{pq_1;x_1t}S_{pq_2;x_2t} = 
{1 \over 2}\sum_{p,t}{\left(\sum_{qx}S_{pq;xt}\right)}^2
\ee
\begin{itemize}
\item A class $q$ should have at most one teacher at a time.
\end{itemize}
\be{form-5}
E_{QT} =
{1 \over 2}\sum_{q,t}\sum_{p_1,x_1}\sum_{p_2,x_2}S_{p_1q;x_1t}S_{p_2q;x_2t} = 
{1 \over 2}\sum_{q,t}{\left(\sum_{px}S_{pq;xt}\right)}^2
\ee
\ed

In (ref: Hop-84, Hop-85), Hopfield introduced a model with neurons
having a graded response. So given $I_i$ a fixed input bias to neuron
$x_i$, $W_{ij}$ the weight of the connection between neurons
$(x_i,x_j)$, and parameter $T$ to denote the ${\it Temperature}$, the
behavior of $x_i$  can be stated as:

\be{form-JH}
 x_i =
 {1 \over 2} + {1 \over 2} \cdot 
 tanh\left(({\sum_{i \neq j}W_{ij}x_j - I_i}) \times {1 \over T}\right)
\ee
When $T \longrightarrow 0$, function ({\ref{form-JH})'s
response takes on a step-wise form. These models of (ref: Hop-85)
are known to be a bit sensitive to parameter tuning for real world
applications, especially those of the type that involves ``finding''
or ``searching for'' $1$ out-of $N$ type problems.  

\bd{don-1}
\ \ At a finite $T$, $x_i$ (of eq. \ref{form-JH})
can be replaced by an MFT variable $\upsilon_{ia}$, then
differentiating to obtain the Potts MFT equations:

\be{form-5.1}
\upsilon_{ia} =
- {1 \over T} \cdot {{\partial E}[{\vec\nu}] \over {\partial \nu_{ia}}}
\ee
\be{form-5.2}
\nu_{ia} =
{e^{\upsilon_{ia}} \over {\sum_b e^{\upsilon_{ib}}}}
\ee

Which enforces

\be{form-5.3}
\sum_{a} \nu_{ia} = 1\,\,\,\, \forall i\, .
\ee
Where $a$ denotes the vector components or all the elements of the
group. Also, we could state that $\nu_{ia}$ $\approx$
Probability(neuron(i) $\in$ state(a)). Equations (\ref{form-5.1},
\ref{form-5.2}) are computed in an iterative manner.
\ed

\bd{dfn-2}
\ \ To carry out the minimization of the energy functions of
definition (\ref{dfn-0}) we need to convert the neural variables
$\{0,1\}$ of $S_i$ to their corresponding {\it mean field variables} 
at $T$. These correspondences are derived (see (ref: Peterson-92)) to
be:

\be{form-6.1}
V_{pq;xt} \approx <S_{pq;xt}>_T\,\,;
\ee
which is the thermal average of the original binary spines, with the
following multi-state $MFT$ equations (see (\ref{form-5.1},
\ref{form-5.2}) for a general format):

\be{form-6}
U_{pq;xt} =
- {1 \over T} \cdot {{\partial E} \over {\partial V_{pq;xt}}}
\ee
\be{form-7}
V_{pq;xt} =
{e^{U_{pq;xt}} \over {\sum_{x^{'}t^{'}}e^{U_{pq;x^{'}t^{'}}}}}
\ee
\ed

The $MFT$ equations of $\upsilon$ and $\nu$ are 
solved {\it iteratively} either synchronously or asynchronously, under 
$T$. Also, these equations give rise to 2 phase  transitions; one in
$x$ (space) and the other in $t$ (time). This leads to factorization
of the system by replacing $S_{pq;xt}$ by {\bf $X$-neurons}
$S^{(X)}_{pq;x}$ and {\bf $Time$-neurons} $S^{(Time)}_{pq;t}$ :
 
\bd{dfn-3} \ \
\be{form-8}
S_{pq;xt} =
S^{(X)}_{pq;x} \cdot S^{(Time)}_{pq;t}
\ee
With separate Potts conditions:

\be{form-9}
\sum_{x}S^{(X)}_{pq;x} = 1
\ee
\be{form-10}
\sum_{t}S^{(Time)}_{pq;t} = 1
\ee

For {\bf $X$-neurons} and {\bf $Time$-neurons}, respectively. 
Rather than dealing with $p$ and $q$ directly, an independent variable
($i$) -acting as an event index- to which $p$ and $q$ are attributes
$p(i)$ and $q(i)$. The value of these exits in a look-up table which
contains all the relevant information to process each event ($i$).
Now, the above (\ref{form-9}, \ref{form-10}) conditions can be 
re-stated as follows:

\be{form-9.1}
\sum_{x}S^{(X)}_{i;x} = 1
\ee
\be{form-10.1}
\sum_{t}S^{(Time)}_{i;t} = 1
\ee

\ed

During annealing, as we go from a higher temperature to a lower one,
a {\it phase transition} is passed at ($T = T_c$), where 
$T_c$\footnote{See pp 29-31, (ref: Hertz-91) for an approximate
value in the case of Ising ferromagnets.} is the {\it critical
temperature}. The value of $T_c$ set the scale or the range of the
value of $T$, and it is estimated by expanding \ref{form-7} around
some trivial fixed-point such as

\be{trivial-form}
V^{(0)}_{pq;xt} = {1 \over {N_x N_t}} \;\text{, at $T$ = 0} .
\ee

Also $T_c$ depends on the self-coupling coefficient $\beta$, and on
the method of the update.  This parameter $\beta$ acts as a
Lagrangian multiplier for the given energy constraint(s). See section 
on the generic Potts algorithm and (ref: Peterson-89).

\bd{dfn-4}\ \
{\em Relative Clamping}: During the revision process of an existing
schedule, sometimes it is desirable to have a certain cluster of
events to stick together in time (or within some specified period of
time). This amounts to making the following notational replacement: 
\be{form-12}
i \longrightarrow (j,k) \;\; ;
\ee
where $(j)$ denotes the cluster and $(k)$ an event within $(j)$. In this
situation, one has a common $t$-neuron $(S^{(Time)}_{j;t})$ for cluster
$j$ but distinct $x$-neurons $(S^{(X)}_{jk;x})$ for the individual events.
\ed

\bd{dfn-3.9}
To deal with events that take place over two consecutive hours, what
(ref: Peterson-89) called effective $time$-neurons
$\tilde{S}^{(Time)}_{i;t}$ are introduced: 
\be{form-11}
\tilde{S}^{(Time)}_{i;t} =
\sum_{k=0}^{g_i - 1} S^{(T)}_{i;t-k} \;\;{\text{;}}
\ee
where $g_i = 1$ for single hours and $= 2$ for double hours. For those 
events that occur between one and two hours, $g_i$ can be adjusted
accordingly to handle these cases.
\ed

Utilizing definition (\ref{dfn-3.9}), constraints (\ref{form-3},
\ref{form-4}, \ref{form-5}) can be put into the following form:

\bd{dfn-3.9.1}

\be{form-11.1}
 E_{XT} =
 {1 \over 2}\sum_{x,t}\sum_{i \neq i^{'}}
 S^{(X)}_{i;x}S^{(X)}_{i^{'};x} \tilde{S}^{(Time)}_{i;t}
 \tilde{S}^{(Time)}_{i^{'};t}
\ee

\be{form-11.2}
 E_{PT} =
 {1 \over 2}\sum_{p,t} \sum_{i \neq i^{'}} \delta_{p(i),p(i^{'})}
 \tilde{S}^{(Time)}_{i;t} \tilde{S}^{(Time)}_{i^{'};t}
\ee

\be{form-11.3}
 E_{QT} =
 {1 \over 2}\sum_{q,t} \sum_{i \neq i^{'}} \delta_{q(i),q(i^{'})}
 \tilde{S}^{(Time)}_{i;t} \tilde{S}^{(Time)}_{i^{'};t}
\ee
\ed

\bd{dfn-6}\ \
The Kronecker $\delta_{q(i),q(i')}$ ensures that only
events $i$ and $i'$ with identical ``primordial'' classes are,  for
example, summed over. Here, $q(i)$ and $q(i')$ refer to classes of the
same type or 2 sections derived from the same original class. 

\be{Kronecker}
\delta_{r,s} =
\twoif {1} {\text{if} r = s}
       {0} {\text{otherwise.}}
\ee
\ed

\bd{dfn-3.9.2}\ \
One important soft constraint that need to be handled is {\it group
formation}. Briefly, it can be stated as follows: for some optional
subjects the classes are broken up into option groups temporarily
forming new classes. In order to account for this, we need a formalism
that allows for the breaking up of these primordial classes and
subsequent recombinations into {\it quasi}-classes. This is a common
practice with those classes where the students have many choices. 
To carry out such an extension to quasi-classes one might
have contributing pairs of events where the primordial
classes are different. Hence one should replace $\delta$ with
a more structured overlap matrix $\Gamma$ which is given by a look-up 
table.
\be{form-11.4}
\delta_{q(i),q(i^{'})} \longrightarrow \Gamma_{q(i),q(i^{'})} \; .
\ee
\ed

\bd{dfn-3.9.3}\ \
From definitions ($\ref{dfn-4}$, $\ref{dfn-3.9.2}$), we need to make
sure that ($\forall k$) events within a given ($j$) cluster 
$p(j,k) \neq 0$. This says that there can be no ``uncovered'' class
period or a class event (that is, periods when no professor has been
assigned to an already scheduled class).
For this purpose we introduce $\tilde{\delta}_{j,k}$
according to:
\be{form-3.9.3}
\tilde{\delta}_{j,k} = \delta_{j,k}(1 - \delta_{j,0}\delta_{k,0}) \; .
\ee
If $j = 0$ then $\forall k$ , $k = 0$ .
\ed

\bd{dfn-3.9.4}\ \
Utilizing the concepts of the above definition ($\ref{dfn-3.9.2}$,
$\ref{dfn-3.9.3}$), the collision terms of definition
($\ref{dfn-3.9.1}$) are extended to ensure that all the hard
constraints are satisfied (or the solutions obtained are legal).
The extension takes the following forms:

\be{form-3.9.4.1.1}
 E_{XT}{\text{(form 1)}} =
 {1 \over 2}\sum_{jk \neq j^{'}k^{'}}\sum_{x}
 S^{(X)}_{jk;x}S^{(X)}_{j^{'}k^{'};x}
\sum_{t}\tilde{S}^{(Time)}_{j;t}\tilde{S}^{(Time)}_{j^{'};t}
\ee

\be{form-3.9.4.1.2}
 E_{XT}{\text{(form 2)}} =
 {1 \over 2}\sum_{j \neq j^{'}}\sum_{k k^{'}}\sum_{x}
 S^{(X)}_{jk;x}S^{(X)}_{j^{'}k^{'};x}
 \sum_{t}\tilde{S}^{(Time)}_{j;t}\tilde{S}^{(Time)}_{j^{'};t}
 + {1 \over 2}\sum_{j}g_j \sum_{k \neq k^{'}} \sum_{x} 
 S^{(X)}_{jk;x}S^{(X)}_{jk^{'};x}
\ee

\be{form-3.9.4.2.1}
 E_{PT}{\text{(form 1)}} =
 {1 \over 2}
\sum_{p,t} \left[ \sum_{jk} \tilde{\delta}_{p,p(j,k)}
 \tilde{S}^{(Time)}_{j;t} \right]
\left[ \sum_{j^{'}\neq j}\sum_{k^{'}}
\tilde{\delta}_{p, p(j^{'},k^{'})} 
 \tilde{S}^{(Time)}_{j^{'};t} \right]
\ee

\be{form-3.9.4.2.2}
 E_{PT}{\text{(form 2)}} =
        {1 \over 2}\sum_{j \neq j^{'}}\sum_{k k^{'}} 
        \tilde{\delta}_{p(j,k), p(j^{'},k^{'})}
        \sum_{t}\tilde{S}^{(Time)}_{j;t} \tilde{S}^{(Time)}_{j^{'};t}
\ee

\be{form-3.9.4.3.1}
 E_{QT}{\text{(form 1)}} =
 {1 \over 2} 
\sum_{qt} \left[ \sum_{jk} \delta_{q,q(j,k)} \tilde{S}^{(Time)}_{j;t}
\right] 
\left[ \sum_{j^{'} \neq j} \sum_{k^{'}} \Gamma_{q,q(j^{'}, k^{'})}
        \tilde{S}^{(Time)}_{j^{'};t} \right]
\ee

\be{form-3.9.4.3.2}
 E_{QT}{\text{(form 2)}} =
        {1 \over 2}
        \sum_{j \neq j^{'}}\sum_{kk^{'}} 
        \Gamma_{q(j,k), q(j^{'},k{'})} \sum_{t} \tilde{S}^{(Time)}_{j;t}
        \tilde{S}^{(Time)}_{j^{'};t}
\ee

Note: Each ``form 2'' above is a derivation from its corresponding
``form 1''. 
\ed

The factorization idea of (dfn \ref{dfn-3.9.4}) can be extended to
{\it mean field} variables
$V^{(X)}_{jk;x}$ and $V^{(Time)}_{j;t}$, corresponding to neurons
$(S^{(X)}_{jk;x})$ and $(S^{(Time)}_{j;t})$, respectively.
$V^{(X)}_{jk;x}$ and $V^{(Time)}_{j;t}$ can be interpreted as being 
the probabilities that the corresponding events occur
at given $x$-value and $t$-value with normalization given by the above
equations (\ref{form-9.1}, \ref{form-10.1}) .

Substituting $V$ for $S$ in (\ref{form-9.1}, \ref{form-10.1}), 
the {\it mean field} equations at a finite 
$T$, are given in terms of the {\it local fields} 
$U^{(X)}_{jk;x}$ and $U^{(Time)}_{j;t}$. Now, the {\it mean field} 
dynamics proceed over the following local and mean field variables:

\bd{dfn-3.99}
\be{form-4.1.3}
U^{(X)}_{jk;x} =
- {1 \over T} \cdot {{\partial E} \over {\partial V^{(X)}_{jk;x}}}
\ee

\be{form-4.1.4}
U^{(Time)}_{j;x} =
- {1 \over T} \cdot {{\partial E} \over {\partial V^{(Time)}_{j;t}}}
\ee

\be{form-4.1.5}
V^{(X)}_{jk;x} =
{e^{U^{(X)}_{jk;x}} \over {\sum_{x^{'}} e^{U^{(X)}_{jk;x^{'}}}}}
\ee

\be{form-4.1.6}
V^{(Time)}_{j;t} =
{e^{U^{(Time)}_{j;t}} \over {\sum_{t^{'}} e^{U^{(Time)}_{j;t^{'}}}}}
\ee

The iterations of the {\it MF} dynamics over equations 
({\ref{form-4.1.3}} - {\ref{form-4.1.6}})
can either be {\bf synchronous} or {\bf serial} (see (ref: Peterson-92)).
\ed

\bd{dfn-4.2}
{\em Synchronous Updating}: \ \
A single sweep consists of first computing all the local fields, and
then updating all the neurons. That means do equations
({\ref{form-4.1.3}},{\ref{form-4.1.4}}) then perform equations
({\ref{form-4.1.5}},{\ref{form-4.1.6}}). 
\ed

\bd{dfn-4.3}
{\em Serial Updating}: \ \
For each neuron $S^{(X)}_{jk;x}$ and/or $S^{(Time)}_{j;t}$, the {\it
local field} is computed immediately before the corresponding {\it
neuron state} is updated. That means either do equation
({\ref{form-4.1.3}}) followed by ({\ref{form-4.1.5}}) and/or do
equation ({\ref{form-4.1.4}}) followed by ({\ref{form-4.1.6}}).
\ed

Next, we need to define 2 important {\bf soft constraints} encountered 
when scheduling classes.


\bd{dfn-5}\ \
\begin{itemize}
\item {\em Spreading}: different classes in a particular area
(i.e. computer science) should be spread over the week in a (MWF) and
(TuTh) type of spread. 
\end{itemize}

First, assume we have access to a subject attribute $s(j,k)$. Second,
we define an effective ``day-neuron'' as follows  
($\tilde{S}^{(Time)}_{j;dh}$ is reduced to a $D$-type neuron by
summing over  all hours of the day):

\be{day-neuron}
S^{(D)}_{i;d} = \sum_h \tilde{S}^{(Time)}_{j;dh} = g_j \cdot
S^{(Time)}_{j;dh}
\ee

So the penalty term $E_{QSD}$ that spreads classes $Q$ in a particular
subject $S$ over different week-days $D$ :

\be{spread-1}
E_{QSD} = 
{1 \over 2} \sum_{qsd}{\left[\sum_{jk}\sum_{h} \delta_{q,q(j,k)}
\delta_{s,s(j,k)} \tilde{S}^{(Time)}_{jk;dh}\right]}^2
\ee

Refining the above a bit more by considering 2 different clusters
$(j)$ and $(\hat{j})$ of classes over $D$ : ($q$ and $s$ are not used
in the summation, since no change occurs to them) 

\be{spread-2}
E_{QSD} =
{1 \over 2} \sum_{d}\sum_{j \neq \hat{j}}\sum_{kh}\sum_{\hat{k}\hat{h}}
\delta_{q(j,k),q(\hat{j},\hat{k})}\delta_{s(j,k), s(\hat{j},\hat{k})}
\hat{S}^{(Time)}_{j;dh}\hat{S}^{(Time)}_{j;d\hat{h}}
\ee

Now for some day $(d)$ of the week:

\be{spread-3}
E_{QSD} =
{1 \over 2} \sum_{j \neq \hat{j}} \sum_{k\hat{k}} 
\delta_{q(j,k),q(\hat{j},\hat{k})}\delta_{s(j,k), s(\hat{j},\hat{k})}
\sum_d S^{(D)}_{j,d} S^{(D)}_{\hat{j},d}
\ee

\begin{itemize}
\item {\em Glueing}: The schedule should have as few ``holes'' as
possible; that is lessons should be glued together. 
\end{itemize}

The penalty term that rewards situations where $(d,h)$-events are
glued to $(d,h-1)$-events can be defined as follows:

\be{glueing-1}
E_{QDH} = -\alpha \sum_{q}\sum_{d,h}
\left[\sum_{j,k}\delta_{q,q(j,k)}\tilde{S}^{(Time)}_{j;dh}\right]
\left[\sum_{\hat{j} \neq j}\Gamma_{q,q(\hat{j},\hat{k})}
\tilde{S}^{(Time)}_{j;d(h-1)}\right]
\ee

The energy term
of (\ref{glueing-1}) can be summarized to say that: {\it if cluster
($j$) with events each symbolized by ($k$) occurs on hour ($h$) of day
($d$) then it would be quite desirable to have another event
($\hat{k}$) from a different cluster ($\hat{j}$) to occur on hour
($h-1$) of that day.} 

Next, since ($q$) is fixed we could re-write the above equation as
follows:

\be{glueing-2}
E_{QDH} = -\alpha\sum_{\hat{j} \neq j}\sum_{k,\hat{k}}\Gamma_{q(j,k),
q(\hat{j},\hat{k})}
\sum_{d,h>1}\tilde{S}^{(Time)}_{j;dh}\tilde{S}^{(Time)}_{\hat{j},d(h-1)}
\ee

The presence of the ($-$) sign in front of (\ref{glueing-1}, 
\ref{glueing-2}) indicates that such energy functions are to be minimized. 
Also, the ($\alpha$) in front of the the above forms governs the strength of 
this reward relative to the energies of equations
(\ref{form-3},\ref{form-4},\ref{form-5}). 
This ($\alpha$) takes on a value between zero and $1$.

Finding {\em legal} solutions with no collisions corresponds to
minimizing the {\em hard} energy:
\be{Hard}
E_{hard} = combine (E_{XT}, E_{PT}, E_{QT})
\ee
to zero. The {\em soft energy} is given by:
\be{Soft}
E_{soft} = combine (E_{QSD}, E_{QDH})\;.
\ee
The total energy $E$ to be minimized is then:
\be{total}
E = combine (E_{hard}, E_{medium}, E_{soft})\;.
\ee
\ed

\subsection{More Definitions and Constraints}

Let $q(i)$ represents a section of class $q$, $O_{q(i)}$ indicates 
that section $q(i)$ is open, and $C_{q(i)}$ is closed. Also, let 
$b(x,t)$ be a balanced interval for some time period $(t)$, $R$ is the
overall set of students, and an arbitrary given subset is  
$r \subseteq R$.

Here, we need to introduce the $R$-neuron: $S^{(R)}_{r(i);(d(t)d(x)}$
where $k= 1,2,\ldots$. Each student $r(i) \in R$ should have a
balanced schedule over the two-group of days $(MWF)$ and
$(TuTh)$. We will also use $r_s$ to denote a single student.

For $X$ and $T$ neurons, let $x(j)$ and $t(j)$ where $j= 1,2,\ldots, n$
stand for places and time slots (or events), respectively. 
Using such notation, an example of the Potts normalization will be:
$S^{(X)}_{p(i)q(i);x(j)}$ and $S^{(Time)}_{p(i)q(i);t(j)}$ respectively,
where $i=1,2,\ldots,n$ is a section index.

\bd{athletes}\ \
{\it The athletes case}: Let $\vec{r} \subseteq R$ denotes the group
of athletes  that are dealt with in this study. 
Let ($r^{*}$,$b$) be playing/training events scheduled during
the weekdays at (or during) times when academic (attending lectures)
events ($r^{+}$,$a$)  are also scheduled, or may possibly occur. Let
($d(x)$, $d(t)$) a place/time during the week at which any of these
events can occur. 

In addition to the above stated Potts neurons, we also will be using :

\begin{enumerate}
\item $S^{(\vec{r},X)}_{r^{*}b;d(x)}$ ,
$S^{(\vec{r},Time)}_{r^{*}b;d(t)}$ : It takes a value of ($1$) if
event ($r^*,b$) takes place in the place and time slot, respectively,
of a weekday $d$, and zero otherwise.

\item $S^{(\vec{r},X)}_{r^{+}a;d(x)}$ ,
$S^{(\vec{r},Time)}_{r^{+}a;d(t)}$ : It takes a value of ($1$) if
event ($r^+,a$) takes place in the place and time slot, respectively,
of a weekday $d$, and zero otherwise.

\item $S^{(\vec{r},X)}_{jk;d(x)}$ : It takes a value of ($1$) if event
($k$) of cluster ($j$) that is associated with the $x$-neuron of
$\vec{r}$ occurs in place $d(x)$, and zero otherwise. 

\item $S^{(\vec{r},X)}_{\hat{j}\hat{k};d(x)}$ : It takes a value of
($1$) if event ($\hat{k}$) of cluster ($\hat{j}$) that is associated
with the $x$-neuron of $\vec{r}$ occurs in place $d(x)$, and zero
otherwise. 

\item $\tilde{S}^{(\vec{r},Time)}_{j;d(t)}$ : 
It takes a value of ($1$) if cluster ($j$) that is associated with
the $t$-neuron of $\vec{r}$ occurs in time $d(t)$, and zero
otherwise.

\item $\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)}$ : 
It takes a value of ($1$) if cluster ($\hat{j}$) that is associated
with the $t$-neuron of $\vec{r}$ occurs in time $d(t)$, and zero
otherwise. 
\end{enumerate}

In the above list, the last 2 neurons act as what Peterson calls {\it
effective} $time$-neurons. Not only will they handle periods of single
hours but also variants of longer and shorter intervals.
\ed

\bd{time-int.dfn}\ \
Set ($\epsilon$) and ($c$) to some desired values depending on the
time period of the set of classes, now, the above {\it effective} 
$t$-neurons can be defined as follows:

\be{time-int.1}
\tilde{S}^{(\vec{r},T)}_{j;d(t)} = 
\sum^{g_j - c}_{k=0} S^{(\vec{r},T)}_{j;d(t-k+\epsilon)}
\ee

$\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)}$ can be defined in a
similar manner. If ($c = 1$) and ($g_i = 2$), then we are dealing with
double-hour intervals so we set ($\epsilon = 0$). If ($c = 1$) and
($g_i = 1$) then ($\epsilon = 0$) and

\be{time-int.2}
\tilde{S}^{(\vec{r},Time)}_{j;d(t)} = S^{(\vec{r},Time)}_{j;d(t)} .
\ee
\ed

\bd{consts-dfn}\ \
Let $S^{(R)}_{q(i)r(i);b(x,t)}$ for balancing the enrollment in a
multi-section courses. Section $q(i)$ of students $r(i)$ are balanced
over space-time $b(x,t)$. 

To enroll student $r(i)$ in a particular section $q(i)$ given the tag
$EC$ to indicate the status of the student,  
we have $S^{(X[Time])}_{q(i)r(i);(d(x)[d(t)])}|EC\in \{0,1\}$  is set
to$0$ if $EC$ is $0$, otherwise, it is set to $1$ and the action
$q(i)r(i)$ occurs during place $d(x)$ and/or time $d(t)$. 

Let $lin(c_i,c_j)$ representing the concept that course $c_i$ is
linked to some other course $c_j$ ; in this case, $c_i$ is a
pre-requisite to $c_j$. The degree of ``linkage'' varies between $0$
and $1$. It takes a value of $1$ if under all circumstances, course
$c_i$ must be taken before $c_j$. 

To deal with this link concept, we assumed in definition (\ref{dfn-5}) 
that we have an attribute $s(j,k)$ to represent a subject. For a
strong linkage, we will assume that we're dealing with subjects
that belong to the same cluster ($j$) (e.g. physics courses), so the 
link between various events of ($j$) is: $lin(s(j,k), s(j,\hat{k}))$.

Other degrees of linkage can involve subjects or courses from
different groups or clusters (e.g. physics and math courses, or 
engineering and math courses), and this connection is in the 
form: $lin(s(j,k), s(\hat{j},\hat{k}))$.

Let $\hat{x}(i)$ be a single location (or a seat) in class
$x(i)$ such that $i=1,2,\ldots, num$ ; where ($num$) is the 
number of class-rooms.
\ed

\subsubsection{Day-Class Balancing}

Each student of any particular section (i.e. $r(k)$) should have a
balanced schedule between ({\it MWF}) and ({\it TuTh}). The balancing
operation is accomplished by spreading as much as possible the courses
of the student over each of ({\it MWF}) and ({\it TuTh}).

Let $D$ be the set $\{\{MWF\}, \{TuTh\}\} = \{ D_1, D_2 \}$. Assume
the existence of the subject attribute $s(j,k)$. We introduce the
energy function $E_{QSD}$ that spreads the classes ($Q$) in a
particular field or subject ($S$) over the  set $\{D_1, D_2\}$.

\bd{balance-dfn.1}
The Potts condition of balancing $S^{({\bf R})}_{q(m)r(m);b(x,t)}$
(take a value of $1$ for balanced and $0$, otherwise) is accomplished
through one (or all) of these energies:

\be{balance.1}
E^{r}_{QSD} =
{1 \over 2}\sum_{j \neq \hat{j}}\sum_{k \hat{k}} 
\delta_{q(j,k),q(\hat{j},\hat{k})}
\delta_{s(j,k),s(\hat{j},\hat{k})}
\sum_{d_1,d_2} S^{(D)}_{j,d_1} S^{(D)}_{\hat{j},d_2}
\ee

\be{balance.2}
E^{r}_{QSD_1} =
{1 \over 2}\sum_{j_1 \neq \hat{j_1}}\sum_{k_1 \hat{k_1}} 
\delta_{q(j_1,k_1),q(\hat{j_1},\hat{k_1})}
\delta_{s(j_1,k_1),s(\hat{j_1},\hat{k_1})}\sum_{d_1}
\sum_{i \neq \hat{i}} S^{(D)}_{i,d_1} 
S^{(D)}_{\hat{i},d_1}
\ee


\be{balance.3}
E^{r}_{QSD_2} =
{1 \over 2}\sum_{j_2 \neq \hat{j_2}}\sum_{k_2 \hat{k_2}} 
\delta_{q(j_2,k_2),q(\hat{j_2},\hat{k_2})}
\delta_{s(j_2,k_2),s(\hat{j_2},\hat{k_2})}\sum_{d_2}
\sum_{i \neq \hat{i}} S^{(D)}_{i,d_2} 
S^{(D)}_{\hat{i},d_2}
\ee
\ed

\subsubsection{The Case of Student Athletes}
The main constraint of this case very much analogous to the {\bf hard}
constraint stated somewhere above : {\it different events should not
occupy the same space-time slot}.

\bd{athletic}

\be{athletic.1}
E^{\vec{r}}_{d(XT)} =
{1 \over 2}\sum_d\sum_{xt}\sum_{\dot{r}b}\sum_{\ddot{r}a} 
S^{(\vec{r},X[Time])}_{\dot{r}b;d(x)[d(t)]} 
S^{(\vec{r},X[Time])}_{\ddot{r}a;d(x)[d(t)]}
\ee

Where $(\dot{r}b)$ and $(\ddot{r}a)$ are some two different events.
After some further refinement, we obtain:

\be{athletic.2}
E^{\vec{r}}_{d(XT)} =
{1 \over 2}\sum_{j \neq \hat{j}}\sum_{k \hat{k}}\sum_{d(x)}
S^{(\vec{r},X)}_{jk;d(x)} S^{(\vec{r},X)}_{\hat{j}\hat{k};d(x)}
\sum_{d(t)} \tilde{S}^{(\vec{r},Time)}_{j;d(t)} 
\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)} + 
{1 \over 2} \sum_j g_j \sum_{k \neq \hat{k}}\sum_{d(t)}
S^{(\vec{r},X)}_{jk;x} S^{(\vec{r},X)}_{j\hat{k};x} 
\ee

Since it is always the case that athletic training facilities 
are kept separate from the academic class-rooms and lecture halls, 
a more refined energy function is obtained:

\be{athletic.3}
E^{\vec{r}}_{d(Time)} =
{1 \over 2}\sum_{j \neq \hat{j}}\sum_{d(t)}
\tilde{S}^{(\vec{r},Time)}_{j;d(t)} 
\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)}
\ee
\ed

\subsubsection{Eligibility for Enrollment}
\bd{The Eligibility Criteria}
Before going through the process of balancing, each student 
($r(i) = r_s$) eligibility for enrolling in a given class need to be
determined.
Let the student's eligibility criteria be denoted by ($EC$), using 
the notation stated above, we could derive the
following energy function ---note: if $EC$ is set to ``No'' then it
takes a $0$ value in the energy function, otherwise, it is a $1$ : 

\be{eligibility-criteria.1}
E_{QR} = 
{1 \over 2}\sum_{qr}\sum_{jk}\sum_{i}
\delta_{q,q(j,k)}\delta_{r,r(i)} S^{X}_{ij;x(i)}|({\it EC\in\{0,1\}})
\ee

\be{eligibility-criteria.2}
\sum_{i} O_{q(i)}=1
\Rightarrow \text{Test:}\left\{ EC(R) + \text{satisfy-constraints}(Hard,Soft)\right\}
\ee
\ed

\subsubsection{Course Pre-requisite and Subsequent Assignment}

Utilizing the link concept (${\bf lin(s,s')}$) between two different
subjects (see definition \ref{consts-dfn}), 
a general energy assignment can be put together as follow: 

\bd{pre-sub.dfn}\ \
Let ($i$) stands for the situation at hand (also can be viewed as 
a space-time event).

\be{pre-requisite.1}
E_{link(1)} =
 -\gamma \sum_{is}{\left[\sum_{jk} \delta_{s,s(j,k)}
\tilde{S}^{(X[Time])}_{ln(s,s');i}\right]}^{2}
\ee

\be{pre-requisite.2}
E_{link(2)} =
 -\gamma \sum_{s}\sum_{j \neq j'}\sum_{k \neq k'} \delta_{s(j,k),s(j',k')}
\sum_{i} \tilde{S}^{(X[Time])}_{ln(s(j,k),s(j',k')); i}
\tilde{S}^{(X[Time])}_{ln(s(j,k),s(j,k')); i}
\ee

\be{pre-requisite.3}
E_{link(3)} =
 -\gamma \sum_{j \neq j'}\sum_{k \neq k'} \delta_{s(j,k),s(j',k')}
\tilde{S}^{(X[Time])}_{ln(s(j,k),s(j',k')); i}
\tilde{S}^{(X[Time])}_{ln(s(j,k),s(j,k')); i}
\ee

If ($\gamma > 0$) then function minimization is the target here,
otherwise maximization will take place.

\ed

\subsubsection{Enroll in all parts  of a course}
There are a number course of different types that not only
have lecture periods but also associated sections that consist
of either labs or recitations or both. Signing up for such a course
entails enrolling in all associated sections of it. 

Let us assume the existence of some subject attribute $s(j,k)$
consisting of 2 components: $\dot{s}$ denoting lectures, and
$\ddot{s}$ stands for labs. Now the energy equation depicting
this assignment can be outlined as follow:

\be{enroll-in-all-parts.1}
E_{enroll(1)} =
{1 \over 2}\sum_{s}\sum_{jk}\sum_{i} \delta_{s,s(j,k)}\delta_{r,r(i)} 
\sum_{t} S^{(Time)}_{jk;t}
\ee

\be{enroll-in-all-parts.2}
E_{enroll(2)} =
{1 \over 2}\sum_{si} \sum_{t\neq t'}\left[\sum_{jk} 
\delta_{s,\dot{s}(j,k)} \delta_{r,r(i)}
S^{(Time)}_{jk;t}\right]\left[\sum_{jk'} \delta_{s,\ddot{s}(j,k')}
\delta_{r,r(i)} S^{(Time)}_{jk';t'}\right]
\ee

\be{enroll-in-all-parts.3}
E_{enroll(3)} =
{1 \over 2}\sum_{k \neq k'}\delta_{s,\dot{s}(j,k)}
\delta_{s,\ddot{s}(j,k')}
\sum_{t \neq t'} S^{(Time)}_{jk;t} S^{(T)}_{jk';t'}
\ee

\subsubsection{Lecture hall capacity}
In order to satisfy this constraint, the following inequality
need to hold.
\be{part-1}
\forall x \sum_{i}||x(i)| - \sum_{j}\sum_{k} S^{(X)}_{jk;x(i)}| \geq 0
\; .
\ee

\subsubsection{Student status and priority mapping}
\be{status-and-priority-mapping}
E_{assign} =
{1 \over 2}\sum_{m, \hat{x}, x}\sum_{i\neq i^{'}}
S^{(X)}_{\hat{x};x} 
\left( S^{(R)}_{r(i);(\hat{x}(m) \in x(m))}|({\it M = Req})
+
S^{(R)}_{r(i^{'});(\hat{x}(m) \in x(m))}|({\it M = \chi})\right)
\ee

\subsubsection{On the Energy Equation of the Potts Model}
The energy function of the Potts spin model is of the form:
\be{enrg}
E = E_{A} + E_{B} + E_{C} + \ldots
\ee

In (ref: Peterson-89), in which the TSP is used as study case,
parameters $\alpha$ and $\beta$ are used to control the strength of
E's restraints. The same parameters can also be used when dealing with
non-Euclidean type CO problems such as scheduling. 

\bo{ob-1}
\ \ For timetabling, $\beta$ is split into two parts ($\beta^{(Time)}$
and $\beta^{(X)}$), one with respect to the time-based neurons and the
other with respect to the space-based neurons. $E_{\beta}$ becomes as
follows: 
\be{diagonal-energy}
- {\beta^{(Time)} \over 2}\sum_j\sum_t {[S_{j;t}^{(T)}]}^2 - 
{\beta^{(X)} \over 2}\sum_{jk}\sum_{x} {[S_{jk;x}^{(X)}]}^2 
\ee
Here, the appearance of $\beta^{(Time)}$ and $\beta^{(X)}$ is only in
the diagonal terms. They play a crucial role in monitoring the
dynamics (phase transitions, etc.). These diagonal terms correspond to
self coupling interactions for the neurons of the network. In our 
simulations, each is set to $-0.1$.
\eo

\subsection{General Comments on the Algorithm}
Let ($N$) denotes the total number of the utilized neurons. At high 
$T$, the mean field solutions will tend to be states near to the
(symmetrical) maximum entropy state: $\forall (i, a) v_{ia} = {1/N}$.
Conversely at low $T$, finding a mean field solution will be
equivalent to using a local optimization method on the internal energy
--- a procedure highly sensitive to the initial conditions and known
to be ineffective (see p. 149 of (ref: Hop-85)).

In addition to the structure of the energy function, there are three
major interdependent issues which arise in completely specifying 
a mean field annealing algorithm for a timetabling problem:

\begin{itemize}
\item The values of the coefficients of terms in the energy function.
\item The types of dynamics used to find solutions of the mean field
 equations at each $T$.
\item The annealing schedule details, i.e. the initial temperature
 $T(0)$, the rules for deciding when to reduce $T$ and
 by how much, and the termination criteria.
\end{itemize}

A quantity introduced by (ref: Peterson-89) is called {\it
saturation}, $\Sigma$, defined as
\be{saturation}
\Sigma = {1 \over N} \sum_{ia} \nu^{2}_{ia} \;\; ,
\ee
to characterize the degree of clustering\footnote{Other measures could
also do this, e.g. an entropy measure $-\sum_{ia}\nu_{ia}ln(\nu_{ia})$,
but $\Sigma$ will be used here for consistency with the original
work.} of events in time and/or space; 
clearly $\Sigma_{min} = {1 \over N}$ and $\Sigma_{max} = 1$.

\subsection{The Algorithm}
\begin{itemize}
\item Choose a problem $\Longrightarrow$ (mapped into) $\{w_{ij}\}$.
\item Find the approximate phase transition temperature by linearizing 
equation (\ref{form-5.2}), see also above section.
\item Add a self-coupling $\beta$-term if necessary. Neural-wise, this
corresponds to a feedback connection from a neuron to itself.
\item Initialize the neurons $\nu_{ia}$ with ${1 \over N} \pm$ {\em
random values}; and set $T(0) = T_c$.
\item Until ($\Sigma \approx {1 \over N}\sum_{ia}\nu^{2}_{ia}\geq 0.99$) do:
\begin{itemize}
  \item At each $T(n)$, update all $\nu_{ia}$ via $\upsilon_{ia}$:
  \begin{itemize}
    \item  $\upsilon_{ia} = - {1 \over {T(n)}} \cdot {{\partial E[\nu]} \over {\partial \\item The types of dynamics used to find solutions of the mean field
nu_{ia}}}.$
    \item  $\nu_{ia} = {e^{\upsilon_{ia}} \over {\sum_b
e^{\upsilon_{ia}}}}. \;\;\; (\equiv f_{ia})$
  \end{itemize}
  \item $T(n+1) = 0.9 \times T(n)$.
\end{itemize}
\item After $\Sigma \geq 0.99$ is reached perform a {\it greedy
heuristics} if needed to account for possible imbalances or rule
violations.
\end{itemize}

\subsection{Some Implementation Details}
The general idea is to initialize the system close to the symmetric  
state at a temperature near $T_c$, then, allowing the system
to converge to a fixed point at each temperature, anneal (i.e. reduce 
$T$) until either the system reaches a stable high saturation
state or some other termination criterion is satisfied.
Specifically, the implementation used here has the following
ingredients\footnote{These follow the method used in
(ref: Peterson-89) but with some additions/improvements done by the
authors -see algorithm above.}: (Martin Simmen (ref: Simmen-93)
investigated in depth these and most of the other issues pertaining to 
the Potts model, using the TSP as a test case.)

\begin{itemize}
\item {\bf Initial State}\ \
Another possible initial value of $\nu^{n=0}_{ia}$ is $(1 + \xi_{ia})/N$,
with the $\xi_{ia}$ being random variables drawn uniformly from 
${[-\xi, \xi]}$. Clearly $\xi$ should be $\ll 1$.

\item {\bf Convergence}\ \
At each $T$ the $\nu_{ia}$ variables are updated as above or by
the following rule\footnote{Thanks to M. Simmen for pointing this out.
Also see (ref: Peterson-89)}:
\be{convergence-rule}
\nu^{(n+1)}_{ia} = \nu^{(n)}_{ia} + \gamma(f_{ia} - \nu^{(n)}_{ia}),
\;\; 0 < \gamma \leq 1\; .
\ee
Setting $\gamma = 1$ recovers the dynamics of (ref: Peterson-89),
and the $\gamma \rightarrow 0$ limit gives the first order Euler
integration method for the differential equation:
\be{diff-equation}
d\nu_{ia} / dt = f_{ia} - \nu_{ia}\;.
\ee
The fixed points (also referred to as the free energy minima)
of this equation constitute the mean field solutions.
The degree of convergence is monitored through some sort
of a ``tolerance'' quantity of the form 
$\sum_{ia}{(\nu_{ia} - f_{ia})}^2$.

\item {\bf Annealing Schedule} \ \
A simple exponential scheme is used, i.e. after reaching a fixed point
at $T$, $T \mapsto T \cdot T_r$, where 
$T_r \in (0,1)$ controls the cooling rate. Roughly speaking,
$T_r$ should be set to $0.9$.

\item {\bf Termination}\ \
One criterion is as outlined above in the algorithm. The other two
involve terminating runs which either fail to reach a fixed point
within a certain number of sweeps at some temperature, or which allow
$T$ to fall below some $T_{min}$, ($T_{min} \ll T_c$), without ever
satisfying the $\Sigma > \Sigma^{thresh}$ 
criterion. $T_{min}$ was set to $10^{-3}$.

\end{itemize}

\subsection{Some Performance Issues of MFA on the TTS}

Given $M_x$ = average number of rooms for each category of classes,
$N_x$ = the overall total number of rooms, and $N_t$ = total number of
slots , putting the figures of the students aside, with a computational
complexity of O($M_x N_x N_t$) per sweep, software implementation of
the MFA are unlikely to be competitive with conventional serial
algorithms such as 3-opt and Lin-Kernighan's, which empirically have
an overall O($N^2$) running time (ref: Johnson-90). Ideally, one would
like to study the likely performance of the method through numerical
simulations using a small step-size $\gamma$ to approximate continuous
time dynamics. However, this is computationally very expensive: not
only do the simulation times grow as $1 \over \gamma$, but also small
values of ($\gamma$) accentuate the difference between the expected
value $T_c^{exp}$ of the critical temperature $T_c$ and the
theoretically computed value $T_c^{theo}$. That is, low $\gamma$ runs
tend to depress $T_c^{exp}$ artificially and thereby degrade solution
quality, and as we observed, result in a number of constraint
violations. Many of these constraints are hard, but the majority are
classified to be soft. This later point can be compensated for by more
stringent annealing parameters, but at the cost of longer simulation
times. No attempts have been make to carry this out.
Also, we observe that the overall quality of the solution degrades as
the problem size increases. This degradation moves from soft
constraint violation to major violations of hard constraints, and
hence, illegal timetables.
From various correspondences with Andrew Gee, He postulated
the following as an observation:

\bo{MFA-problem}\ \
We observe that the NP-completeness of the problem (the timetabling,
that is) is manifested in the fact that the constraint polytope is
non-integral. Furthermore, for a decent number of constraints, there
will be many (probably in the thousands) time more non-integral
vertices than 0-1 vertices. [Note here, since there is no objective,
we just have to find any 0-1 point within the polytope to get a valid
timetable]. So, chances are that as we use some annealing mechanism to
derive the state vector towards the boundary of the polytope , it's
going to settle into a non-integral vertex, which is no use to anyone.
(Actually, Peterson alludes to these states in his paper with Gislen,
claiming that they are easily decomposed into two good solutions!)
Basically, there is no getting around this. Thinking of it from a
polyhedral perspective, we're trying to find the one integral vertex
in several thousands non-integral ones, and we're just driving the
state vector in a random direction towards the polytope boundary:
well, chances that it's not going to work.
\eo

One thing that remains unclear is to what extent very slow annealing can
bring MFA close to finding the optimal (not only legal) schedule.
Whereas simulated annealing, given certain conditions, guarantees
convergence to the optimal solutions for asymptotically slow
annealing, such results is not true for MFA algorithms, because MFA is
an approximation.

\subsection{Mapping Timetabling onto the Potts Neural Nets}

\subsubsection{The Relationship Between the Key Timetabling Entities}

Our problem consist of the following sets: teachers ($P$), classes
($C$), students ($S$), rooms ($R$), and time periods (${\cal I}$).
\vspace{0.1cm}
\begin{figure}[htb]
\begin{center}
\begin{picture}(230,230)
\psfig{figure=scheduling-net.eps,height=230pt,width=230pt}
\end{picture}
\end{center}
\caption{A one-way arrow indicates that the ``bottom'' entity depends
on the ``top'' one. A two-way arrow indicates a mutual dependency
between the two entities.}
\label{relations}
\end{figure}

\subsubsection{Neural Mapping}

Involve mapping the problem onto the Potts model.
\begin{itemize}
\item These multi-state neurons, forming the network, represent
different activities of the timetabling problem. 

\item Construct the constraints of the problem using these neurons.

\item Derive those MF equations controlling the dynamics of the net:

\begin{itemize}
\item Introduce the MF variables , for example, 

\be{rel.1}
V_{pq;xt} = <S_{pq;xt}>_T\;\;.
\ee

\item Start off with the partition function:
\be{part}
Z = {\displaystyle \sum_{[\vec{s}]}} e^{-E[\vec{s}] / T} \;.
\ee

\item Re-write $Z$ in terms of a multidimensional integral
over continuous vectors ($\vec{u}_i, \vec{v}_i$).

\item The MFT approximation to $<\vec{s}_i>$ is given by the value 
$\vec{u}_i$ at a saddle-point equations:\\

\be{mft.1}
{{\partial E} \over {\partial v_{ia}}} = 0
\ee
and
\be{mft.2}
{{\partial E} \over {\partial u_{ia}}} = 0\;,
\ee
generating a set of self-consistency equations in the MF variables:

\be{mfv.1}
u_{ia} = 
- {1 \over T} {{\partial E[\vec{v}]} \over {\partial v_{ia}}}\;,
\ee
and
\be{mfv.2}
v_{ia} = {e^{u_{ia}} \over {{\displaystyle \sum_b} e^{u_{ib}}}}\;.
\ee

Where $a$ denotes the vector components.

For events (pq, xt), we get the following MF variables:

\be{mfv.3}
u_{pq,xt} = 
- {1 \over T} {{\partial E} \over {\partial v_{pq,xt}}} \;,
\ee
and
\be{mfv.4}
v_{pq,xt} = 
        {e^{u_{pq,xt}} \over {{\displaystyle \sum_{x' t'}}
e^{u_{pq,x't'}}}} \;.
\ee

From the later follows, that ${\displaystyle \sum_{a}} v_{ia} =
1\;\;$. One can think of the mean field $v_{ia}$ as the probability
for the Potts neuron $i$ to be in state $a$.
\end{itemize}
\item Iterate $v_{pq, xt}$ equation with annealing, i.e. starting at a
high temperature, and successively lowering it in the course of the
process.  A phase transition is passed at $T = T_c$.

\item Estimate the critical temperature $T_c$, which sets the scale of
$T$, by expanding the above $v_{pq,xt}$ equation around the trivial
fix-point : at $T = 0$
\be{fixed}
v^{(0)}_{pq,xt} = {1 \over N_x N_t}\;.
\ee
\end{itemize}

\section{Overview of MFA}

One of the potential drawbacks of using simulated annealing for hard 
optimization problems is that finding a good solution can often take
an unacceptably long time. Mean-field annealing (MFA) attempts to
avoid this problem by using a deterministic approximation to simulated
annealing, by attempting to average over the statistics of the
annealing process.  
The result is improved execution speed at the expense of solution
quality. Although not strictly a continuous descent technique, MFA is
closely related to the Hopfield neural network~\cite{HKP,32}. 

Mean-field annealing has been successfully applied to high school
class scheduling~\cite{28}. For scheduling, it is advantageous to use
a Potts neural encoding to specify discrete neural variables (or
neurons) for the problem. This is defined in its simplest form as a
mapping of events onto space-time slots, for example an event $i$, in
this case a professor-class pair ($p,q$), is mapped onto a space-time
slot $a$, in this case a classroom-timeslot pair ($x,t$).  
Now, the Potts neurons $S_{ia}$ are defined to be $1$ if event $i$
takes place in space-time slot $a$, and $0$ otherwise.
In this way, the constraints involved can be embedded in the neural
net in terms of the weights $w_{i,j}$ of the neural network, which
encode a Potts normalization condition such as $\sum_{a} S_{ia} = 1$.

For a full derivation of the mean-field
annealing algorithm from its roots in statistical physics, see 
Hertz {\it et al.}~\cite{HKP} or  
Peterson {\it et al.}~\cite{26.1}. 

Here we will just give a brief overview of the method.
The basic idea is that it is possible to approximate the actual cost  
or energy function $E$, which is a function of discrete neural
variables $S_{ia}$, by an effective energy function $E^{'}$ that can
be represented in terms of continuous variables $U_{ia}$ and
$V_{ia}$. These are known as mean field variables, since $V_{ia}$ is
an approximation to the average value of $S_{ia}$ at a given
temperature $T$.

This approach effectively smooths out the energy function and makes it
easier to find the minimum value, which is obtained by solving 
the saddle point equations 
${\partial E^{'} \over \partial V_{ia}} = 0$ and  
${\partial E^{'} \over \partial U_{ia}} = 0$, 
which generate a set of self-consistent mean field theory (MFT)
equations in terms of the mean field variables $U$ and $V$:

\be{spe-1}
U_{ia} = - {1 \over T} {\partial E \over \partial V_{ia}}
\ee

\be{spe-2}
V_{ia} = {e^{U_{ia}} \over \sum_b e^{U_{ib}}}. 
\ee

The MFA algorithm involves solving equations \ref{spe-1} and
\ref{spe-2} at a series of progressively lower temperatures $T$: this
process is known as temperature annealing. 
The critical temperature $T_c$, which sets the scale of $T$, is
estimated by expanding equation \ref{spe-2} around the trivial
fixed-point~\cite{27, 28} $V^{(0)}_{ia} = {1 \over N_a}$, where $N_a$
is the number of possible states of each of the network neurons.
For example, for the events defined by professor-class pairs
$(p,q)$ mapped onto classroom-timeslots $(x,t)$, we have $N_p N_q$
neurons, each of which has $N_x N_t$ possible states, in which case
$V^{(0)}_{pq;xt} = {1 \over N_x N_t}$.

Equations \ref{spe-1} and \ref{spe-2} can be solved iteratively using 
either synchronous or serial updating. The iterative dynamics to
evolve the mean field variables toward a self-consistent solution 
is explained in detail by Peterson {\it et al.}~\cite{26}. 
The solutions correspond to stable states of the Hopfield
network~\cite{32}.
% so both techniques appear to be seeking the same solution points.
Observe from equation \ref{spe-2} that any solution to the MF
equations respects a continuous version of the Potts condition 

\be{potts-cond}
\sum_a V_{ia} = 1  \;\;\;\forall\;i .
\ee


\subsection{The Algorithm}

The generic MFA algorithm appears in Figure~\ref{MFA-algorithm}. 
At high temperatures $T$, the mean-field solutions will be states  
near the fixed-point symmetrical maximum entropy state $V_{ia} =
1/N_a$. 
At low temperatures, finding a mean-field solution will be equivalent
to using the Hopfield model, which is highly sensitive to the initial  
conditions and known to be ineffective for hard problems~\cite{32}.
MFA improves over the Hopfield model by using annealing to slowly 
decrease the temperature in order to sidestep these problems.
 
These characteristics are similar to those of simulated annealing,
which is no surprise since both it and the mean-field method compute 
thermal averages over Gibbs distributions of discrete states, the
former stochastically and the latter through a deterministic
approximation. It is therefore natural to couple the mean-field method
with the concept of annealing from high to low temperatures. 

In addition to the structure of the energy function, there are three
major interdependent issues which arise in completely specifying 
a mean-field annealing algorithm for a timetabling problem:

\begin{itemize}
\item The values of the coefficients of terms in the energy function.
\item The types of dynamics used to find solutions of the MFT
 equations at each $T$.
\item The annealing schedule details, i.e. the initial temperature
 $T(0)$, the rules for deciding when to reduce $T$ and
 by how much, and the termination criteria.
\end{itemize}

Peterson {\it et al.}~\cite{28} introduced a quantity 
called {\it saturation}, $\Sigma$, defined as
\be{saturation}
\Sigma = {1 \over N_i} \sum_{ia} V^{2}_{ia} \;\; ,
\ee
where $N_{i}$ is the number of events (in this case the number of 
professor-class pairs).
This characterizes the degree of clustering of events in time and/or
space, $\Sigma_{min} = {1 \over N_a}$ corresponds to high temperature,
whereas $\Sigma_{max} = 1$ means that all the $V_{ia}$ have converged
to 0 or 1 values, indicating that each event has been assigned to a
space-time slot.


\begin{figure}[tbp]
\begin{center}
\begin{enumerate}
\item Choose a problem and encode the constraints into weights
  $\{w_{ij}\}$. 
\item Find the approximate phase transition temperature by linearizing
  equation (\ref{spe-2}).
\item Add a self-coupling $\beta$-term if necessary. In a neural net, 
  this corresponds to a feedback connection from a neuron to itself.
\item Initialize the neurons $V_{ia}$ to high temperature values
  ${1 \over N_a}$ plus a small random term such as ${rand[-1,1] \times
    0.001}$; and set $T(0) = T_c$.
\item Until ($\Sigma \geq 0.99$) do:
\begin{itemize}
  \item At each $T(n)$, update all $U_{ia}$ and $V_{ia}$ by iterating
    to a solution of the mean field equations.
  \item $T(n+1) = \alpha T(n)$, we chose $\alpha = 0.9$
\end{itemize}
\item The discrete values $S_{ia}$ that specify the schedule are
  obtained by rounding the mean field values $V_{ia}$ to the nearest
  integer (0 or 1). 
\item Perform {\it greedy heuristics} if needed to account for possible 
   imbalances or rule violations.
\end{enumerate}
\end{center}
\caption{The Generic Mean-Field  Annealing Algorithm}
\label{MFA-algorithm}
\end{figure}

The first step of Figure~\ref{MFA-algorithm} is to map the constraints
of the problem into the neural net connection weights.
In our implementation, at each $T_n$ the MFA algorithm
(Figure~\ref{MFA-algorithm})  performs one update per neural variable 
(defined as one sweep) 
with sequential updating using equations \ref{spe-1} and \ref{spe-2}. 
After reaching a saturation value close to 1 (we chose $\Sigma = 0.99$)
we check whether the
obtained solutions are valid, i.e. $E_{hard} = 0$. If this is not the
case the network is re-initialized and is allowed to
resettle. We repeat this procedure a number of times until the best 
solution is found. 
A similar procedure was carried out on high school
scheduling by Peterson {\it et al.}~\cite{28}.

The MFA implementation was a little more complicated 
than the implementation of simulated annealing and the expert system,
since it had many more parameters to handle, and it was 
often more difficult to find optimal values for these parameters.  
For example, one complication is the computation of the critical
temperature $T_c$, which involved an iterative procedure of a
linearized dynamic system.  On the other hand, we observed that the 
convergence time was indeed much less than any of the convergence
times of the simulated annealing using the three annealing schedules
studied.
For more details on our MFA implementation, see Ref.~\cite{SalehMFA}.






















