Statistiques
| Révision :

root / Ising / Cython / Ising2D.py @ 98

Historique | Voir | Annoter | Télécharger (3,95 ko)

1 18 equemene
#!/usr/bin/env python
2 18 equemene
#
3 18 equemene
# Ising2D model in serial mode
4 18 equemene
#
5 18 equemene
# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr>
6 18 equemene
7 18 equemene
import sys
8 18 equemene
import numpy
9 18 equemene
from PIL import Image
10 18 equemene
from math import exp
11 18 equemene
from random import random
12 18 equemene
import time
13 18 equemene
import getopt
14 18 equemene
import matplotlib.pyplot as plt
15 18 equemene
import Metropolis
16 18 equemene
from Metropolis import Metropolis
17 18 equemene
18 18 equemene
def ImageOutput(sigma,prefix):
19 18 equemene
    Max=sigma.max()
20 18 equemene
    Min=sigma.min()
21 18 equemene
22 18 equemene
    # Normalize value as 8bits Integer
23 18 equemene
    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
24 18 equemene
    image = Image.fromarray(SigmaInt)
25 18 equemene
    image.save("%s.jpg" % prefix)
26 18 equemene
27 18 equemene
def Magnetization(sigma,M):
28 18 equemene
    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
29 18 equemene
30 18 equemene
def Energy(sigma,J):
31 18 equemene
    # Copier et caster
32 18 equemene
    E=numpy.copy(sigma).astype(numpy.float32)
33 18 equemene
34 18 equemene
    # Appel par slice
35 18 equemene
    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
36 18 equemene
                                  E[1:-1,:-2]+E[1:-1,2:])
37 18 equemene
38 18 equemene
    # Bien nettoyer la peripherie
39 18 equemene
    E[:,0]=0
40 18 equemene
    E[:,-1]=0
41 18 equemene
    E[0,:]=0
42 18 equemene
    E[-1,:]=0
43 18 equemene
44 18 equemene
    Energy=numpy.sum(E)
45 18 equemene
46 18 equemene
    return(Energy/(E.shape[0]*E.shape[1]*1.0))
47 18 equemene
48 18 equemene
def DisplayCurves(T,E,M,J,B):
49 18 equemene
50 18 equemene
    plt.xlabel("Temperature")
51 18 equemene
    plt.ylabel("Energy")
52 18 equemene
53 18 equemene
    Experience,=plt.plot(T,E,label="Energy")
54 18 equemene
55 18 equemene
    plt.legend()
56 18 equemene
    plt.show()
57 18 equemene
58 18 equemene
if __name__=='__main__':
59 18 equemene
60 18 equemene
    # Set defaults values
61 18 equemene
    # Coupling factor
62 18 equemene
    J=1.
63 18 equemene
    # Magnetic Field
64 18 equemene
    B=0.
65 18 equemene
    # Size of Lattice
66 18 equemene
    Size=256
67 18 equemene
    # Default Temperatures (start, end, step)
68 18 equemene
    Tmin=0.1
69 18 equemene
    Tmax=5
70 18 equemene
    Tstep=0.1
71 18 equemene
    # Default Number of Iterations
72 18 equemene
    Iterations=Size*Size
73 18 equemene
74 18 equemene
    # Curves is True to print the curves
75 18 equemene
    Curves=False
76 18 equemene
77 18 equemene
    try:
78 18 equemene
        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep="])
79 18 equemene
    except getopt.GetoptError:
80 18 equemene
        print '%s -j <Coupling Factor> -b <Magnetic Field> -z <Size of Lattice> -i <Iterations> -s <Minimum Temperature> -e <Maximum Temperature> -p <steP Temperature> -c (Print Curves)' % sys.argv[0]
81 18 equemene
        sys.exit(2)
82 18 equemene
83 18 equemene
84 18 equemene
    for opt, arg in opts:
85 18 equemene
        if opt == '-h':
86 18 equemene
            print '%s -j <Coupling Factor> -b <Magnetic Field> -z <Size of Lattice> -i <Iterations> -s <Minimum Temperature> -e <Maximum Temperature> -p <steP Temperature> -c (Print Curves)' % sys.argv[0]
87 18 equemene
            sys.exit()
88 18 equemene
        elif opt == '-c':
89 18 equemene
            Curves=True
90 18 equemene
        elif opt in ("-j", "--coupling"):
91 18 equemene
            J = float(arg)
92 18 equemene
        elif opt in ("-b", "--magneticfield"):
93 18 equemene
            B = float(arg)
94 18 equemene
        elif opt in ("-s", "--tempmin"):
95 18 equemene
            Tmin = float(arg)
96 18 equemene
        elif opt in ("-e", "--tempmax"):
97 18 equemene
            Tmax = arg
98 18 equemene
        elif opt in ("-p", "--tempstep"):
99 18 equemene
            Tstep = numpy.uint32(arg)
100 18 equemene
        elif opt in ("-i", "--iterations"):
101 18 equemene
            Iterations = int(arg)
102 18 equemene
        elif opt in ("-z", "--size"):
103 18 equemene
            Size = int(arg)
104 18 equemene
105 18 equemene
    print "Coupling Factor : %s" % J
106 18 equemene
    print "Magnetic Field :  %s" % B
107 18 equemene
    print "Size of lattice : %s" % Size
108 18 equemene
    print "Iterations : %s" % Iterations
109 18 equemene
    print "Temperature on start : %s" % Tmin
110 18 equemene
    print "Temperature on end : %s" % Tmax
111 18 equemene
    print "Temperature step : %s" % Tstep
112 18 equemene
113 18 equemene
    LAPIMAGE=False
114 18 equemene
115 18 equemene
    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
116 18 equemene
117 18 equemene
    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
118 18 equemene
119 18 equemene
    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
120 18 equemene
121 18 equemene
    E=[]
122 18 equemene
    M=[]
123 18 equemene
124 18 equemene
    for T in Trange:
125 18 equemene
        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
126 18 equemene
        sigma=numpy.copy(sigmaIn)
127 18 equemene
        duration=Metropolis(sigma,J,B,T,Iterations)
128 18 equemene
        E=numpy.append(E,Energy(sigma,J))
129 18 equemene
        M=numpy.append(M,Magnetization(sigma,B))
130 18 equemene
        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
131 18 equemene
132 18 equemene
        print "CPU Time : %f" % (duration)
133 18 equemene
        print "Total Energy at Temperature %f : %f" % (T,E[-1])
134 18 equemene
        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
135 18 equemene
136 18 equemene
    if Curves:
137 18 equemene
        DisplayCurves(Trange,E,M,J,B)
138 18 equemene
139 18 equemene
    # Save output
140 18 equemene
    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))