diff --git a/smcRowCol.py b/smcRowCol.py new file mode 100755 index 0000000..e9a4062 --- /dev/null +++ b/smcRowCol.py @@ -0,0 +1,135 @@ +"""smc.py + +2D constrained channel using fully adapted SMC. + +""" +import argparse +import numpy as np +from scipy import misc +from numpy.random import random_sample +import time + +def main(): + # Parse command-line arguments + parser = argparse.ArgumentParser(usage=__doc__) + parser.add_argument("--nx", type=int, help="number of rows") + parser.add_argument("--ny", type=int, help="number of columns") + parser.add_argument("--iter", type=int, help="iterations to run") + parser.add_argument("--particles", type=str, help="list with particles N") + args = parser.parse_args() + + # Number of particles - N (design variable) + particles = args.particles.split(',') + iterations = args.iter + CM = np.zeros( iterations ) + + # Initialize arrays/matrices + nx = args.nx + ny = args.ny + xDomain = np.arange( 2 ) + interactionPotentials = np.ones( (2, 2) ) + interactionPotentials[1, 1] = 0 + + fileName = str(nx)+'x'+str(ny)+'informationRateTEST.txt' + + # New file, print initial info, first line + f = open(fileName, 'w') + f.write('nx ny\n') + f.write(str(nx) + ' ' + str(ny) + '\n') + f.write('particles informationRateEst \n') + f.close() + + for iIter in particles: + print iIter + for i in np.arange(0,iterations): + N = int(iIter) + print N + # SMC initializations + trajectorySMC = np.zeros( (N, nx, ny) ) + tempWeights = np.ones( N ) + ancestors = np.zeros( N, np.int ) + tempParticleDist = np.zeros( len(xDomain) ) + + # BP initializations + messages = np.zeros( (N, nx, len(xDomain)) ) + + # --------------- + # SMC + # --------------- + # SMC first iteration, forward filtering + for iForward in np.arange(0,nx-1): + for iParticle in np.arange(0,N): + if iForward > 0 and iForward < nx-1: + messages[iParticle,iForward,:] = np.dot(interactionPotentials, messages[iParticle,iForward-1,:]) + else: + messages[iParticle,iForward,:] = np.dot( interactionPotentials, np.ones(2) ) + tempWeights = np.sum( messages[:,nx-2,:], axis=1 ) + CM[i] += np.log2(np.mean( tempWeights )) + + for iBackward in np.arange(0,nx)[::-1]: + for iParticle in np.arange(0,N): + tempParticleDist = np.ones( len(xDomain) ) + if iBackward < nx-1: + tempParticleDist = tempParticleDist * interactionPotentials[trajectorySMC[iParticle, iBackward+1, 0],:] + if iBackward > 0: + tempParticleDist = tempParticleDist * messages[iParticle, iBackward-1,:] + tempParticleDist = tempParticleDist / np.sum( tempParticleDist ) + ind = discreteSampling( tempParticleDist, xDomain, 1) + trajectorySMC[iParticle, iBackward, 0] = xDomain[ind] + + # SMC MAIN LOOP + for iSMC in np.arange(1,ny): + # Forward filtering + for iForward in np.arange(0,nx): + for iParticle in np.arange(0,N): + if iForward > 0 and iForward < nx-1: + messages[iParticle,iForward,:] = np.dot(interactionPotentials, interactionPotentials[trajectorySMC[iParticle, iForward, iSMC-1],:] * messages[iParticle,iForward-1,:]) + elif iForward == 0: + messages[iParticle,iForward,:] = np.dot(interactionPotentials, interactionPotentials[trajectorySMC[iParticle, iForward, iSMC-1],:]) + else: + messages[iParticle,iForward,:] = interactionPotentials[trajectorySMC[iParticle, iForward, iSMC-1],:] * messages[iParticle,iForward-1,:] + + tempWeights = np.sum( messages[:,nx-1,:], axis=1) + CM[i] += np.log2( np.mean( tempWeights ) ) + + if np.sum( tempWeights ) == 0: + raw_input() + # Sample ancestors + ancestors = discreteSampling( tempWeights / np.sum( tempWeights ), np.arange(N), N) + + # Backward sampling + for iBackward in np.arange(0,nx)[::-1]: + for iParticle in np.arange(0,N): + tempParticleDist = interactionPotentials[trajectorySMC[ancestors[iParticle], iBackward, iSMC-1],:] + if iBackward < nx-1: + tempParticleDist = tempParticleDist * interactionPotentials[trajectorySMC[iParticle, iBackward+1, iSMC],:] + if iBackward > 0: + tempParticleDist = tempParticleDist * messages[ancestors[iParticle],iBackward-1,:] + tempParticleDist = tempParticleDist / np.sum(tempParticleDist) + ind = discreteSampling( tempParticleDist, xDomain, 1 ) + trajectorySMC[iParticle, iBackward, iSMC] = xDomain[ind] + + trajectorySMC[:, :, :iSMC] = trajectorySMC[ancestors.astype(int), :, :iSMC].reshape( (N, nx, iSMC) ) + + CM[i] = CM[i] / (nx*ny) + f = open(fileName, 'a') + f.write(str(iIter) + ' ') + np.savetxt(f, CM.reshape( (1,iterations) )) + f.close() + #print 'MI est:',CM + +def discreteSampling(weights, domain, nrSamples): + bins = np.cumsum(weights) + return domain[np.digitize(random_sample(nrSamples), bins)] + +def ravel_multi_index(coord, shape): + return coord[0] * shape[1] + coord[1] + +def unravel_index(coord, shape): + iy = np.remainder(coord, shape[1]) + ix = (coord - iy) / shape[1] + return ix, iy + +if __name__ == "__main__": + main() + diff --git a/smcRowStrips.py b/smcRowStrips.py new file mode 100755 index 0000000..9add304 --- /dev/null +++ b/smcRowStrips.py @@ -0,0 +1,187 @@ +"""smc.py + +2D constrained channel using fully adapted SMC adding several columns at a time. + +""" +from optparse import OptionParser +import numpy as np +from scipy import misc +from numpy.random import random_sample +import time + +def main(): + # Parse command-line arguments + parser = OptionParser() + parser.add_option("--nx", type=int, help="number of rows") + parser.add_option("--ny", type=int, help="number of columns") + parser.add_option("--iter", type=int, help="iterations to run") + parser.add_option("--particles", type=str, help="list with particles N") + parser.add_option("--stripWidth", type=int, help="width of strip") + (args, options) = parser.parse_args() + + # Number of particles - N (design variable) + particles = args.particles.split(',') + iterations = args.iter + CM = np.zeros( iterations ) + + # Initialize arrays/matrices + nx = args.nx + ny = args.ny + stripWidth = args.stripWidth + nrInd = int(ny/stripWidth) + # Interaction potentials + interactionPotentials = np.ones( (2, 2) ) + interactionPotentials[1, 1] = 0 + intPot = np.ones( (2, 2) ) + intPot[1, 1] = 0 + nrComb = 2**stripWidth + xDomain = np.zeros((nrComb,stripWidth), dtype=np.uint64) + for iIter in range(nrComb): + tempStr = np.binary_repr(iIter, width=stripWidth) + for iStrip in range(stripWidth): + xDomain[iIter, iStrip] = int(tempStr[iStrip]) + fileName = str(nx)+'x'+str(ny)+'informationRateROWstripW' + str(stripWidth) + '.txt' + + # New file, print initial info, first line + f = open(fileName, 'w') + f.write('nx ny\n') + f.write(str(nx) + ' ' + str(ny) + '\n') + f.write('particles timeElapsed informationRateEst \n') + f.close() + + for iIter in particles: + print iIter + startSMC = time.time() + CM = np.zeros( iterations ) + for i in np.arange(0,iterations): + N = int(iIter) + # SMC initializations + trajectorySMC = np.zeros( (N, nx, ny) ) + tempWeights = np.ones( N ) + ancestors = np.zeros( N, np.int ) + + # BP initializations + messages = np.zeros( (N, nrComb, nx-1) ) + normConstMessages = np.zeros( (N, nx) ) + + # --------------- + # SMC + # --------------- + # CSMC first iteration, forward filtering + for iParticle in range(N): + for iRow in range(nx-1): + for iCur in range(nrComb): + tempDist = np.ones(nrComb) + for iPrev in range(nrComb): + for iStrip in range(stripWidth): + if iStrip != stripWidth-1: + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iPrev, iStrip+1]] + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iCur, iStrip]] + if iRow > 0: + tempDist[iPrev] *= messages[iParticle,iPrev, iRow-1] + messages[iParticle,iCur,iRow] = np.sum(tempDist) + normConstMessages[iParticle,iRow] = np.sum(messages[iParticle,:,iRow]) + messages[iParticle,:,iRow] = messages[iParticle,:,iRow] / normConstMessages[iParticle,iRow] + # Column sum + tempDist = np.ones(nrComb) + for iCur in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + tempDist[iCur] *= messages[iParticle, iCur, nx-2] + normConstMessages[iParticle,-1] = np.sum( tempDist ) + tempWeights = np.prod(normConstMessages, axis=1) + CM[i] += np.log2(np.mean( tempWeights )) + + # First iteration, Backward sampling + for iParticle in range(N): + for iRow in range(nx)[::-1]: + tempDist = np.ones( nrComb ) + for iCur in range(nrComb): + for iStrip in range(stripWidth): + if iStrip != stripWidth-1: + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + if iRow < nx-1: + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], trajectorySMC[iParticle,iRow+1, iStrip]] + if iRow > 0: + tempDist[iCur] *= messages[iParticle, iCur, iRow-1] + tempDist = tempDist / np.sum(tempDist) + curInd = discreteSampling( tempDist, range(nrComb), 1 ) + for iStrip in range(stripWidth): + trajectorySMC[iParticle,iRow,iStrip] = xDomain[curInd,iStrip] + + # SMC MAIN LOOP + for iSMC in np.arange(1,nrInd): + # BP initializations + messages = np.zeros( (N, nrComb, nx-1) ) + normConstMessages = np.zeros( (N, nx) ) + # Forward filtering + for iParticle in range(N): + for iRow in range(nx-1): + for iCur in range(nrComb): + tempDist = np.ones(nrComb) + for iPrev in range(nrComb): + for iStrip in range(stripWidth): + if iStrip != stripWidth-1: + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iPrev, iStrip+1]] + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iCur, iStrip]] + tempDist[iPrev] *= intPot[xDomain[iPrev,0], trajectorySMC[iParticle, iRow, iSMC*stripWidth-1]] + if iRow > 0: + tempDist[iPrev] *= messages[iParticle,iPrev, iRow-1] + messages[iParticle,iCur,iRow] = np.sum(tempDist) + normConstMessages[iParticle,iRow] = np.sum(messages[iParticle,:,iRow]) + messages[iParticle,:,iRow] = messages[iParticle,:,iRow] / normConstMessages[iParticle,iRow] + # Column sum + tempDist = np.ones(nrComb) + for iCur in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + tempDist[iCur] *= intPot[xDomain[iCur,0], trajectorySMC[iParticle, nx-1, iSMC*stripWidth-1]] + tempDist[iCur] *= messages[iParticle, iCur, nx-2] + normConstMessages[iParticle,nx-1] = np.sum( tempDist ) + tempWeights = np.prod(normConstMessages, axis=1) + CM[i] += np.log2(np.mean( tempWeights )) + + # Sample ancestors + ancestors = discreteSampling( tempWeights / np.sum( tempWeights ), np.arange(N), N) + + # Backward sampling + for iParticle in range(N): + for iRow in range(nx)[::-1]: + tempDist = np.ones( nrComb ) + for iCur in range(nrComb): + for iStrip in range(stripWidth): + if iStrip != stripWidth-1: + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + if iRow < nx-1: + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], trajectorySMC[iParticle,iRow+1, iSMC*stripWidth+iStrip]] + tempDist[iCur] *= intPot[xDomain[iCur,0], trajectorySMC[ancestors[iParticle], iRow, iSMC*stripWidth-1]] + if iRow > 0: + tempDist[iCur] *= messages[ancestors[iParticle], iCur, iRow-1] + tempDist = tempDist / np.sum(tempDist) + # Sampling + curInd = discreteSampling( tempDist, range(nrComb), 1 ) + for iStrip in range(stripWidth): + trajectorySMC[iParticle,iRow,iSMC*stripWidth+iStrip] = xDomain[curInd,iStrip] + trajectorySMC[:, :, :iSMC*stripWidth] = trajectorySMC[ancestors.astype(int), :, :iSMC*stripWidth] + CM[i] = CM[i] / (nx*ny) + f = open(fileName, 'a') + f.write(str(iIter) + ' ' + str(time.time() - startSMC) + ' ') + np.savetxt(f, CM.reshape( (1,iterations) )) + f.close() + #print 'MI est:',CM + +def discreteSampling(weights, domain, nrSamples): + bins = np.cumsum(weights) + return domain[np.digitize(random_sample(nrSamples), bins)] + +def ravel_multi_index(coord, shape): + return coord[0] * shape[1] + coord[1] + +def unravel_index(coord, shape): + iy = np.remainder(coord, shape[1]) + ix = (coord - iy) / shape[1] + return ix, iy + +if __name__ == "__main__": + main() + diff --git a/treeSampler.py b/treeSampler.py new file mode 100755 index 0000000..5dd3349 --- /dev/null +++ b/treeSampler.py @@ -0,0 +1,145 @@ +"""treeSampler.py + +Sample from 2D constrained channel using tree sampling. + +""" +from optparse import OptionParser +import numpy as np +from scipy import misc +from numpy.random import random_sample +import time + +def main(): + # Parse command-line arguments + parser = OptionParser() + parser.add_option("--nx", type=int, help="number of rows") + parser.add_option("--ny", type=int, help="number of columns") + parser.add_option("--iter", type=int, help="iterations to run") + parser.add_option("--chains", type=int, help="number of independent chains") + (args, options) = parser.parse_args() + + # Initialize arrays/matrices + nx = args.nx + ny = args.ny + nrChains = args.chains + + iterNr = args.iter + output = np.zeros( (nrChains, nx, ny) ) + #output[0,:,:] = (np.random.randint(2,size=nx*ny)).reshape((nx,ny)) + fileName = str(nx) + 'x' + str(ny) + 'TSTIME.txt' + # Block structure + temp = np.arange(ny) + block1 = temp[::2] + block2 = np.setdiff1d(temp, block1) + + # Interaction potentials + intPot = np.ones( (2,2), dtype=np.uint64 ) + intPot[1,1] = 0 + xDomain = np.arange(2) + + f = open(fileName, 'w') + f.write('nx ny\n') + f.write(str(nx) + ' ' + str(ny) + '\n') + f.write('iterNr timeElapsed zA zB\n') + f.close() + + # Main loop + for iIter in np.arange(iterNr): + print iIter + startTS = time.time() + fMarginalA = np.zeros( (nrChains, ny), dtype=np.uint64 ) + fMarginalB = np.zeros( (nrChains, ny), dtype=np.uint64 ) + for iChain in range(nrChains): + # -------------- + # First tree + # -------------- + #print output + for iCol in block1: + messages = np.ones( (len(xDomain), nx), dtype=np.uint64 ) + # Forward filtering + for iRow in range(nx): + tempDist = np.ones( 2, dtype=np.uint64 ) + if iCol > 0: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol-1]] + if iCol < ny-1: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol+1]] + if iRow > 0 and iRow < nx-1: + messages[:,iRow] = np.dot(intPot, tempDist * messages[:,iRow-1]) + elif iRow == 0: + messages[:,iRow] = np.dot(intPot, tempDist) + else: + messages[:,iRow] = tempDist * messages[:,iRow-1] + fMarginalA[iChain,iCol] = np.sum( messages[:,iRow] ) + + # Backward sampling + for iRow in range(nx)[::-1]: + tempDist = np.ones( 2 ) + if iCol > 0: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol-1]] + if iCol < ny-1: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol+1]] + if iRow > 0: + tempDist = tempDist * messages[:,iRow-1] + if iRow < nx-1: + tempDist = tempDist * intPot[:,output[iChain,iRow+1, iCol]] + tempDist = tempDist / np.sum(tempDist) + output[iChain,iRow, iCol] = discreteSampling( tempDist, xDomain, 1 ) + #print output + # -------------- + # Second block + # ------------- + for iCol in block2: + messages = np.ones( (len(xDomain), nx), dtype=np.uint64 ) + # Forward filtering + for iRow in range(nx): + tempDist = np.ones( 2 ) + if iCol > 0: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol-1]] + if iCol < ny-1: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol+1]] + if iRow > 0 and iRow < nx-1: + messages[:,iRow] = np.dot(intPot, tempDist * messages[:,iRow-1]) + elif iRow == 0: + messages[:,iRow] = np.dot(intPot, tempDist) + else: + messages[:,iRow] = tempDist * messages[:,iRow-1] + fMarginalB[iChain,iCol] = np.sum( messages[:,iRow] ) + + # Backward sampling + for iRow in range(nx)[::-1]: + tempDist = np.ones( 2 ) + if iCol > 0: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol-1]] + if iCol < ny-1: + tempDist = tempDist * intPot[:, output[iChain,iRow,iCol+1]] + if iRow > 0: + tempDist = tempDist * messages[:,iRow-1] + if iRow < nx-1: + tempDist = tempDist * intPot[:,output[iChain, iRow+1, iCol]] + tempDist = tempDist / np.sum(tempDist) + output[iChain, iRow, iCol] = discreteSampling( tempDist, xDomain, 1 ) + #print fMarginal + #raw_input() + fMarginal = np.zeros( (nrChains,ny) , dtype=np.uint64) + fMarginal[:,:len(block1)] = fMarginalA[:,block1] + fMarginal[:,len(block1):len(block1)+len(block2)] = fMarginalB[:,block2] + f = open(fileName, 'a') + f.write(str(iIter) + ' ' + str(time.time()-startTS) + ' ') + np.savetxt(f, fMarginal.reshape( (1,nrChains*ny) ), fmt='%u') + f.close() + +def discreteSampling(weights, domain, nrSamples): + bins = np.cumsum(weights) + return domain[np.digitize(random_sample(nrSamples), bins)] + +def ravel_multi_index(coord, shape): + return coord[0] * shape[1] + coord[1] + +def unravel_index(coord, shape): + iy = np.remainder(coord, shape[1]) + ix = (coord - iy) / shape[1] + return ix, iy + + +if __name__ == "__main__": + main() diff --git a/treeSamplerStrips.py b/treeSamplerStrips.py new file mode 100755 index 0000000..e1ce600 --- /dev/null +++ b/treeSamplerStrips.py @@ -0,0 +1,221 @@ +"""treeSampler.py + +Sample from 2D constrained channel using tree sampling + larger blocks/strips. + +""" +from optparse import OptionParser +import numpy as np +from scipy import misc +from numpy.random import random_sample +import time + +def main(): + # Parse command-line arguments + parser = OptionParser() + parser.add_option("--nx", type=int, help="number of rows") + parser.add_option("--ny", type=int, help="number of columns") + parser.add_option("--iter", type=int, help="iterations to run") + parser.add_option("--chains", type=int, help="number of independent chains") + parser.add_option("--stripWidth", type=int, help="width of the strips") + (args, options) = parser.parse_args() + + # Initialize arrays/matrices + nx = args.nx + ny = args.ny + nrChains = args.chains + stripWidth = args.stripWidth + iterNr = args.iter + output = np.zeros( (nrChains, nx, ny) ) + fileName = str(nx) + 'x' + str(ny) + 'TSstripw' + str(stripWidth) + '.txt' + # Block structure + index = np.arange(ny) + nrInd1 = int(np.ceil(float(ny)/float(stripWidth)/2)) + nrInd2 = ny/stripWidth - nrInd1 + #print nrInd1, nrInd2 + block1 = np.zeros( (nrInd1,stripWidth), dtype=np.uint64 ) + block2 = np.zeros( (nrInd2,stripWidth), dtype=np.uint64 ) + temp = 0 + for iIndex in range(nrInd1): + block1[iIndex,:] = index[stripWidth*iIndex+temp:stripWidth*(iIndex+1)+temp] + temp = temp + stripWidth + temp = stripWidth + for iIndex in range(nrInd2): + block2[iIndex,:] = index[stripWidth*iIndex+temp:stripWidth*(iIndex+1)+temp] + temp = temp + stripWidth + #print block1, block2 + + # Interaction potentials + intPot = np.ones( (2,2), dtype=np.uint64 ) + intPot[1,1] = 0 + nrComb = 2**stripWidth + xDomain = np.zeros((nrComb,stripWidth), dtype=np.uint64) + for iIter in range(nrComb): + tempStr = np.binary_repr(iIter, width=stripWidth) + for iStrip in range(stripWidth): + xDomain[iIter, iStrip] = int(tempStr[iStrip]) + #nrDimArray = '(' + '2,'*stripWidth + str(nx) + ')' + #print xDomain + + # Create file + f = open(fileName, 'w') + f.write('nx ny\n') + f.write(str(nx) + ' ' + str(ny) + '\n') + f.write('iterNr timeElapsed zA zB\n') + f.close() + + # Main loop + for iIter in np.arange(iterNr): + print iIter + startTS = time.time() + fMarginalA = np.zeros( (nrChains, nrInd1) ) + fMarginalB = np.zeros( (nrChains, nrInd2) ) + for iChain in range(nrChains): + # -------------- + # First tree + # -------------- + #print output[iChain,:,:].reshape( (nx,ny) ) + for iBlock in range(nrInd1): + messages = np.ones( (nrComb, nx-1) ) + normConstMessages = np.zeros( nx ) + # Forward filtering + for iRow in range(nx-1): + for iCur in range(nrComb): + tempDist = np.ones(nrComb) + for iPrev in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iPrev, iStrip+1]] + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iCur, iStrip]] + tempDist[iPrev] *= intPot[xDomain[iPrev,stripWidth-1], xDomain[iCur, stripWidth-1]] + if block1[iBlock, 0] > 0: + tempDist[iPrev] *= intPot[xDomain[iPrev,0], output[iChain, iRow, block1[iBlock, 0]-1]] + if block1[iBlock, stripWidth-1] < ny-1: + tempDist[iPrev] *= intPot[xDomain[iPrev,-1], output[iChain, iRow, block1[iBlock, stripWidth-1]+1]] + if iRow > 0: + tempDist[iPrev] *= messages[iPrev, iRow-1] + messages[iCur,iRow] = np.sum(tempDist) + normConstMessages[iRow] = np.sum(messages[:,iRow]) + messages[:,iRow] = messages[:,iRow] / normConstMessages[iRow] + # Column sum + tempDist = np.ones(nrComb) + for iCur in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + if block1[iBlock, 0] > 0: + tempDist[iCur] *= intPot[xDomain[iCur,0], output[iChain, nx-1, block1[ +iBlock, 0]-1]] + if block1[iBlock, -1] < ny-1: + tempDist[iCur] *= intPot[xDomain[iCur,-1], output[iChain, nx-1, block1[iBlock, -1]+1]] + tempDist[iCur] *= messages[iCur, nx-2] + normConstMessages[-1] = np.sum( tempDist ) + fMarginalA[iChain,iBlock] = np.sum( np.log(normConstMessages) ) + #raw_input() + + # Backward sampling + for iRow in range(nx)[::-1]: + tempDist = np.ones( nrComb ) + for iCur in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + if iRow < nx-1: + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], output[iChain,iRow+1, block1[iBlock, iStrip]]] + if iRow < nx-1: + tempDist[iCur] *= intPot[xDomain[iCur,-1], output[iChain,iRow+1, block1[iBlock, stripWidth-1]]] + if block1[iBlock, 0] > 0: + tempDist[iCur] *= intPot[xDomain[iCur,0], output[iChain, iRow, block1[iBlock, 0]-1]] + if block1[iBlock, stripWidth-1] < ny-1: + tempDist[iCur] *= intPot[xDomain[iCur,-1], output[iChain, iRow, block1[iBlock, stripWidth-1]+1]] + if iRow > 0: + tempDist[iCur] *= messages[iCur, iRow-1] + tempDist = tempDist / np.sum(tempDist) + curInd = discreteSampling( tempDist, range(nrComb), 1 ) + for iStrip in range(stripWidth): + output[iChain,iRow,block1[iBlock, iStrip]] = xDomain[curInd,iStrip] + + # -------------- + # Second block + # ------------- + for iBlock in range(nrInd2): + messages = np.ones( (nrComb, nx-1) ) + normConstMessages = np.zeros( nx ) + # Forward filtering + for iRow in range(nx-1): + for iCur in range(nrComb): + tempDist = np.ones(nrComb) + for iPrev in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iPrev, iStrip+1]] + tempDist[iPrev] *= intPot[xDomain[iPrev,iStrip], xDomain[iCur, iStrip]] + tempDist[iPrev] *= intPot[xDomain[iPrev,-1], xDomain[iCur, -1]] + if block2[iBlock, 0] > 0: + tempDist[iPrev] *= intPot[xDomain[iPrev,0], output[iChain, iRow, block2[iBlock, 0]-1]] + if block2[iBlock, -1] < ny-1: + tempDist[iPrev] *= intPot[xDomain[iPrev,-1], output[iChain, iRow, block2[iBlock, -1]+1]] + if iRow > 0: + tempDist[iPrev] *= messages[iPrev, iRow-1] + messages[iCur,iRow] = np.sum(tempDist) + normConstMessages[iRow] = np.sum(messages[:,iRow]) + messages[:,iRow] = messages[:,iRow] / normConstMessages[iRow] + + #print messages[:,iRow] + + tempDist = np.ones(nrComb) + for iCur in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + if block2[iBlock, 0] > 0: + tempDist[iCur] *= intPot[xDomain[iCur,0], output[iChain, nx-1, block2[iBlock, 0]-1]] + if block2[iBlock, -1] < ny-1: + tempDist[iCur] *= intPot[xDomain[iCur,-1], output[iChain, nx-1, block2[iBlock, -1]+1]] + tempDist[iCur] *= messages[iCur, nx-2] + normConstMessages[-1] = np.sum( tempDist ) + fMarginalB[iChain,iBlock] = np.sum( np.log(normConstMessages) ) + #print fMarginalB[iChain,iBlock] + #raw_input() + + # Backward sampling + for iRow in range(nx)[::-1]: + tempDist = np.ones( nrComb ) + for iCur in range(nrComb): + for iStrip in range(stripWidth-1): + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], xDomain[iCur, iStrip+1]] + if iRow < nx-1: + tempDist[iCur] *= intPot[xDomain[iCur,iStrip], output[iChain,iRow+1, block2[iBlock, iStrip]]] + if iRow < nx-1: + tempDist[iCur] *= intPot[xDomain[iCur,-1], output[iChain,iRow+1, block2[iBlock, -1]]] + if block2[iBlock, 0] > 0: + tempDist[iCur] *= intPot[xDomain[iCur,0], output[iChain, iRow, block2[iBlock, 0]-1]] + if block2[iBlock, -1] < ny-1: + tempDist[iCur] *= intPot[xDomain[iCur,-1], output[iChain, iRow, block2[iBlock, -1]+1]] + if iRow > 0: + tempDist[iCur] *= messages[iCur, iRow-1] + tempDist = tempDist / np.sum(tempDist) + curInd = discreteSampling( tempDist, range(nrComb), 1 ) + for iStrip in range(stripWidth): + output[iChain,iRow, block2[iBlock, iStrip]] = xDomain[curInd,iStrip] + #print fMarginal + #raw_input() + #print output[0,:,:].reshape( (nx,ny) ) + #raw_input() + fMarginal = np.zeros( (nrChains,nrInd1+nrInd2) ) + fMarginal[:,:len(block1)] = fMarginalA + fMarginal[:,len(block1):len(block1)+len(block2)] = fMarginalB + f = open(fileName, 'a') + f.write(str(iIter) + ' ' + str(time.time()-startTS) + ' ') + np.savetxt(f, fMarginal.reshape( (1,nrChains*(nrInd1+nrInd2)) ) ) + f.close() + +def discreteSampling(weights, domain, nrSamples): + bins = np.cumsum(weights) + return domain[np.digitize(random_sample(nrSamples), bins)] + +def ravel_multi_index(coord, shape): + return coord[0] * shape[1] + coord[1] + +def unravel_index(coord, shape): + iy = np.remainder(coord, shape[1]) + ix = (coord - iy) / shape[1] + return ix, iy + + +if __name__ == "__main__": + main()