#!/usr/bin/python -Wall

# ================================================================
# John Kerl
# kerl at math dot arizona dot edu
# 2008-02-07
# ================================================================
# This program supports my paper
#   "The Metropolis-Hastings algorithm by example"
# which may be found at
#   http://math.arizona.edu/~kerl/doc/mhcoin.pdf
# ================================================================
# This software is released under the terms of the GNU GPL.
# Please see LICENSE.txt in the same directory as this file.
# ================================================================

from math import *
import sackmat_m
import random
import copy
import sys

from mhcoin_states  import *
from mhcoin_util    import *
from mhcoin_S       import *
from mhcoin_stochmx import *

# ----------------------------------------------------------------
# State energy for the independent (coin) model.
def compute_indep_energy(spins, Ps, beta):
	P = Ps[state_to_index(spins)]
	return -log(P) / beta

# ----------------------------------------------------------------
def compute_indep_energy_delta(spins, k, ps, beta):
	pk = ps[k]
	return spins[k] * log(pk/(1-pk)) / beta

# ----------------------------------------------------------------
def test_indep_energy_delta(n, ps, Ps, beta):
	states = all_states(n)
	num_states = len(states)

	print "Testing independent energy delta"
	print "ps = ",
	sackmat_m.print_row_vector(ps, "%11.7f")

	for i in range(0, num_states):
		spins = states[i]
		for k in range(0, n):

			short_way = compute_indep_energy_delta(spins, k, ps, beta)

			E  = compute_indep_energy(spins, Ps, beta)
			spins[k] = -spins[k]
			Ep = compute_indep_energy(spins, Ps, beta)
			spins[k] = -spins[k]
			long_way = Ep-E

			print ">> i=%d k=%d E=%11.7f E'=%11.7f E'-E = %11.7f DE = %11.7f diff=%f" % \
				(i, k, E,Ep,long_way, short_way, long_way-short_way)
	print

# ================================================================
# State energy for the Ising model.

def compute_Ising_energy(spins, S, h):
	n = len(spins)
	E = 0.0
	for i in range(0, n):
		for j in range(0, n):
			E -= S[i][j] * spins[i] * spins[j]
		E += h[i] * spins[i]
	return E

# ----------------------------------------------------------------
# See the paper for a derivation of this formula.
# Delta E =
#    2 sum_{j!=k} S_kj T_kT_j
#  + 2 sum_{i!=k} S_ik T_iT_k
#  - 2 h_k T_k.

def compute_Ising_energy_delta(spins, k, S, h):
	n = len(spins)

	DE = 0.0
	for j in range(0, n):
		if (j != k):
			DE += S[k][j] * spins[k] * spins[j]
	for i in range(0, n):
		if (i != k):
			DE += S[i][k] * spins[i] * spins[k]
	DE -= h[k] * spins[k]

	return 2*DE

# ----------------------------------------------------------------
def test_Ising_energy_delta(n, S, h):
	states = all_states(n)
	num_states = len(states)

	print "Testing Ising energy delta"

	for i in range(0, num_states):
		spins = states[i]
		for k in range(0, n):

			short_way = compute_Ising_energy_delta(spins, k, S, h)

			E  = compute_Ising_energy(spins, S, h)
			spins[k] = -spins[k]
			Ep = compute_Ising_energy(spins, S, h)
			spins[k] = -spins[k]
			long_way = Ep-E

			print ">> i=%d k=%d E=%11.7f E'=%11.7f E'-E = %11.7f DE = %11.7f diff=%f" % \
				(i, k, E,Ep,long_way, short_way, long_way-short_way)
	print

# ----------------------------------------------------------------
# Shannon-Gibbs entropy:  - sum_i P_i log2(P_i).

def compute_entropy(Ps):
	num_states = len(Ps)
	H = 0.0
	for i in range(0, num_states):
		P = Ps[i]
		H -= P * log(P) / log(2)
		#print "i = %2d P = %11.7f log2(P) = %11.7f" % \
		#	(i, P, log(P)/log(2))
	return H

# ================================================================
def get_IID_fair_Ps(n, verbose=0):
	if (verbose == 2):
		print "IID, fair:"

	ps = [0.5] * n
	Ps = []
	sum = 0.0
	for state in all_states(n):
		P = 1.0/2**n
		sum += P
		Ps.append(P)

		if (verbose == 2):
			print "  [",
			sackmat_m.print_row_vector_no_cr(state, "%2d")
			print "]",
			sackmat_m.print_row_vector_no_cr([P])
			print

	if (verbose == 2):
		print "sum P = %11.7f" % (sum)
		print

	return [ps, Ps]

# ----------------------------------------------------------------
def get_IID_p_Ps(p, n, verbose=0):
	if (verbose == 2):
		print "IID, p = %.1f:" % (p)

	ps = [p] * n
	Ps = []
	sum = 0.0
	for state in all_states(n):
		P = 1.0
		for ei in state:
			if (ei == +1): # +1, heads
				ni = 0
			else:         # -1, tails
				ni = 1
			P *= p**(1-ni) * (1-p)**ni
		sum += P
		Ps.append(P)

		if (verbose == 2):
			print "  [",
			sackmat_m.print_row_vector_no_cr(state, "%2d")
			print "]",
			sackmat_m.print_row_vector_no_cr([P])
			print
	if (verbose == 2):
		print "sum P = %11.7f" % (sum)
		print

	return [ps, Ps]

# ----------------------------------------------------------------
def get_indep_Ps(plo, phi, n, verbose=0):
	ps = lin_range(plo, phi, n)
	if (verbose == 2):
		print "IID, p's = [",
		sackmat_m.print_row_vector_no_cr(ps)
		print "]"

	Ps = []
	sum = 0.0
	for state in all_states(n):
		P = 1.0
		n = len(state)
		for i in range(0, n):
			ei = state[i]
			pi = ps[i]
			if (ei == +1): # +1, heads
				ni = 0
			else:         # -1, tails
				ni = 1
			P *= pi**(1-ni) * (1-pi)**ni
		sum += P
		Ps.append(P)

		if (verbose == 2):
			print "  [",
			sackmat_m.print_row_vector_no_cr(state, "%2d")
			print "]",
			sackmat_m.print_row_vector_no_cr([P])
			print
	if (verbose == 2):
		print "sum P = %11.7f" % (sum)
		print

	return [ps, Ps]

# ================================================================
def compute_Ising_Es_Z_and_Ps(states, S, h, beta, verbose=0):
	Z   = 0.0
	Eis = []
	Pis = []

	for i in range(0, num_states):
		Ei = compute_Ising_energy(states[i], S, h)
		if (verbose == 2):
			print "i = %3d E(" % (i),
			sackmat_m.print_row_vector_no_cr(states[i], "%2d")
			print ") = %11.7f   exp(-beta Ei) = %11.7f" % \
				(Ei, exp(-beta*Ei))
			#print "E(%2d) = %11.7f   exp(-beta Ei) = %11.7f" % \
			#	(i, Ei, exp(-beta*Ei))
		Z += exp(-beta * Ei)
		Eis.append(Ei)
	if (verbose == 2):
		print
		print "Z   = %11.7f" % (Z)
		print

	for i in range(0, num_states):
		Ei = Eis[i]
		Pi = exp(-beta*Ei) / Z
		Pis.append(Pi)

	return [Z, Eis, Pis]

# ================================================================
def get_indep_Markov_matrix(n, states, ps, beta):
	N = len(states)
	M = sackmat_m.make_zero_matrix(N, N)

	for old_state in states:
		i = state_to_index(old_state)
		sum_p = 0.0

		# Accessible states are those which differ by one spin, along with the
		# same state.
		for k in range(0, n):
			new_state = copy.copy(old_state)
			new_state[k] = -new_state[k]
			j = state_to_index(new_state)

			Delta_E = compute_indep_energy_delta(old_state, k, ps, beta)

			p = min(1.0, exp(-beta * Delta_E)) / n
			M[i][j] = p
			sum_p += p

		M[i][i] = 1.0 - sum_p
	return M

# ----------------------------------------------------------------
def get_Ising_Markov_matrix(n, states, S, h, Ising_Z):
	n = len(h)
	states = all_states(n)
	N = 2 ** n
	M = sackmat_m.make_zero_matrix(N, N)

	for old_state in states:
		i = state_to_index(old_state)
		sum_p = 0.0

		# Accessible states are those which differ by one spin, along with the
		# same state.
		for k in range(0, n):
			new_state = copy.copy(old_state)
			new_state[k] = -new_state[k]
			j = state_to_index(new_state)

			old_E = compute_Ising_energy(old_state, S, h)
			new_E = compute_Ising_energy(new_state, S, h)
			Delta_E = new_E - old_E

			p = min(1.0, exp(-beta * Delta_E)) / n
			M[i][j] = p
			sum_p += p

		M[i][i] = 1.0 - sum_p
	return M

# ================================================================
# CDF(k) = P(state index <= k)

def Ps_to_CDF(Ps):
	num_states = len(Ps)
	CDF = [0.0] * num_states
	sum = 0.0
	for k in range(0, num_states):
		sum += Ps[k]
		CDF[k] = sum
	return CDF

# ----------------------------------------------------------------
def print_Ps_and_CDF(states, Ps, CDF, name):
	print "Probabilities and CDF for %s:" % (name)
	num_states = len(Ps)
	sum = 0.0
	for i in range(0, num_states):
		Pi = Ps[i]
		state = states[i]
		print "  [",
		sackmat_m.print_row_vector_no_cr(state, "%2d")
		print "]",
		print "%11.7f" % (Pi),
		print "%11.7f" % (CDF[i]),
		print
		sum += Pi
	print "sum P = %11.7f" % (sum)
	H = compute_entropy(Ps)
	print "Shannon-Gibbs entropy = %11.7f volume = %11.7f" % (H, 2**H)
	print

# ----------------------------------------------------------------
def print_Es_and_Pstars(states, Ps, Es, Z, name):
	print "Energies and probability numerators for %s:" % (name)
	num_states = len(Ps)
	for i in range(0, num_states):
		state = states[i]
		print "  [",
		sackmat_m.print_row_vector_no_cr(state, "%2d")
		print "]",
		print "%11.7f" % (Es[i]),
		print "%11.7f" % (Ps[i]*Z), # P = P*/Z
		print
	print

# ================================================================
# Sample using independent coin flips.

def get_sample_Qs_indep(ps, n, nsamp, verbosity=0):
	num_states = 2**n
	counts = [0.0] * num_states
	state = [0] * n

	if (verbosity == 2):
		print "Getting %d samples using independence sampling" % (nsamp)
		print "ps = ",
		sackmat_m.print_row_vector(ps, "%11.7f")

	for sampi in range(0, nsamp):

		# xxx make an independent-random-state function.
		for i in range(0, n):
			# Example:  p = 0.99 of heads, which is +1.
			# So if the unit-pseudorandom value u is 0 <= u < 0.99
			# then we want +1.
			u = random.uniform(0.0, 1.0)
			if (u < ps[i]):
				state[i] = +1
			else:
				state[i] = -1

		if (verbosity == 2):
			print "  Sample %5d state = " % (sampi),
			sackmat_m.print_row_vector(state, "%2d")
		index = state_to_index(state)
		counts[index] += 1

	if (verbosity == 2):
		print

	# Scale the sample counts down to sample proportions.
	for k in range(0, num_states):
		counts[k] *= 1.0/nsamp
	return counts

# ----------------------------------------------------------------
# CDF sampling.
#
# CDF(k) = P(state # <= k)
# Example:  There are four states each with probability 0.25.
#
#  k   P(state == k)   P(state <= k)
# ---  -------------   -------------
#  0       0.25             0.25
#  1       0.25             0.50
#  2       0.25             0.75
#  3       0.25             1.00
#
# Get a pseudorandom number uniformly distributed on [0.0, 1.0).
# * if (u < CDF(0)): select state 0
# * if (u < CDF(1)): select state 1
# * if (u < CDF(2)): select state 2
# * if (u < CDF(3)): select state 3

def get_sample_Qs_from_CDF(CDF, n, num_states, nsamp, verbosity=0):
	num_states = len(CDF)
	counts = [0.0] * num_states

	if (verbosity == 2):
		print "Getting %d samples using CDF sampling" % (nsamp)

	for sampi in range(0, nsamp):
		u = random.uniform(0.0, 1.0)
		index = -1
		for k in range(0, num_states):
			if (u < CDF[k]):
				index = k
				break
		if (index == -1):
			print >> sys.stderr, \
				"Coding error detected in get_sample_Qs_from_CDF."
			sys.exit(1)
		if (verbosity == 2):
			print "  Sample %5d state = " % (sampi),
			sackmat_m.print_row_vector(index_to_state(index, n), "%2d")
		counts[index] += 1

	if (verbosity == 2):
		print

	# Scale the sample counts down to sample proportions.
	for k in range(0, num_states):
		counts[k] *= 1.0/nsamp
	return counts

# ----------------------------------------------------------------
def get_random_initial_spins(n, p_of_up):
	spins = [0] * n
	for k in range(0, n):
		u = random.uniform(0.0, 1.0)
		if (u < p_of_up):
			spins[k] =  1
		else:
			spins[k] = -1
	return spins

# ================================================================
def Metro_indep_step(spins, k, ps, Ps, beta, verbosity):
	N = len(spins)
	flipped = 0

	Delta_E = compute_indep_energy_delta(spins, k, ps, beta)

	p = min(1.0, exp(-beta * Delta_E))
	u = random.uniform(0.0, 1.0)
	if (u < p):
		spins[k] = -spins[k]
		flipped = 1

	if (verbosity >= 3):
		print "   Metro step at site %2d:" % (k),
		sackmat_m.print_row_vector_no_cr(spins, "%2d")
		print " ( %4d )" % (state_to_index(spins)),
		if (flipped):
			print " *",
		else:
			print " .",
		print "DE=%11.7f" % (Delta_E),
		print "E=%11.7f" % (compute_indep_energy(spins, Ps, beta))

# ----------------------------------------------------------------
def Metro_indep_sweep(spins, ps, Ps, beta, rand_metro_k, verbosity):
	n = len(spins)
	for j in range(0, n):
		if (rand_metro_k):
			k = random.randint(0,n-1)
		else:
			k = j
		Metro_indep_step(spins, k, ps, Ps, beta, verbosity)

# ----------------------------------------------------------------
# spins is the initial state.  It is modified.

def get_sample_Qs_from_Metro_indep(spins, ps, Ps, beta, num_sweeps, \
rand_metro_k, metro_verbosity, q_verbosity=0):

	n = len(spins)
	num_states = 2**n

	counts = [0] * num_states

	# Pass 1: thermalization
	# Pass 2: data accumulation
	descs = ["thermalization", "data accumulation"]
	for iter in [0, 1]:
		for sweepi in range(0, num_sweeps):
			#if (qv == 2): print spins
			if (metro_verbosity >= 2):
				print " Metro sweep %2d (%s):         " % (sweepi, descs[iter]),
				sackmat_m.print_row_vector_no_cr(spins, "%2d")
				print " ( %4d )" % (state_to_index(spins)),
				print "  ",
				print "%11.7f" % (compute_indep_energy(spins, Ps, beta))
			Metro_indep_sweep(spins, ps, Ps, beta, rand_metro_k, metro_verbosity)
			index = state_to_index(spins)
			if (iter == 1):
				counts[index] += 1
			if (metro_verbosity >= 3):
				print
		if (metro_verbosity >= 2):
			print

	# Scale the sample counts down to sample proportions.
	for k in range(0, num_states):
		counts[k] *= 1.0/num_sweeps
	return counts

# ================================================================
def Metro_Ising_step(spins, k, S, h, beta, verbosity):
	flipped = 0

	Delta_E = compute_Ising_energy_delta(spins, k, S, h)

	p = min(1.0, exp(-beta * Delta_E))
	u = random.uniform(0.0, 1.0)
	if (u < p):
		spins[k] = -spins[k]
		flipped = 1

	if (verbosity >= 3):
		print "  ",
		sackmat_m.print_row_vector_no_cr(spins, "%2d")
		print " site= %2d:" % (k),
		print "idx=  %4d " % (state_to_index(spins)),
		print "DE=%11.7f" % (Delta_E),
		print "U=%11.7f" % (u),
		print "V=%11.7f" % (p),
		if (flipped):
			print " acc",
		else:
			print " rej",
		print "E=%11.7f" % (compute_Ising_energy(spins, S, h))

# ----------------------------------------------------------------
def Metro_Ising_sweep(spins, S, h, beta, rand_metro_k, verbosity):
	n = len(spins)
	for j in range(0, n):
		if (rand_metro_k):
			k = random.randint(0,n-1)
		else:
			k = j
		Metro_Ising_step(spins, k, S, h, beta, verbosity)

# ----------------------------------------------------------------
# xxx examine this theoretically.
def Metro_Ising_wide_sweep(spins, S, h, beta, rand_metro_k, verbosity):
	n = len(spins)
	new_spins = unif_rand_state(n)

	old_E = compute_Ising_energy(    spins, S, h)
	new_E = compute_Ising_energy(new_spins, S, h)
	Delta_E = new_E - old_E

	flipped = 0
	p = min(1.0, exp(-beta * Delta_E))
	u = random.uniform(0.0, 1.0)
	if (u < p):
		#spins = new_spins
		for i in range(0, n):
			spins[i] = new_spins[i]
		flipped = 1

	if (verbosity >= 3):
		print "   Metro wide sweep:"
		sackmat_m.print_row_vector_no_cr(spins, "%2d")
		print " ( %4d )" % (state_to_index(spins)),
		if (flipped):
			print " *",
		else:
			print " .",
		print "DE=%11.7f" % (Delta_E),
		print "E=%11.7f" % (compute_Ising_energy(spins, S, h))

# ----------------------------------------------------------------
# spins is the initial state.  It is modified.

def get_sample_Qs_from_Metro_Ising(spins, S, h, beta, num_sweeps, \
rand_metro_k, metro_verbosity, q_verbosity=0):

	n = len(spins)
	num_states = 2**n
	counts = [0] * num_states

	# Pass 1: thermalization
	# Pass 2: data accumulation
	descs = ["thermalization", "data accumulation"]
	for iter in [0, 1]:
		for sweepi in range(0, num_sweeps):
			#if (qv == 2): print spins
			if (metro_verbosity >= 2):
				print " Metro sweep %2d (%s):         " % (sweepi, descs[iter]),
				sackmat_m.print_row_vector_no_cr(spins, "%2d")
				print " ( %4d )" % (state_to_index(spins)),
				print "  ",
				print "%11.7f" % (compute_Ising_energy(spins, S, h))

			#xxx
			Metro_Ising_sweep(spins, S, h, beta, rand_metro_k, metro_verbosity)
			#Metro_Ising_wide_sweep(spins, S, h, beta, rand_metro_k, metro_verbosity)

			index = state_to_index(spins)
			if (iter == 1):
				counts[index] += 1
			if (metro_verbosity >= 3):
				print
			if (metro_verbosity == 1):
				E = compute_Ising_energy(spins, S, h)
				print "%11.7f %11.7f # E" % (E, E*iter)

		if (metro_verbosity >= 2):
			print

	# Scale the sample counts down to sample proportions.
	for k in range(0, num_states):
		counts[k] *= 1.0/num_sweeps
	return counts

# ----------------------------------------------------------------
def raw_Metro_Ising(spins, S, h, beta, num_sweeps, \
rand_metro_k, metro_verbosity, q_verbosity=0):

	# Pass 1: thermalization
	# Pass 2: data accumulation
	descs = ["thermalization", "data accumulation"]
	counter = 0
	for iter in [0, 1]:
		for sweepi in range(0, num_sweeps):
			#if (qv == 2): print spins
			if (metro_verbosity >= 2):
				print " Metro sweep %2d (%s):         " % (sweepi, descs[iter]),
				sackmat_m.print_row_vector_no_cr(spins, "%2d")
				print " ( %4d )" % (state_to_index(spins)),
				print "  ",
				print "%11.7f" % (compute_Ising_energy(spins, S, h))

			Metro_Ising_sweep(spins, S, h, beta, rand_metro_k, metro_verbosity)

			if (metro_verbosity >= 3):
				print
			if (metro_verbosity == 1):
				E = compute_Ising_energy(spins, S, h)
				print "%6d %11.7f %11.7f # E" % (counter, E, E*iter)
			counter += 1

		if (metro_verbosity >= 2):
			print

# ================================================================
def report_Ps_and_Qs(states, Ps, Qs, nsamp, sample_style, name):
	num_states = len(states)
	print "P's and Q's: %d samples, %s sampling, %s:" \
		% (nsamp, sample_style, name)
	for k in range(0, num_states):
		state = states[k]
		print "  [",
		sackmat_m.print_row_vector_no_cr(state, "%2d")
		print "]",
		sackmat_m.print_row_vector_no_cr([Ps[k], Qs[k]])
		print "  diff =",
		sackmat_m.print_row_vector_no_cr([Ps[k] - Qs[k]])
		print
	print

# ================================================================
def usage():
	print >> sys.stderr, "Usage: %s [options ...]" % (sys.argv[0])
	print >> sys.stderr, "Options:"
	print >> sys.stderr, ""
	print >> sys.stderr, "  nis:         non-interacting"
	print >> sys.stderr, "  sis:         self-interacting"
	print >> sys.stderr, "  mfs:         mean-field"
	print >> sys.stderr, "  simfs:       self-interacting mean-field"
	print >> sys.stderr, "  ann:         aligning nearest-neighbor"
	print >> sys.stderr, "  aann:        anti-aligning nearest-neighbor"
	print >> sys.stderr, "  annp:        aligning nearest-neighbor periodic"
	print >> sys.stderr, "  aannp:       anti-aligning nearest-neighbor periodic"
	print >> sys.stderr, ""
	print >> sys.stderr, "  n=...:       number of spins"
	print >> sys.stderr, "  beta=...:    inverse temperature"
	print >> sys.stderr, "  p=...:       p for IID case"
	print >> sys.stderr, "  plo=...:     plo for independent case"
	print >> sys.stderr, "  phi=...:     phi for independent case"
	print >> sys.stderr, "  hl=...:      magnetic-field strength at left"
	print >> sys.stderr, "  hr=...:      magnetic-field strength at right"
	print >> sys.stderr, "  sc=...:      coupling strength"
	print >> sys.stderr, "  nsamp=...:   number of statistical samples"
	print >> sys.stderr, ""
	print >> sys.stderr, "  pup=[0.0-1.0]:   probability of up for Metro init state"
	print >> sys.stderr, "  nsweep=...:  number of Metropolis sweeps"
	print >> sys.stderr, "  rmk={0,1}:   randomize Metropolis sites"
	print >> sys.stderr, "               (else sweeps are sequential)"
	print >> sys.stderr, ""
	print >> sys.stderr, "  hv={0,1}:    print magnetic-field values"
	print >> sys.stderr, "  pv={0,1,2}:  print probabilities"
	print >> sys.stderr, "  qv={0,1}:    print sample proportions"
	print >> sys.stderr, "  sv={0,1}:    print connection matrix"
	print >> sys.stderr, "  zv={0,1}:    print Ei's and Z"
	print >> sys.stderr, "  mv={0-3}:    Metropolis verbosity"
	print >> sys.stderr, "  betav={0,1}: print beta"
	print >> sys.stderr, "  mmv={0,1}:   print Markov matrix and evolution of PMFs"
	print >> sys.stderr, "  mkp_reps=...: number of evolved PMFs to print"
	sys.exit(1)

# ================================================================
# ================================================================
# ================================================================
# Start of script

n     = 3
beta  = 1.0

zv    = 0
pv    = 1
qv    = 1
sv    = 0
hv    = 0
mv    = 0
betav = 0
mmv   = 0
mkp_reps = 10

p   = 0.9
plo = 0.6
phi = 0.8
pup = 1.0
rmk = 0

tde = 0

non_interacting_S                         = 1
self_interacting_S                        = 2
mean_field_S                              = 3
self_interacting_mean_field_S             = 4
aligning_nearest_neighbor_S               = 5
anti_aligning_nearest_neighbor_S          = 6
aligning_nearest_neighbor_periodic_S      = 7
anti_aligning_nearest_neighbor_periodic_S = 8

S_style = aligning_nearest_neighbor_S

sc = 1.0
hl = 0.0
hr = 0.0
nsamp = 10000
nsweep = 10000

metro_ising_only = 0

# ----------------------------------------------------------------
# Parse the command line.  Syntax: not -n 4 but n=4.
# This is crude (it's not syntax-error tolerant) but it's an easy-to-code
# hack.  It works in any language with an eval/exec (e.g. Python, Perl).
argc = len(sys.argv)
for argi in range(1, argc):
	if ((sys.argv[argi] == '-h') or (sys.argv[argi] == '--help')):
		usage()
	elif (sys.argv[argi] == "nis"):
		S_style = non_interacting_S
	elif (sys.argv[argi] == "sis"):
		S_style = self_interacting_S
	elif (sys.argv[argi] == "mfs"):
		S_style = mean_field_S
	elif (sys.argv[argi] == "simfs"):
		S_style = self_interacting_mean_field_S
	elif (sys.argv[argi] == "ann"):
		S_style = aligning_nearest_neighbor_S
	elif (sys.argv[argi] == "aann"):
		S_style = anti_aligning_nearest_neighbor_S
	elif (sys.argv[argi] == "annp"):
		S_style = aligning_nearest_neighbor_periodic_S
	elif (sys.argv[argi] == "aannp"):
		S_style = anti_aligning_nearest_neighbor_periodic_S
	else:
		v = -1
		exec(sys.argv[argi])

		if (v == 1):
			zv    = 1
			pv    = 1
			qv    = 1
			sv    = 1
			hv    = 1
			mv    = 1
			betav = 1

		if (v == 0):
			zv    = 0
			pv    = 0
			qv    = 0
			sv    = 0
			hv    = 0
			mv    = 0
			betav = 0
	pass

# ----------------------------------------------------------------
# Interaction matrix

if   S_style == non_interacting_S:
	S = get_non_interacting_S(n, sc)
elif S_style == self_interacting_S:
	S = get_self_interacting_S(n, sc)
elif S_style == mean_field_S:
	S = get_mean_field_S(n, sc)
elif S_style == self_interacting_mean_field_S:
	S = get_self_interacting_mean_field_S(n, sc)
elif S_style == aligning_nearest_neighbor_S:
	S = get_aligning_nearest_neighbor_S(n, sc)
elif S_style == anti_aligning_nearest_neighbor_S:
	S = get_anti_aligning_nearest_neighbor_S(n, sc)
elif S_style == aligning_nearest_neighbor_periodic_S:
	S = get_aligning_nearest_neighbor_periodic_S(n, sc)
elif S_style == anti_aligning_nearest_neighbor_periodic_S:
	S = get_anti_aligning_nearest_neighbor_periodic_S(n, sc)
else:
	print >> sys.stderr, "%s:  coding error:  S_style unrecognized." \
		% (sys.argv[0])
	sys.exit(1)

# Magnetic field

h = lin_range(hl,hr,n)

# Descriptions of the various ensembles, for display purposes.
d1 = "IID fair"
d2 = "IID p = %f" % (p)
#d3 = "indep. p = %f to %f" % (plo, phi)
d3 = "indep. w/ varying p"
d4 = "Ising model"

# ----------------------------------------------------------------
if metro_ising_only == 1:
	raw_Metro_Ising(get_random_initial_spins(n, pup), \
		S, h, beta, nsweep, rmk, mv, qv)
	sys.exit(1)

# ----------------------------------------------------------------
# Cases 1 through 3, which are easier due to independence.
[IID_fair_ps, IID_fair_Ps] = get_IID_fair_Ps(n, verbose=pv)
[IID_p_ps,    IID_p_Ps]    = get_IID_p_Ps(p, n, verbose=pv)
[indep_ps,    indep_Ps]    = get_indep_Ps(plo, phi, n, verbose=pv)

# ----------------------------------------------------------------
if (tde):
	test_indep_energy_delta(n, indep_ps, indep_Ps, beta)
	test_Ising_energy_delta(n, S, h)
	sys.exit(0)

# ----------------------------------------------------------------
# Ising 1D model.

if (betav):
	print "beta = %11.7f" % (beta)
	print
if (sv):
	print "S = "
	S.printf("%4.1f")
	print
if (hv):
	print "h = ",
	sackmat_m.print_row_vector(h)
	print

states = all_states(n)
num_states = len(states)

[Ising_Z, Ising_Es, Ising_Ps] = \
	compute_Ising_Es_Z_and_Ps(states, S, h, beta, verbose=zv)

# ----------------------------------------------------------------
if (mmv):

	print "Markov matrix for %s:" % (d1)
	M1 = get_indep_Markov_matrix(n, states, IID_fair_ps, beta)
	print_matrix_with_row_sums(M1)
	print

	print "Markov matrix for %s:" % (d2)
	M2 = get_indep_Markov_matrix(n, states, IID_p_ps, beta)
	print_matrix_with_row_sums(M2)
	print

	print "Markov matrix for %s:" % (d3)
	M3 = get_indep_Markov_matrix(n, states, indep_ps, beta)
	print_matrix_with_row_sums(M3)
	print

	print "Markov matrix for %s:" % (d4)
	M4 = get_Ising_Markov_matrix(n, states, S, h, Ising_Z); print
	print_matrix_with_row_sums(M4)
	print

	iP0 = sackmat_m.make_zero_matrix(1, num_states)
	iP0[0][0] = 1.0

	iP1 = sackmat_m.make_zero_matrix(1, num_states)
	iP1[0][-1] = 1.0

	iPU = sackmat_m.make_zero_matrix(1, num_states)
	for j in range(0, num_states):
		iPU[0][j] = 1.0/num_states

	mkp_Ms     = [ M1, M2, M3, M4]
	mkp_Mdescs = [ d1, d2, d3, d4 ]
	mkp_target_Ps = [IID_fair_Ps, IID_p_Ps, indep_Ps, Ising_Ps]
	mkp_iPs     = [ iP0, iP1, iPU ]
	mkp_iPdescs = [ "initial spins up", "initial spins down", "initial spins uniform random"]
	mkp_nMs = len(mkp_Ms)
	mkp_niPs = len(mkp_iPs)

	for i in range(0, mkp_nMs):

		print "================================================================"
		print "Target stable distribution: %s" % (mkp_Mdescs[i])
		print "----------------------------------------------------------------"
		print
		sackmat_m.print_row_vector(mkp_target_Ps[i])
		print

		print "================================================================"
		print "Markov-matrix iteration: %s" % (mkp_Mdescs[i])
		print "----------------------------------------------------------------"
		# Print powers of the Markov matrices themselves.
		M  = copy.deepcopy(mkp_Ms[i])
		Mk = copy.deepcopy(M)
		print
		for k in range(0, mkp_reps):
			print "k =", k
			Mk.printf(); print
			Mk  = Mk * M
		print

		# Print evolution of the PMFs under right multiplication by the Markov
		# matrices.
		for j in range(0, mkp_niPs):
			print "================================================================"
			print "PMF iteration: %s, %s" % (mkp_Mdescs[i], mkp_iPdescs[j])
			print "----------------------------------------------------------------"
			M  = copy.deepcopy(mkp_Ms[i])
			Pk = copy.deepcopy(mkp_iPs[j])

			print
			for k in range(0, mkp_reps):
				Pk.printf();
				Pk = Pk * M
			print

print "================================================================"

# ----------------------------------------------------------------
IID_fair_CDF = Ps_to_CDF(IID_fair_Ps)
if (pv):
	print_Ps_and_CDF(states, IID_fair_Ps, IID_fair_CDF, d1)

IID_p_CDF    = Ps_to_CDF(IID_p_Ps)
if (pv):
	print_Ps_and_CDF(states, IID_p_Ps,    IID_p_CDF,    d2)

indep_CDF    = Ps_to_CDF(indep_Ps)
if (pv):
	print_Ps_and_CDF(states, indep_Ps,    indep_CDF,    d3)

Ising_CDF    = Ps_to_CDF(Ising_Ps)
if (pv):
	print_Ps_and_CDF(states, Ising_Ps,    Ising_CDF,    d4)
	print_Es_and_Pstars(states, Ising_Ps, Ising_Es, Ising_Z, d4)
	print "----------------------------------------------------------------"
	print

# ================================================================
# Sample from the first three distributions using independent coin tosses.

IID_fair_indep_Qs = get_sample_Qs_indep(IID_fair_ps, n, nsamp, qv)
if (qv):
	report_Ps_and_Qs(states, IID_fair_Ps, IID_fair_indep_Qs, nsamp, "indep.", d1)

IID_p_indep_Qs    = get_sample_Qs_indep(IID_p_ps,    n, nsamp, qv)
if (qv):
	report_Ps_and_Qs(states, IID_p_Ps,    IID_p_indep_Qs,    nsamp, "indep.", d2)

indep_indep_Qs    = get_sample_Qs_indep(indep_ps,    n, nsamp, qv)
if (qv):
	report_Ps_and_Qs(states, indep_Ps,    indep_indep_Qs,    nsamp, "indep.", d3)
	print
	print "----------------------------------------------------------------"
	print

# ----------------------------------------------------------------
# Sample from all four distributions using CDF sampling.

IID_fair_CDF_Qs = get_sample_Qs_from_CDF(IID_fair_CDF, n, num_states, nsamp, qv)
if (qv):
	report_Ps_and_Qs(states, IID_fair_Ps, IID_fair_CDF_Qs, nsamp, "CDF", d1)

IID_p_CDF_Qs    = get_sample_Qs_from_CDF(IID_p_CDF,    n, num_states, nsamp, qv)
if (qv):
	report_Ps_and_Qs(states, IID_p_Ps,    IID_p_CDF_Qs,    nsamp, "CDF", d2)

indep_CDF_Qs    = get_sample_Qs_from_CDF(indep_CDF,    n, num_states, nsamp, qv)
if (qv):
	report_Ps_and_Qs(states, indep_Ps,    indep_CDF_Qs,    nsamp, "CDF", d3)

Ising_CDF_Qs    = get_sample_Qs_from_CDF(Ising_CDF,    n, num_states, nsamp, qv)
if (qv):
	report_Ps_and_Qs(states, Ising_Ps,    Ising_CDF_Qs,    nsamp, "CDF", d4)
	print
	print "----------------------------------------------------------------"
	print

# ----------------------------------------------------------------
# Sample from all four distributions using Metropolis sampling.

IID_fair_Metro_Qs = get_sample_Qs_from_Metro_indep(get_random_initial_spins(n, pup), \
	IID_fair_ps, IID_fair_Ps, beta, nsweep, rmk, mv, qv)
if (qv):
	report_Ps_and_Qs(states, IID_fair_Ps, IID_fair_Metro_Qs, nsamp, "Metro", d1)

IID_p_Metro_Qs    = get_sample_Qs_from_Metro_indep(get_random_initial_spins(n, pup), \
	IID_p_ps, IID_p_Ps, beta, nsweep, rmk, mv, qv)
if (qv):
	report_Ps_and_Qs(states, IID_p_Ps,    IID_p_Metro_Qs,    nsamp, "Metro", d2)

indep_Metro_Qs    = get_sample_Qs_from_Metro_indep(get_random_initial_spins(n, pup), \
	indep_ps, indep_Ps, beta, nsweep, rmk, mv, qv)
if (qv):
	report_Ps_and_Qs(states, indep_Ps,    indep_Metro_Qs,    nsamp, "Metro", d3)

Ising_Metro_Qs = get_sample_Qs_from_Metro_Ising(get_random_initial_spins(n, pup), \
	S, h, beta, nsweep, rmk, mv, qv)
if (qv):
	report_Ps_and_Qs(states, Ising_Ps,    Ising_Metro_Qs,    nsweep, "Metro", d4)
	print
	print "----------------------------------------------------------------"
	print