Working copy of 30 May, 2016

Author: gdurin

Diff for: contours.py

@@ -0,0 +1,2 @@
+import numpy as np
+

Diff for: getAxyDistributions.py

@@ -0,0 +1,112 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar 20 10:59:08 2013
+
+@author: gf
+"""
+
+import os
+import pickle
+import numpy as np
+import itertools
+import getLogDistributions as gLD
+reload(gLD)
+import matplotlib.pyplot as plt
+import tables
+
+NN = {1:(1,1),
+      2:(2,1),
+      3:(3,1),
+      4:(2,2),
+      5:(3,2),
+      6:(3,2),
+      7:(3,3),      
+      8:(3,3),
+      9:(3,3),      
+      10:(5,2),
+      11:(4,3),
+      12:(4,3),
+      13:(4,4),
+      14:(4,4),
+      15:(5,3),
+      16:(4,4)
+      }
+
+mainDir = "/run/media/gf/DATA/meas/Barkh/Films/NiO_Fe/"
+f5filename = os.path.join(mainDir, "NiO_Fe_results.h5")
+
+f5 = tables.openFile(f5filename, mode="r")
+
+runNode = f5.root.NiO80.M10x.run21
+
+area = 1000.*1000.
+log_step = 1./10
+
+Axy_labels = []
+# Collect Axy labels first
+#for seqNode in f5.listNodes(runNode):
+    #AxyNode = seqNode.clus
+    #for label in f5.listNodes(AxyNode):
+        #Axy_labels.append(label.name)
+#Axy_labels = set(Axy_labels)
+
+aval_clus = 'clus'
+Axy = {}
+Axy[aval_clus] = {}
+
+for seqNode in f5.listNodes(runNode):
+    clusNode = f5.getNode(seqNode, aval_clus)
+    for lb in f5.listNodes(clusNode):
+        Axy[aval_clus][lb.name] = np.concatenate((
+            Axy[aval_clus].get(lb,np.array([])),lb.read()))
+f5.close()
+
+fig = plt.figure()
+fig.set_size_inches((20,10),forward=True)
+rows, cols = NN[len(Axy[aval_clus])]
+for n, key in enumerate(Axy[aval_clus]):
+    selected_Axy = Axy[aval_clus][key]
+    print "==============="
+    print "Cluster of type %s" % key
+    N_S, bins = np.histogram(selected_Axy, np.arange(1, max(selected_Axy)+2), 
+                             normed=False)
+    S = bins[:-1]
+    # Compress the data to non-zero values
+    #cond = N_S!=0
+    #N_S = np.compress(cond, N_S)
+    #S = np.compress(cond, S)i
+    yAverage = []
+    # Get log bins
+    SlogBins, logBins = gLD.getLogBins(bins[0], bins[-1], log_step)
+    for i, j in zip(logBins[:-1],logBins[1:]):
+        q1, q2 = np.greater_equal(S, i), np.less(S, j)
+        q = np.logical_and(q1, q2)
+        if sum(q) == 0:
+            averageValue = np.NaN
+        else:
+            #allElements = [val for val in itertools.chain(*A_S[q])]
+            averageValue = np.sum(S[q]*N_S[q])/area/(j-i)
+        yAverage.append(averageValue)
+    yAverage =  np.asanyarray(yAverage)
+    # Check if there are NaN values
+    iNan = np.isnan(yAverage)
+    x = SlogBins[~iNan]
+    y = yAverage[~iNan]
+    
+    plt.subplot(rows, cols, n+1)
+    plt.loglog(x, y, 'o', label=key)
+    # Plot the fit
+    q = np.logical_and(np.greater_equal(x, 10), np.less(x, 1000))
+    if sum(q)!=0:
+        X, Y = np.log10(x[q]), np.log10(y[q])
+        m, k = np.polyfit(X, Y, 1)
+        lin_x = np.linspace(10,1000)
+        new_y = 10**k*lin_x**m
+        plt.plot(lin_x, new_y, '-')
+        print "exponent: %f" % m
+    x, PS, PSerror = gLD.logDistribution(selected_Axy, log_step, normed=False)
+    plt.loglog(x, x*PS/area, '^')
+    plt.loglog(x, PS, 'v', label="P(S)")
+    plt.errorbar(x, PS, PSerror,fmt=None)
+    plt.legend()
+plt.show()

Diff for: getAxyLabels.py

@@ -47,7 +47,7 @@ def getAxyLabels(labels, imageDir="Left_to_right", edgeThickness=1, fraction=Non
     elif imageDir == "Top_to_bottom":
         lrbt = top, bottom, left, right
     else:
-        raise ValueError, "avalanche direction not defined"
+        raise ValueError("avalanche direction not defined")
     # Find the sets
     # Select the labels contained in the edges
     #setLrbt = set(left+right+bottom+top)
@@ -82,7 +82,7 @@ def getAxyLabels(labels, imageDir="Left_to_right", edgeThickness=1, fraction=Non
     #a = a.repeat(100,axis=0)
     #a = a.repeat(66,axis=1)
     labels, n = nd.label(a, structure)
-    print labels
+    print(labels)
     list_sizes = nd.sum(a, labels, range(1,n+1))
     array_sizes = scipy.array(list_sizes,dtype='int16')
     array_Axy = getAxyLabels(labels,'Bottom_to_top', edgeThickness=1, fraction=0.5)
@@ -91,5 +91,5 @@ def getAxyLabels(labels, imageDir="Left_to_right", edgeThickness=1, fraction=Non
         d[Axy] = scipy.concatenate((d.get(Axy,a0),sizes))
 
     for i, results in enumerate(zip(array_sizes,array_Axy)):
-        print i+1, results
-    print "Done in %.3f s" % (time.time()-startTime)
+        print(i+1, results)
+    print("Done in %.3f s" % (time.time()-startTime))

Diff for: getLogDistributions.py

@@ -48,6 +48,11 @@ def logDistribution(listValues, log_step=0.2, first_point=None,
                     last_point=None, normed=True):
     """
     Calculate the distribution in log scale from a list of values
+    Returns:
+    --------------
+    X : array of the x values at the center (in log-scale) of the bin
+    Y : array of the y values at the center (in log-scale) of the bin
+    Yerr :  array of the y errors, calculated by a binomial distribution
     """
     # Check the list of Values
     if not checkIfVoid(listValues):
@@ -56,13 +61,17 @@ def logDistribution(listValues, log_step=0.2, first_point=None,
         first_point = scipy.amin(listValues)
     if not last_point:
         last_point = scipy.amax(listValues)
-    xbins, bins = getLogBins(first_point, last_point, log_step)
-    yhist = scipy.stats.stats.histogram2(listValues, bins)
-    deltas = bins[1:]-bins[:-1]
-    yhist = yhist[:-1]/deltas
+    X, Xbins = getLogBins(first_point, last_point, log_step)
+    Y = scipy.stats.stats.histogram2(listValues, Xbins)
+    deltas = Xbins[1:]-Xbins[:-1]        # Check if zeros occurs
+    bool0 = Y!=0
+    for ary in [X,Y,deltas]:
+        ary = scipy.compress(bool0, ary)      
+    Y = Y[:-1]/deltas
+    Yerr = (Y * (1. - Y/sum(Y)))**0.5 / (deltas*sum(Y))
     if normed:
-        yhist = yhist/scipy.sum(yhist)
-    return xbins, yhist
+        Y = Y/scipy.sum(Y)
+    return X, Y, Yerr
 
 def averageLogDistribution(values, log_step=0.2, density=False,
                            first_point=None, last_point=None):

Diff for: getRofImages.py

@@ -0,0 +1,70 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 19 09:06:32 2013
+
+@author: gf
+"""
+
+import os, re, glob, sys
+import rof
+import numpy as np
+import Image
+import scipy.misc as misc
+
+def denoise(im):
+    U,T = rof.denoise(im,im)
+    return np.asarray(U, dtype='int32')
+    
+
+import skimage.io as im_io
+
+class Imread_convert():
+
+    def __call__(self, f):
+        im = np.array(Image.open(f), dtype='int32')
+        #            imageList = list(im.getdata())
+        #            sizeX, sizeY = im.size
+        #            return np.asanyarray(imageList).reshape(sizeY, sizeX)
+        return im
+
+
+def natural_key(string_):
+    """See http://www.codinghorror.com/blog/archives/001018.html"""
+    return [int(s) if s.isdigit() else s for s in re.split(r'(\d+)', string_)]
+    
+mainDir, pattern = "/media/DATA/meas/Barkh/Films/NiO-Fe/NiO80/NiO80 run19 10x pinLin 10iter/", "Data1-*.tif"
+
+s = "(%s|%s)" % tuple(pattern.split("*"))
+patternCompiled = re.compile(s)
+# Load all the image filenames
+imageFileNames = glob.glob1(mainDir, pattern)
+# Sort it with natural keys
+imageFileNames.sort(key=natural_key)
+
+if not len(imageFileNames):
+    print "ERROR, no images in %s" % mainDir
+    sys.exit()
+else:
+    print "Found %d images in %s" % (len(imageFileNames), mainDir)
+    
+# Prepare the rof subdir
+rofDir = os.path.join(mainDir, 'rof')
+if not os.path.exists(rofDir):
+    os.mkdir(rofDir)
+
+    
+imread_convert = Imread_convert()
+# Load the images
+print "Loading images: "
+load_pattern = [os.path.join(mainDir,ifn) for ifn in imageFileNames]
+
+imageCollection = im_io.ImageCollection(load_pattern, load_func=imread_convert)
+
+for i, im in enumerate(imageCollection):
+    print "Image n. %i" % i
+    U = denoise(im)
+    fileName = imageFileNames[i]
+    fileName = os.path.join(mainDir, 'rof', fileName)
+    q = misc.toimage(U)
+    q.save(fileName)
+    

Diff for: gpuSwitchtime.py

@@ -193,7 +193,7 @@ def get_gpuSwitchTime(stackImages, kernel, device=0):
     dim_z = 80
     a = np.random.randn(dim_z,dim_y,dim_x)
     a = a.astype(np.int32)
-    print "Loading %.2f MB of data" % (a.nbytes/1e6)
+    print("Loading %.2f MB of data" % (a.nbytes/1e6))
     # Call the GPU kernel
     kernel = np.array([-1]*5+[1]*5)
     t0 = time.time()
@@ -203,16 +203,15 @@ def get_gpuSwitchTime(stackImages, kernel, device=0):
     step = kernel
     cpuswitch=np.zeros((dim_y,dim_x),dtype=np.int32)
     cpulevels=np.zeros((dim_y,dim_x),dtype=np.int32)
-    print "Loading %.2f MB of data" % (2*cpuswitch.nbytes/1e6)
+    print("Loading %.2f MB of data" % (2*cpuswitch.nbytes/1e6))
     t3=time.time()
     for i in range(0,dim_x):
         for j in range(0,dim_y):
             indice=(nd.convolve1d(a[:,j,i],step,mode='reflect')).argmin()      
             cpuswitch[j,i]=indice
     timeCpu = time.time()-t3
-    print ("GPU calculus done = %.4f s" %timeGpu)
-    print ("CPU calculus done = %.4f s" %timeCpu)
-    print "Difference on switch : \n"
-    print gpuswitch-cpuswitch
-    
-    print ("\nGPU is %d times faster than CPU " %(timeCpu/timeGpu))
+    print("GPU calculus done = %.4f s" %timeGpu)
+    print("CPU calculus done = %.4f s" %timeCpu)
+    print("Difference on switch : \n")
+    print(gpuswitch-cpuswitch)
+    print("\nGPU is %d times faster than CPU " %(timeCpu/timeGpu))