root / Ising / Cython / Ising2D.py
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))
|