#!/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