Statistiques
| Révision :

root / Ising / Cython / Ising2D.py @ 89

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

1
#!/usr/bin/env python
2
#
3
# Ising2D model in serial mode
4
#
5
# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
6

    
7
import sys
8
import numpy
9
from PIL import Image
10
from math import exp
11
from random import random
12
import time
13
import getopt
14
import matplotlib.pyplot as plt
15
import Metropolis
16
from Metropolis import Metropolis
17

    
18
def ImageOutput(sigma,prefix):
19
    Max=sigma.max()
20
    Min=sigma.min()
21
    
22
    # Normalize value as 8bits Integer
23
    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
24
    image = Image.fromarray(SigmaInt)
25
    image.save("%s.jpg" % prefix)
26

    
27
def Magnetization(sigma,M):
28
    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
29

    
30
def Energy(sigma,J):
31
    # Copier et caster 
32
    E=numpy.copy(sigma).astype(numpy.float32)
33
    
34
    # Appel par slice
35
    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
36
                                  E[1:-1,:-2]+E[1:-1,2:])
37
    
38
    # Bien nettoyer la peripherie
39
    E[:,0]=0
40
    E[:,-1]=0
41
    E[0,:]=0
42
    E[-1,:]=0
43
    
44
    Energy=numpy.sum(E)
45

    
46
    return(Energy/(E.shape[0]*E.shape[1]*1.0))
47

    
48
def DisplayCurves(T,E,M,J,B):
49

    
50
    plt.xlabel("Temperature")
51
    plt.ylabel("Energy")
52

    
53
    Experience,=plt.plot(T,E,label="Energy") 
54

    
55
    plt.legend()
56
    plt.show()
57

    
58
if __name__=='__main__':
59

    
60
    # Set defaults values
61
    # Coupling factor
62
    J=1.
63
    # Magnetic Field
64
    B=0.
65
    # Size of Lattice
66
    Size=256
67
    # Default Temperatures (start, end, step)
68
    Tmin=0.1
69
    Tmax=5
70
    Tstep=0.1
71
    # Default Number of Iterations
72
    Iterations=Size*Size
73

    
74
    # Curves is True to print the curves
75
    Curves=False
76

    
77
    try:
78
        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep="])
79
    except getopt.GetoptError:
80
        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
        sys.exit(2)
82
    
83
 
84
    for opt, arg in opts:
85
        if opt == '-h':
86
            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
            sys.exit()
88
        elif opt == '-c':
89
            Curves=True
90
        elif opt in ("-j", "--coupling"):
91
            J = float(arg)
92
        elif opt in ("-b", "--magneticfield"):
93
            B = float(arg)
94
        elif opt in ("-s", "--tempmin"):
95
            Tmin = float(arg)
96
        elif opt in ("-e", "--tempmax"):
97
            Tmax = arg
98
        elif opt in ("-p", "--tempstep"):
99
            Tstep = numpy.uint32(arg)
100
        elif opt in ("-i", "--iterations"):
101
            Iterations = int(arg)
102
        elif opt in ("-z", "--size"):
103
            Size = int(arg)
104
      
105
    print "Coupling Factor : %s" % J
106
    print "Magnetic Field :  %s" % B
107
    print "Size of lattice : %s" % Size
108
    print "Iterations : %s" % Iterations
109
    print "Temperature on start : %s" % Tmin
110
    print "Temperature on end : %s" % Tmax
111
    print "Temperature step : %s" % Tstep
112

    
113
    LAPIMAGE=False
114

    
115
    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
116

    
117
    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
118

    
119
    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
120

    
121
    E=[]
122
    M=[]
123

    
124
    for T in Trange:
125
        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
126
        sigma=numpy.copy(sigmaIn)
127
        duration=Metropolis(sigma,J,B,T,Iterations)
128
        E=numpy.append(E,Energy(sigma,J))
129
        M=numpy.append(M,Magnetization(sigma,B))
130
        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
131
    
132
        print "CPU Time : %f" % (duration)
133
        print "Total Energy at Temperature %f : %f" % (T,E[-1])
134
        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
135

    
136
    if Curves:
137
        DisplayCurves(Trange,E,M,J,B)
138

    
139
    # Save output
140
    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
141