Wolff-Ising-rough.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import base64\n",
    "try:\n",
    "    base64.encodebytes\n",
    "except AttributeError:\n",
    "    base64.encodebytes = base64.encodestring\n",
    "\n",
    "%pylab inline\n",
    "\n",
    "import pandas as pd\n",
    "import itertools as it\n",
    "from joblib import Parallel, delayed\n",
    "from matplotlib import animation as an\n",
    "#from IPython import display\n",
    "from tempfile import NamedTemporaryFile\n",
    "\n",
    "from numba import jit,b1, u2, i8, void\n",
    "#import numba as nb\n",
    "\n",
    "#@jit #(b1(b1[:,:], i8[:])) #, nopython=True)\n",
    "def aget(lattice, coords):\n",
    "    return lattice.item(*coords.tolist())\n",
    "\n",
    "#@jit #(void(b1[:,:], i8[:], b1)) #, nopython=True)\n",
    "def aset(lattice, coords, val):\n",
    "    l = coords.tolist()\n",
    "    l.append(val)\n",
    "    lattice.itemset(*l)\n",
    "\n",
    "@jit(nopython=True)\n",
    "def nextpos(current, dims, dim, delta):\n",
    "    nxt = current.copy()\n",
    "    nxt[dim] = (nxt[dim]+delta)%dims[dim]\n",
    "    return nxt\n",
    "\n",
    "#@jit#(nopython=True)\n",
    "def cluster(lattice, prob):\n",
    "    dims = np.array(lattice.shape)\n",
    "    start = mlab.amap(randint,dims)\n",
    "    state = not aget(lattice, start)\n",
    "    aset(lattice,start,state)\n",
    "    #print >>sys.stderr, \"$$$\",\n",
    "    ndim = lattice.ndim\n",
    "    n = 0\n",
    "    S = [start]\n",
    "    while S:\n",
    "        n+=1\n",
    "        current = S.pop()\n",
    "        for dim in range(ndim):\n",
    "            for delta in [-1,1]:\n",
    "                nxt = nextpos(current, dims, dim, delta)\n",
    "                if  (np.random.random() < prob) and (aget(lattice,nxt) ^ state):\n",
    "                    aset(lattice,nxt,state)\n",
    "                    #print >>sys.stderr, nxt,\n",
    "                    S.append(nxt)\n",
    "    #print >>sys.stderr\n",
    "    return n*(2*state-1)\n",
    "\n",
    "#@jit\n",
    "def run(lattice, prob):\n",
    "    yield lattice.sum() - lattice.size/2\n",
    "    while True:\n",
    "        yield cluster(lattice, prob)\n",
    "\n",
    "\n",
    "Tc = 2/log(1 + sqrt(2))\n",
    "Pc = 1 - 1/(1+sqrt(2))\n",
    "#xi = .3603 ## value given by ...\n",
    "Xi = 1.25 # value that actually seems optimal here\n",
    "\n",
    "def Tfin(dim, xi=Xi):\n",
    "    return Tc*(1 + float(xi)/dim)\n",
    "\n",
    "def Prob(T):\n",
    "    return 1-exp(-2/T)\n",
    "\n",
    "\n",
    "# In[7]:\n",
    "\n",
    "deltas0=[]\n",
    "d = 2**9\n",
    "A = zeros((d,d), dtype=bool)\n",
    "#A = np.random.random((d,d)) < 0.5\n",
    "Pfin = Prob(Tfin(d))\n",
    "print Pfin, A.size\n",
    "print Tc, Pc, Prob(Tc), Tfin(d), Prob(Tfin(d))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "\n",
    "############################\n",
    "\n",
    "# deltas0=fromiter(run(A, Pfin), int, A.size)\n",
    "# exit(0)\n",
    "\n",
    "###########################################################\n",
    "\n",
    "# def nrun(xi):\n",
    "#     AA = zeros((d,d), dtype=bool)\n",
    "#     return fromiter(run(AA, Prob(Tfin(d, xi))), int, 2*AA.size)\n",
    "\n",
    "# deltas = pd.DataFrame(Parallel(n_jobs=-1, verbose=3)(amap(delayed(nrun), linspace(0, 2.0, 201)))).T\n",
    "\n",
    "# pd.set_option(\"display.width\", 256)\n",
    "\n",
    "# print abs(deltas).describe()\n",
    "\n",
    "# print abs(deltas.cumsum()).describe()\n",
    "\n",
    "# print (deltas.cumsum()**2).describe()\n",
    "\n",
    "# deltas.to_csv(\"zzz2.csv\")\n",
    "# exit(0)\n",
    "\n",
    "###########################################################\n",
    "\n",
    "\n",
    "from sklearn import decomposition as dec\n",
    "def cmpr(a, k):\n",
    "    model = dec.NMF(k)\n",
    "    print model\n",
    "    print a\n",
    "    print ma.is_masked(a)\n",
    "    a = a.astype(float)\n",
    "    print a\n",
    "    print ma.is_masked(a)\n",
    "    a=model.fit_transform(a)\n",
    "    print a\n",
    "    print ma.is_masked(a)\n",
    "    return model.inverse_transform(a)\n",
    "    #return a.dot(model.components_)\n",
    "\n",
    "from scipy import ndimage\n",
    "from scipy import signal as sig\n",
    "from skimage import restoration as rst\n",
    "\n",
    "#@jit(nopython=True)\n",
    "def smooth(a):\n",
    "    a = a.astype(float)\n",
    "    #orig = a\n",
    "    #a = a + normal(0,2**-1,a.shape)\n",
    "    #a = a+standard_cauchy(a.shape)\n",
    "    #a = ndimage.median_filter(a, 2**1, mode='wrap')\n",
    "    #a = ndimage.uniform_filter(a, 2**1, mode='wrap')\n",
    "    #a = ndimage.fourier_uniform(a, 2**5)\n",
    "    #a = 1.1**a\n",
    "    #a = cmpr(a, 2**5)\n",
    "    #a = fftcmp2(a, 2**12)\n",
    "    a=denoise(a, 2**2, 'bior3.3')\n",
    "    #a = rst.denoise_bilateral(a, 2**5)\n",
    "    #a = sig.wiener(a,6)\n",
    "    #a = ndimage.median_filter(a, 2**1, mode='wrap')\n",
    "    #a = sig.hilbert2(a)\n",
    "    #err = a - orig\n",
    "    #print >>sys.stderr, a.min(), a.max(), a.mean(), a.std(), orig.std(), sqrt((err**2).mean())\n",
    "    #print >>sys.stderr, abs(orig - orig.mean()).mean(), abs(a - a.mean()).mean(), abs(err).mean()\n",
    "    #a -= a.mean()\n",
    "    #a -= 0.15\n",
    "    #a /= a.std()\n",
    "    return a\n",
    "\n",
    "import pywt\n",
    "\n",
    "from sklearn import preprocessing as pre\n",
    "\n",
    "#@jit(nopython=True)\n",
    "def denoise(data, noiseSigma, wavenm, mode='ppd'):\n",
    "    wavelet = pywt.Wavelet(wavenm)\n",
    "    WC = pywt.wavedec2(data,wavelet,mode=mode)\n",
    "    #print wavenm\n",
    "    #print data.mean(), data.std(), data.min(), data.max()\n",
    "    threshold=noiseSigma*sqrt(2*log2(data.size))\n",
    "    WC = [WC[0]] + map(lambda x: pywt.thresholding.soft(x,threshold), WC[1:])\n",
    "    #WC = [WC[0]] + map(lambda x: pywt.thresholding.hard(x,threshold), WC[1:])\n",
    "    #WC = [WC[0]] + map(lambda x: pywt.thresholding.greater(x,-threshold,-threshold), WC[1:])\n",
    "    #WC = [WC[0]] + map(lambda x: pywt.thresholding.less(x,threshold,threshold), WC[1:])\n",
    "    #WC = map(lambda x: pywt.thresholding.soft(x,threshold), WC)\n",
    "    #WC = map(lambda x: pywt.thresholding.hard(x,threshold), WC)\n",
    "    #WC = [WC[0]]+map(lambda x: (abs(asarray(x, float64))**noiseSigma) * sign(x), WC[1:])\n",
    "    #WC = map(lambda x: (abs(asarray(x, float64))**noiseSigma) * sign(x), WC)\n",
    "    #ZZZ = hstack(map(lambda x: ravel(asarray(x)), NWC))\n",
    "    #print >>sys.stderr, \">>>\", ZZZ.mean(), ZZZ.std(), ZZZ.min(), ZZZ.max()\n",
    "    #print >>sys.stderr, len(NWC), map(lambda x: x.shape, NWC)\n",
    "    #print WC\n",
    "    #print NWC\n",
    "    res = pywt.waverec2(WC, wavelet, mode=mode)\n",
    "    #print res.mean(), res.std(), res.min(), res.max()\n",
    "    #print\n",
    "    return pre.minmax_scale(pywt.waverec2(WC, wavelet, mode=mode).ravel()).reshape(data.shape)\n",
    "\n",
    "\n",
    "from scipy import fftpack\n",
    "def get_2D_dct(img):\n",
    "    \"\"\" Get 2D Cosine Transform of Image\n",
    "    \"\"\"\n",
    "    return fftpack.dct(fftpack.dct(img.T, norm='ortho').T, norm='ortho')\n",
    "\n",
    "def get_2d_idct(coefficients):\n",
    "    \"\"\" Get 2D Inverse Cosine Transform of Image\n",
    "    \"\"\"\n",
    "    return fftpack.idct(fftpack.idct(coefficients.T, norm='ortho').T, norm='ortho')\n",
    "\n",
    "def fftcmp(pixels, ii):\n",
    "    dct_size = pixels.shape[0]\n",
    "    dct = get_2D_dct(pixels)\n",
    "    reconstructed_images = []\n",
    "    dct_copy = dct.copy()\n",
    "    dct_copy[ii:,:] = 0\n",
    "    dct_copy[:,ii:] = 0\n",
    "    return get_2d_idct(dct_copy);\n",
    "\n",
    "def fftcmp2(pixels, k):\n",
    "    dct = get_2D_dct(pixels)\n",
    "    mask = ((-abs(dct.ravel())).argsort().argsort() < k).reshape(dct.shape)\n",
    "    dct2 = dct*mask.astype(float)\n",
    "    idct = get_2d_idct(dct2)\n",
    "    return idct\n",
    "\n",
    "\n",
    "import sys\n",
    "#@jit\n",
    "def animate(img, prob, frames=60):\n",
    "    lattice = img.get_array()\n",
    "    def init():\n",
    "        img.set_array(lattice)\n",
    "        return img\n",
    "    def frm(i):\n",
    "        cluster(lattice, prob)\n",
    "        img.set_array(smooth(lattice))\n",
    "        print >>sys.stderr, i,\n",
    "        return img\n",
    "    return an.FuncAnimation(img.get_figure(), init_func=init, func=frm, frames=frames, repeat=False, blit=False, interval=2**-20)\n",
    "\n",
    "\n",
    "\n",
    "r, c = A.shape\n",
    "dpi = gcf().dpi\n",
    "close()\n",
    "figure(figsize=(r/dpi, c/dpi))\n",
    "Img=figimage(A, cmap=\"gray\",\n",
    "             vmin=-0., vmax=1.)\n",
    "#Img=imshow(A,cmap=\"gray\", vmin=0.0, vmax=1.0)\n",
    "\n",
    "print dir(Img)\n",
    "crf=16\n",
    "fps=18\n",
    "secs = 4*60\n",
    "nfrm = secs*fps\n",
    "fs = 2**17\n",
    "bv = \"14k\"\n",
    "anim = animate(Img, Pfin, nfrm)\n",
    "show()\n",
    "\n",
    "# with NamedTemporaryFile(dir='.', prefix='isingwolff-d%s-xi%s-crf%s-fps%s-' % (d,Xi,crf, fps), suffix='.mp4', delete=False) as f:\n",
    "#     anim.save(f.name, fps=fps, dpi=Img.get_figure().get_dpi(), extra_args=['-vcodec', 'libx264', \"-crf\", str(crf)], writer='ffmpeg')\n",
    "\n",
    "# with NamedTemporaryFile(dir='.', prefix='isingwolff-d%s-xi%s-nfrm%s-fs%s-' % (d, Xi, nfrm, fs), suffix='.mp4', delete=False) as f:\n",
    "#     anim.save(f.name, fps=fps, dpi=Img.get_figure().get_dpi(), extra_args=['-vcodec', 'libx264', \"-fs\", str(fs)], writer='ffmpeg')\n",
    "\n",
    "# with NamedTemporaryFile(dir='.', prefix='isingwolff-d%s-xi%s-secs%s-bv%s-fps%s-' % (d, Xi, secs, bv, fps), suffix='.mp4', delete=False) as f:\n",
    "#     anim.save(f.name, fps=fps, dpi=Img.get_figure().get_dpi(), extra_args=['-vcodec', 'libx264', \"-b:v\", bv], writer='ffmpeg')\n",
    "\n",
    "exit(0)\n",
    "\n",
    "###########################################################\n",
    "\n",
    "\n",
    "print isinteractive()\n",
    "grp = pool.Pool(100)\n",
    "\n",
    "#imgx=figimage(A, cmap=\"gray\")\n",
    "ZZ=0\n",
    "#@jit\n",
    "#def vis():\n",
    "    #figure()\n",
    "    #fig, ax = plt.subplots()\n",
    "    #figure(figsize=(18,18))\n",
    "def vis(x):\n",
    "        global ZZ, imgx\n",
    "        if ZZ == 0:\n",
    "           print \".\",\n",
    "           A[0,0] = 1\n",
    "           # delaxes()\n",
    "           #ion()\n",
    "           imgx=figimage(A, cmap=\"gray\")\n",
    "           pause(2**-32)\n",
    "           #fig = plt.gcf()\n",
    "           #draw()\n",
    "           # ax = plt.Axes(fig, [0., 0., 1., 1.])\n",
    "           # ax.set_axis_off()\n",
    "           # fig.add_axes(ax)\n",
    "           #axis(\"off\")\n",
    "           #show()\n",
    "           #clim()   # clamp the color limits\n",
    "           #title(\"Boring slide show\")\n",
    "        else:\n",
    "           pass\n",
    "           imgx.set_array(A)\n",
    "           #draw()\n",
    "           #gevent.spawn(lambda : imgx.set_data(A)).start()\n",
    "           #suptitle(\"   \".join(map(str,[ZZ, x, A.sum()])))\n",
    "           if ZZ % 5 == 0: pause(10**-30)\n",
    "           #if ZZ % 5 == 0: grp.apply(pause,[2**-32])\n",
    "           #gevent.spawn(lambda : pause(2**-32)).start()\n",
    "        ZZ+=1\n",
    "        return x\n",
    "#return vis2\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#figure(figsize=(18,18))\n",
    "#axis(\"off\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "deltas0=fromiter(it.imap(vis, run(A, Prob(Tfin(d)))), int, A.size)\n",
    "\n",
    "exit(0)\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#%time anim=animate(A,Pfin,deltas0.append,d**2)\n",
    "#%time display_animation(anim)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "deltas=array(deltas0)\n",
    "sig=cumsum(deltas)\n",
    "fluct=abs(deltas)\n",
    "\n",
    "pwr=abs(sig) # + A.size/2\n",
    "#plot(signal.medfilt(sig/(d**2.),1))\n",
    "#plot(pwr)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(5, 5))\n",
    "\n",
    "plot(sorted(log2(fluct), reverse=True))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "plot(log2(arange(len(pwr))), log2(sorted((pwr), reverse=True)))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "plot(log2(bincount(fluct/(d))))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "flucttbl=bincount(fluct/(d))\n",
    "figure(figsize=(7,7))\n",
    "plot(log2(arange(len(flucttbl))+0), log2(flucttbl))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "pwrtbl=bincount((pwr/d))\n",
    "plot((arange(pwrtbl.size)), (pwrtbl))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(15,3))\n",
    "plot(deltas[-1000:])\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(15,3))\n",
    "plot(sig[-1000:]**2)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "scatter(arcsinh((sig[:-1])/float(d**1.3)),arcsinh((sig[1:])/float(d**1.3)),alpha=2**-4, marker=\".\")\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "scatter(arcsinh(deltas[1:]/float(d)),arcsinh(deltas[:-1]/float(d)),alpha=2**-4, marker=\".\")\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "scatter((sig/float(d**0)),(deltas/float(d**0)),alpha=2**-3, marker=\".\")\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "from scipy import signal\n",
    "#plot(log2(arange(len(deltas))+1), signal.medfilt((abs(cumsum(deltas))), 8*d-1))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#plot(signal.medfilt((fluct), 8*d-1))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#plot(log2(arange(len(fluct))+1),signal.medfilt(log2(fluct),8*d-1))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "plot(log2(bincount(amap(lambda x: int(log2(x)), fluct[len(fluct)/2:]))))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#figure(figsize(10,10))\n",
    "imshow(A)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "z=psd((sig[d:]**2)-mean(sig[d:]**2),NFFT=d, noverlap=d/2)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "plot(log2(z[1]),log2(z[0]))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "from scipy import stats\n",
    "print stats.kurtosis(vstack([deltas,fluct,sig,pwr, pwr**2]), axis = 1, bias=False)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "print std(vstack([deltas,fluct,sig,pwr, pwr**2])/double(A.size), axis = 1)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "print 2*median(pwr)/double(A.size)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figimage\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "subplots()[0]\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "gcf().figimage\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "x=figimage(A);\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "x\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "y=imshow(A)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "y\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "#from gevent import monkey\n",
    "#monkey.patch_all()\n",
    "#import gevent\n",
    "#from gevent import pool\n",
    "\n",
    "import base64\n",
    "try:\n",
    "    base64.encodebytes\n",
    "except AttributeError:\n",
    "    base64.encodebytes = base64.encodestring\n",
    "\n",
    "# In[2]:\n",
    "\n",
    "from pylab import *\n",
    "import pandas as pd\n",
    "import itertools as it\n",
    "from joblib import Parallel, delayed\n",
    "from matplotlib import animation as an\n",
    "#from IPython import display\n",
    "from tempfile import NamedTemporaryFile\n",
    "\n",
    "\n",
    "# ## Video emedding\n",
    "\n",
    "# In[3]:\n",
    "\n",
    "from numba import jit,b1, u2, i8, void\n",
    "#import numba as nb\n",
    "\n",
    "#@jit #(b1(b1[:,:], i8[:])) #, nopython=True)\n",
    "def aget(lattice, coords):\n",
    "    return lattice.item(*coords.tolist())\n",
    "\n",
    "#@jit #(void(b1[:,:], i8[:], b1)) #, nopython=True)\n",
    "def aset(lattice, coords, val):\n",
    "    l = coords.tolist()\n",
    "    l.append(val)\n",
    "    lattice.itemset(*l)\n",
    "\n",
    "@jit(nopython=True)\n",
    "def nextpos(current, dims, dim, delta):\n",
    "    nxt = current.copy()\n",
    "    nxt[dim] = (nxt[dim]+delta)%dims[dim]\n",
    "    return nxt\n",
    "\n",
    "#@jit#(nopython=True)\n",
    "def cluster(lattice, prob):\n",
    "    dims = np.array(lattice.shape)\n",
    "    start = mlab.amap(randint,dims)\n",
    "    state = not aget(lattice, start)\n",
    "    aset(lattice,start,state)\n",
    "    #print >>sys.stderr, \"$$$\",\n",
    "    ndim = lattice.ndim\n",
    "    n = 0\n",
    "    S = [start]\n",
    "    while S:\n",
    "        n+=1\n",
    "        current = S.pop()\n",
    "        for dim in range(ndim):\n",
    "            for delta in [-1,1]:\n",
    "                nxt = nextpos(current, dims, dim, delta)\n",
    "                if  (np.random.random() < prob) and (aget(lattice,nxt) ^ state):\n",
    "                    aset(lattice,nxt,state)\n",
    "                    #print >>sys.stderr, nxt,\n",
    "                    S.append(nxt)\n",
    "    #print >>sys.stderr\n",
    "    return n*(2*state-1)\n",
    "\n",
    "#@jit\n",
    "def run(lattice, prob):\n",
    "    yield lattice.sum() - lattice.size/2\n",
    "    while True:\n",
    "        yield cluster(lattice, prob)\n",
    "\n",
    "\n",
    "# ## Finite lattice corrrections\n",
    "\n",
    "# In[6]:\n",
    "\n",
    "Tc = 2/log(1 + sqrt(2))\n",
    "Pc = 1 - 1/(1+sqrt(2))\n",
    "#xi = .3603 ## value given by ...\n",
    "Xi = 1.25 # value that actually seems optimal here\n",
    "\n",
    "def Tfin(dim, xi=Xi):\n",
    "    return Tc*(1 + float(xi)/dim)\n",
    "\n",
    "def Prob(T):\n",
    "    return 1-exp(-2/T)\n",
    "\n",
    "\n",
    "# In[7]:\n",
    "\n",
    "deltas0=[]\n",
    "d = 2**9\n",
    "A = zeros((d,d), dtype=bool)\n",
    "#A = np.random.random((d,d)) < 0.5\n",
    "Pfin = Prob(Tfin(d))\n",
    "print Pfin, A.size\n",
    "print Tc, Pc, Prob(Tc), Tfin(d), Prob(Tfin(d))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "\n",
    "############################\n",
    "\n",
    "# deltas0=fromiter(run(A, Pfin), int, A.size)\n",
    "# exit(0)\n",
    "\n",
    "###########################################################\n",
    "\n",
    "# def nrun(xi):\n",
    "#     AA = zeros((d,d), dtype=bool)\n",
    "#     return fromiter(run(AA, Prob(Tfin(d, xi))), int, 2*AA.size)\n",
    "\n",
    "# deltas = pd.DataFrame(Parallel(n_jobs=-1, verbose=3)(amap(delayed(nrun), linspace(0, 2.0, 201)))).T\n",
    "\n",
    "# pd.set_option(\"display.width\", 256)\n",
    "\n",
    "# print abs(deltas).describe()\n",
    "\n",
    "# print abs(deltas.cumsum()).describe()\n",
    "\n",
    "# print (deltas.cumsum()**2).describe()\n",
    "\n",
    "# deltas.to_csv(\"zzz2.csv\")\n",
    "# exit(0)\n",
    "\n",
    "###########################################################\n",
    "\n",
    "\n",
    "from sklearn import decomposition as dec\n",
    "def cmpr(a, k):\n",
    "    model = dec.NMF(k)\n",
    "    print model\n",
    "    print a\n",
    "    print ma.is_masked(a)\n",
    "    a = a.astype(float)\n",
    "    print a\n",
    "    print ma.is_masked(a)\n",
    "    a=model.fit_transform(a)\n",
    "    print a\n",
    "    print ma.is_masked(a)\n",
    "    return model.inverse_transform(a)\n",
    "    #return a.dot(model.components_)\n",
    "\n",
    "from scipy import ndimage\n",
    "from scipy import signal as sig\n",
    "from skimage import restoration as rst\n",
    "\n",
    "#@jit(nopython=True)\n",
    "def smooth(a):\n",
    "    a = a.astype(float)\n",
    "    #orig = a\n",
    "    #a = a + normal(0,2**-1,a.shape)\n",
    "    #a = a+standard_cauchy(a.shape)\n",
    "    #a = ndimage.median_filter(a, 2**1, mode='wrap')\n",
    "    #a = ndimage.uniform_filter(a, 2**1, mode='wrap')\n",
    "    #a = ndimage.fourier_uniform(a, 2**5)\n",
    "    #a = 1.1**a\n",
    "    #a = cmpr(a, 2**5)\n",
    "    #a = fftcmp2(a, 2**12)\n",
    "    a=denoise(a, 2**2, 'bior3.3')\n",
    "    #a = rst.denoise_bilateral(a, 2**5)\n",
    "    #a = sig.wiener(a,6)\n",
    "    #a = ndimage.median_filter(a, 2**1, mode='wrap')\n",
    "    #a = sig.hilbert2(a)\n",
    "    #err = a - orig\n",
    "    #print >>sys.stderr, a.min(), a.max(), a.mean(), a.std(), orig.std(), sqrt((err**2).mean())\n",
    "    #print >>sys.stderr, abs(orig - orig.mean()).mean(), abs(a - a.mean()).mean(), abs(err).mean()\n",
    "    #a -= a.mean()\n",
    "    #a -= 0.15\n",
    "    #a /= a.std()\n",
    "    return a\n",
    "\n",
    "import pywt\n",
    "\n",
    "from sklearn import preprocessing as pre\n",
    "\n",
    "#@jit(nopython=True)\n",
    "def denoise(data, noiseSigma, wavenm, mode='ppd'):\n",
    "    wavelet = pywt.Wavelet(wavenm)\n",
    "    WC = pywt.wavedec2(data,wavelet,mode=mode)\n",
    "    #print wavenm\n",
    "    #print data.mean(), data.std(), data.min(), data.max()\n",
    "    threshold=noiseSigma*sqrt(2*log2(data.size))\n",
    "    WC = [WC[0]] + map(lambda x: pywt.thresholding.soft(x,threshold), WC[1:])\n",
    "    #WC = [WC[0]] + map(lambda x: pywt.thresholding.hard(x,threshold), WC[1:])\n",
    "    #WC = [WC[0]] + map(lambda x: pywt.thresholding.greater(x,-threshold,-threshold), WC[1:])\n",
    "    #WC = [WC[0]] + map(lambda x: pywt.thresholding.less(x,threshold,threshold), WC[1:])\n",
    "    #WC = map(lambda x: pywt.thresholding.soft(x,threshold), WC)\n",
    "    #WC = map(lambda x: pywt.thresholding.hard(x,threshold), WC)\n",
    "    #WC = [WC[0]]+map(lambda x: (abs(asarray(x, float64))**noiseSigma) * sign(x), WC[1:])\n",
    "    #WC = map(lambda x: (abs(asarray(x, float64))**noiseSigma) * sign(x), WC)\n",
    "    #ZZZ = hstack(map(lambda x: ravel(asarray(x)), NWC))\n",
    "    #print >>sys.stderr, \">>>\", ZZZ.mean(), ZZZ.std(), ZZZ.min(), ZZZ.max()\n",
    "    #print >>sys.stderr, len(NWC), map(lambda x: x.shape, NWC)\n",
    "    #print WC\n",
    "    #print NWC\n",
    "    res = pywt.waverec2(WC, wavelet, mode=mode)\n",
    "    #print res.mean(), res.std(), res.min(), res.max()\n",
    "    #print\n",
    "    return pre.minmax_scale(pywt.waverec2(WC, wavelet, mode=mode).ravel()).reshape(data.shape)\n",
    "\n",
    "\n",
    "from scipy import fftpack\n",
    "def get_2D_dct(img):\n",
    "    \"\"\" Get 2D Cosine Transform of Image\n",
    "    \"\"\"\n",
    "    return fftpack.dct(fftpack.dct(img.T, norm='ortho').T, norm='ortho')\n",
    "\n",
    "def get_2d_idct(coefficients):\n",
    "    \"\"\" Get 2D Inverse Cosine Transform of Image\n",
    "    \"\"\"\n",
    "    return fftpack.idct(fftpack.idct(coefficients.T, norm='ortho').T, norm='ortho')\n",
    "\n",
    "def fftcmp(pixels, ii):\n",
    "    dct_size = pixels.shape[0]\n",
    "    dct = get_2D_dct(pixels)\n",
    "    reconstructed_images = []\n",
    "    dct_copy = dct.copy()\n",
    "    dct_copy[ii:,:] = 0\n",
    "    dct_copy[:,ii:] = 0\n",
    "    return get_2d_idct(dct_copy);\n",
    "\n",
    "def fftcmp2(pixels, k):\n",
    "    dct = get_2D_dct(pixels)\n",
    "    mask = ((-abs(dct.ravel())).argsort().argsort() < k).reshape(dct.shape)\n",
    "    dct2 = dct*mask.astype(float)\n",
    "    idct = get_2d_idct(dct2)\n",
    "    return idct\n",
    "\n",
    "\n",
    "import sys\n",
    "#@jit\n",
    "def animate(img, prob, frames=60):\n",
    "    lattice = img.get_array()\n",
    "    def init():\n",
    "        img.set_array(lattice)\n",
    "        return img\n",
    "    def frm(i):\n",
    "        cluster(lattice, prob)\n",
    "        img.set_array(smooth(lattice))\n",
    "        print >>sys.stderr, i,\n",
    "        return img\n",
    "    return an.FuncAnimation(img.get_figure(), init_func=init, func=frm, frames=frames, repeat=False, blit=False, interval=2**-20)\n",
    "\n",
    "\n",
    "\n",
    "r, c = A.shape\n",
    "dpi = gcf().dpi\n",
    "close()\n",
    "figure(figsize=(r/dpi, c/dpi))\n",
    "Img=figimage(A, cmap=\"gray\",\n",
    "             vmin=-0., vmax=1.)\n",
    "#Img=imshow(A,cmap=\"gray\", vmin=0.0, vmax=1.0)\n",
    "\n",
    "print dir(Img)\n",
    "crf=16\n",
    "fps=18\n",
    "secs = 4*60\n",
    "nfrm = secs*fps\n",
    "fs = 2**17\n",
    "bv = \"14k\"\n",
    "anim = animate(Img, Pfin, nfrm)\n",
    "show()\n",
    "\n",
    "# with NamedTemporaryFile(dir='.', prefix='isingwolff-d%s-xi%s-crf%s-fps%s-' % (d,Xi,crf, fps), suffix='.mp4', delete=False) as f:\n",
    "#     anim.save(f.name, fps=fps, dpi=Img.get_figure().get_dpi(), extra_args=['-vcodec', 'libx264', \"-crf\", str(crf)], writer='ffmpeg')\n",
    "\n",
    "# with NamedTemporaryFile(dir='.', prefix='isingwolff-d%s-xi%s-nfrm%s-fs%s-' % (d, Xi, nfrm, fs), suffix='.mp4', delete=False) as f:\n",
    "#     anim.save(f.name, fps=fps, dpi=Img.get_figure().get_dpi(), extra_args=['-vcodec', 'libx264', \"-fs\", str(fs)], writer='ffmpeg')\n",
    "\n",
    "# with NamedTemporaryFile(dir='.', prefix='isingwolff-d%s-xi%s-secs%s-bv%s-fps%s-' % (d, Xi, secs, bv, fps), suffix='.mp4', delete=False) as f:\n",
    "#     anim.save(f.name, fps=fps, dpi=Img.get_figure().get_dpi(), extra_args=['-vcodec', 'libx264', \"-b:v\", bv], writer='ffmpeg')\n",
    "\n",
    "exit(0)\n",
    "\n",
    "###########################################################\n",
    "\n",
    "\n",
    "print isinteractive()\n",
    "grp = pool.Pool(100)\n",
    "\n",
    "#imgx=figimage(A, cmap=\"gray\")\n",
    "ZZ=0\n",
    "#@jit\n",
    "#def vis():\n",
    "    #figure()\n",
    "    #fig, ax = plt.subplots()\n",
    "    #figure(figsize=(18,18))\n",
    "def vis(x):\n",
    "        global ZZ, imgx\n",
    "        if ZZ == 0:\n",
    "           print \".\",\n",
    "           A[0,0] = 1\n",
    "           # delaxes()\n",
    "           #ion()\n",
    "           imgx=figimage(A, cmap=\"gray\")\n",
    "           pause(2**-32)\n",
    "           #fig = plt.gcf()\n",
    "           #draw()\n",
    "           # ax = plt.Axes(fig, [0., 0., 1., 1.])\n",
    "           # ax.set_axis_off()\n",
    "           # fig.add_axes(ax)\n",
    "           #axis(\"off\")\n",
    "           #show()\n",
    "           #clim()   # clamp the color limits\n",
    "           #title(\"Boring slide show\")\n",
    "        else:\n",
    "           pass\n",
    "           imgx.set_array(A)\n",
    "           #draw()\n",
    "           #gevent.spawn(lambda : imgx.set_data(A)).start()\n",
    "           #suptitle(\"   \".join(map(str,[ZZ, x, A.sum()])))\n",
    "           if ZZ % 5 == 0: pause(10**-30)\n",
    "           #if ZZ % 5 == 0: grp.apply(pause,[2**-32])\n",
    "           #gevent.spawn(lambda : pause(2**-32)).start()\n",
    "        ZZ+=1\n",
    "        return x\n",
    "#return vis2\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#figure(figsize=(18,18))\n",
    "#axis(\"off\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "deltas0=fromiter(it.imap(vis, run(A, Prob(Tfin(d)))), int, A.size)\n",
    "\n",
    "exit(0)\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#%time anim=animate(A,Pfin,deltas0.append,d**2)\n",
    "#%time display_animation(anim)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "deltas=array(deltas0)\n",
    "sig=cumsum(deltas)\n",
    "fluct=abs(deltas)\n",
    "\n",
    "pwr=abs(sig) # + A.size/2\n",
    "#plot(signal.medfilt(sig/(d**2.),1))\n",
    "#plot(pwr)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(5, 5))\n",
    "\n",
    "plot(sorted(log2(fluct), reverse=True))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "plot(log2(arange(len(pwr))), log2(sorted((pwr), reverse=True)))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "plot(log2(bincount(fluct/(d))))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "flucttbl=bincount(fluct/(d))\n",
    "figure(figsize=(7,7))\n",
    "plot(log2(arange(len(flucttbl))+0), log2(flucttbl))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "pwrtbl=bincount((pwr/d))\n",
    "plot((arange(pwrtbl.size)), (pwrtbl))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(15,3))\n",
    "plot(deltas[-1000:])\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(15,3))\n",
    "plot(sig[-1000:]**2)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "scatter(arcsinh((sig[:-1])/float(d**1.3)),arcsinh((sig[1:])/float(d**1.3)),alpha=2**-4, marker=\".\")\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "scatter(arcsinh(deltas[1:]/float(d)),arcsinh(deltas[:-1]/float(d)),alpha=2**-4, marker=\".\")\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "scatter((sig/float(d**0)),(deltas/float(d**0)),alpha=2**-3, marker=\".\")\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "from scipy import signal\n",
    "#plot(log2(arange(len(deltas))+1), signal.medfilt((abs(cumsum(deltas))), 8*d-1))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#plot(signal.medfilt((fluct), 8*d-1))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#plot(log2(arange(len(fluct))+1),signal.medfilt(log2(fluct),8*d-1))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "plot(log2(bincount(amap(lambda x: int(log2(x)), fluct[len(fluct)/2:]))))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "#figure(figsize(10,10))\n",
    "imshow(A)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "z=psd((sig[d:]**2)-mean(sig[d:]**2),NFFT=d, noverlap=d/2)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figure(figsize=(7,7))\n",
    "plot(log2(z[1]),log2(z[0]))\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "from scipy import stats\n",
    "print stats.kurtosis(vstack([deltas,fluct,sig,pwr, pwr**2]), axis = 1, bias=False)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "print std(vstack([deltas,fluct,sig,pwr, pwr**2])/double(A.size), axis = 1)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "print 2*median(pwr)/double(A.size)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "figimage\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "subplots()[0]\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "gcf().figimage\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "x=figimage(A);\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "x\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "y=imshow(A)\n",
    "\n",
    "\n",
    "# In[ ]:\n",
    "\n",
    "y\n",
    "\n",
    "\n",
    "# In[ ]:\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}