Index: Ising/Cython/Metropolis.c
===================================================================
--- Ising/Cython/Metropolis.c	(revision 0)
+++ Ising/Cython/Metropolis.c	(revision 18)
@@ -0,0 +1,2208 @@
+/* Generated by Cython 0.15.1 on Sun May 26 16:17:32 2013 */
+
+#define PY_SSIZE_T_CLEAN
+#include "Python.h"
+#ifndef Py_PYTHON_H
+    #error Python headers needed to compile C extensions, please install development version of Python.
+#else
+
+#include <stddef.h> /* For offsetof */
+#ifndef offsetof
+#define offsetof(type, member) ( (size_t) & ((type*)0) -> member )
+#endif
+
+#if !defined(WIN32) && !defined(MS_WINDOWS)
+  #ifndef __stdcall
+    #define __stdcall
+  #endif
+  #ifndef __cdecl
+    #define __cdecl
+  #endif
+  #ifndef __fastcall
+    #define __fastcall
+  #endif
+#endif
+
+#ifndef DL_IMPORT
+  #define DL_IMPORT(t) t
+#endif
+#ifndef DL_EXPORT
+  #define DL_EXPORT(t) t
+#endif
+
+#ifndef PY_LONG_LONG
+  #define PY_LONG_LONG LONG_LONG
+#endif
+
+#if PY_VERSION_HEX < 0x02040000
+  #define METH_COEXIST 0
+  #define PyDict_CheckExact(op) (Py_TYPE(op) == &PyDict_Type)
+  #define PyDict_Contains(d,o)   PySequence_Contains(d,o)
+#endif
+
+#if PY_VERSION_HEX < 0x02050000
+  typedef int Py_ssize_t;
+  #define PY_SSIZE_T_MAX INT_MAX
+  #define PY_SSIZE_T_MIN INT_MIN
+  #define PY_FORMAT_SIZE_T ""
+  #define PyInt_FromSsize_t(z) PyInt_FromLong(z)
+  #define PyInt_AsSsize_t(o)   __Pyx_PyInt_AsInt(o)
+  #define PyNumber_Index(o)    PyNumber_Int(o)
+  #define PyIndex_Check(o)     PyNumber_Check(o)
+  #define PyErr_WarnEx(category, message, stacklevel) PyErr_Warn(category, message)
+#endif
+
+#if PY_VERSION_HEX < 0x02060000
+  #define Py_REFCNT(ob) (((PyObject*)(ob))->ob_refcnt)
+  #define Py_TYPE(ob)   (((PyObject*)(ob))->ob_type)
+  #define Py_SIZE(ob)   (((PyVarObject*)(ob))->ob_size)
+  #define PyVarObject_HEAD_INIT(type, size) \
+          PyObject_HEAD_INIT(type) size,
+  #define PyType_Modified(t)
+
+  typedef struct {
+     void *buf;
+     PyObject *obj;
+     Py_ssize_t len;
+     Py_ssize_t itemsize;
+     int readonly;
+     int ndim;
+     char *format;
+     Py_ssize_t *shape;
+     Py_ssize_t *strides;
+     Py_ssize_t *suboffsets;
+     void *internal;
+  } Py_buffer;
+
+  #define PyBUF_SIMPLE 0
+  #define PyBUF_WRITABLE 0x0001
+  #define PyBUF_FORMAT 0x0004
+  #define PyBUF_ND 0x0008
+  #define PyBUF_STRIDES (0x0010 | PyBUF_ND)
+  #define PyBUF_C_CONTIGUOUS (0x0020 | PyBUF_STRIDES)
+  #define PyBUF_F_CONTIGUOUS (0x0040 | PyBUF_STRIDES)
+  #define PyBUF_ANY_CONTIGUOUS (0x0080 | PyBUF_STRIDES)
+  #define PyBUF_INDIRECT (0x0100 | PyBUF_STRIDES)
+
+#endif
+
+#if PY_MAJOR_VERSION < 3
+  #define __Pyx_BUILTIN_MODULE_NAME "__builtin__"
+#else
+  #define __Pyx_BUILTIN_MODULE_NAME "builtins"
+#endif
+
+#if PY_MAJOR_VERSION >= 3
+  #define Py_TPFLAGS_CHECKTYPES 0
+  #define Py_TPFLAGS_HAVE_INDEX 0
+#endif
+
+#if (PY_VERSION_HEX < 0x02060000) || (PY_MAJOR_VERSION >= 3)
+  #define Py_TPFLAGS_HAVE_NEWBUFFER 0
+#endif
+
+#if PY_MAJOR_VERSION >= 3
+  #define PyBaseString_Type            PyUnicode_Type
+  #define PyStringObject               PyUnicodeObject
+  #define PyString_Type                PyUnicode_Type
+  #define PyString_Check               PyUnicode_Check
+  #define PyString_CheckExact          PyUnicode_CheckExact
+#endif
+
+#if PY_VERSION_HEX < 0x02060000
+  #define PyBytesObject                PyStringObject
+  #define PyBytes_Type                 PyString_Type
+  #define PyBytes_Check                PyString_Check
+  #define PyBytes_CheckExact           PyString_CheckExact
+  #define PyBytes_FromString           PyString_FromString
+  #define PyBytes_FromStringAndSize    PyString_FromStringAndSize
+  #define PyBytes_FromFormat           PyString_FromFormat
+  #define PyBytes_DecodeEscape         PyString_DecodeEscape
+  #define PyBytes_AsString             PyString_AsString
+  #define PyBytes_AsStringAndSize      PyString_AsStringAndSize
+  #define PyBytes_Size                 PyString_Size
+  #define PyBytes_AS_STRING            PyString_AS_STRING
+  #define PyBytes_GET_SIZE             PyString_GET_SIZE
+  #define PyBytes_Repr                 PyString_Repr
+  #define PyBytes_Concat               PyString_Concat
+  #define PyBytes_ConcatAndDel         PyString_ConcatAndDel
+#endif
+
+#if PY_VERSION_HEX < 0x02060000
+  #define PySet_Check(obj)             PyObject_TypeCheck(obj, &PySet_Type)
+  #define PyFrozenSet_Check(obj)       PyObject_TypeCheck(obj, &PyFrozenSet_Type)
+#endif
+#ifndef PySet_CheckExact
+  #define PySet_CheckExact(obj)        (Py_TYPE(obj) == &PySet_Type)
+#endif
+
+#define __Pyx_TypeCheck(obj, type) PyObject_TypeCheck(obj, (PyTypeObject *)type)
+
+#if PY_MAJOR_VERSION >= 3
+  #define PyIntObject                  PyLongObject
+  #define PyInt_Type                   PyLong_Type
+  #define PyInt_Check(op)              PyLong_Check(op)
+  #define PyInt_CheckExact(op)         PyLong_CheckExact(op)
+  #define PyInt_FromString             PyLong_FromString
+  #define PyInt_FromUnicode            PyLong_FromUnicode
+  #define PyInt_FromLong               PyLong_FromLong
+  #define PyInt_FromSize_t             PyLong_FromSize_t
+  #define PyInt_FromSsize_t            PyLong_FromSsize_t
+  #define PyInt_AsLong                 PyLong_AsLong
+  #define PyInt_AS_LONG                PyLong_AS_LONG
+  #define PyInt_AsSsize_t              PyLong_AsSsize_t
+  #define PyInt_AsUnsignedLongMask     PyLong_AsUnsignedLongMask
+  #define PyInt_AsUnsignedLongLongMask PyLong_AsUnsignedLongLongMask
+#endif
+
+#if PY_MAJOR_VERSION >= 3
+  #define PyBoolObject                 PyLongObject
+#endif
+
+#if PY_VERSION_HEX < 0x03020000
+  typedef long Py_hash_t;
+  #define __Pyx_PyInt_FromHash_t PyInt_FromLong
+  #define __Pyx_PyInt_AsHash_t   PyInt_AsLong
+#else
+  #define __Pyx_PyInt_FromHash_t PyInt_FromSsize_t
+  #define __Pyx_PyInt_AsHash_t   PyInt_AsSsize_t
+#endif
+
+
+#if PY_MAJOR_VERSION >= 3
+  #define __Pyx_PyNumber_Divide(x,y)         PyNumber_TrueDivide(x,y)
+  #define __Pyx_PyNumber_InPlaceDivide(x,y)  PyNumber_InPlaceTrueDivide(x,y)
+#else
+  #define __Pyx_PyNumber_Divide(x,y)         PyNumber_Divide(x,y)
+  #define __Pyx_PyNumber_InPlaceDivide(x,y)  PyNumber_InPlaceDivide(x,y)
+#endif
+
+#if (PY_MAJOR_VERSION < 3) || (PY_VERSION_HEX >= 0x03010300)
+  #define __Pyx_PySequence_GetSlice(obj, a, b) PySequence_GetSlice(obj, a, b)
+  #define __Pyx_PySequence_SetSlice(obj, a, b, value) PySequence_SetSlice(obj, a, b, value)
+  #define __Pyx_PySequence_DelSlice(obj, a, b) PySequence_DelSlice(obj, a, b)
+#else
+  #define __Pyx_PySequence_GetSlice(obj, a, b) (unlikely(!(obj)) ? \
+        (PyErr_SetString(PyExc_SystemError, "null argument to internal routine"), (PyObject*)0) : \
+        (likely((obj)->ob_type->tp_as_mapping) ? (PySequence_GetSlice(obj, a, b)) : \
+            (PyErr_Format(PyExc_TypeError, "'%.200s' object is unsliceable", (obj)->ob_type->tp_name), (PyObject*)0)))
+  #define __Pyx_PySequence_SetSlice(obj, a, b, value) (unlikely(!(obj)) ? \
+        (PyErr_SetString(PyExc_SystemError, "null argument to internal routine"), -1) : \
+        (likely((obj)->ob_type->tp_as_mapping) ? (PySequence_SetSlice(obj, a, b, value)) : \
+            (PyErr_Format(PyExc_TypeError, "'%.200s' object doesn't support slice assignment", (obj)->ob_type->tp_name), -1)))
+  #define __Pyx_PySequence_DelSlice(obj, a, b) (unlikely(!(obj)) ? \
+        (PyErr_SetString(PyExc_SystemError, "null argument to internal routine"), -1) : \
+        (likely((obj)->ob_type->tp_as_mapping) ? (PySequence_DelSlice(obj, a, b)) : \
+            (PyErr_Format(PyExc_TypeError, "'%.200s' object doesn't support slice deletion", (obj)->ob_type->tp_name), -1)))
+#endif
+
+#if PY_MAJOR_VERSION >= 3
+  #define PyMethod_New(func, self, klass) ((self) ? PyMethod_New(func, self) : PyInstanceMethod_New(func))
+#endif
+
+#if PY_VERSION_HEX < 0x02050000
+  #define __Pyx_GetAttrString(o,n)   PyObject_GetAttrString((o),((char *)(n)))
+  #define __Pyx_SetAttrString(o,n,a) PyObject_SetAttrString((o),((char *)(n)),(a))
+  #define __Pyx_DelAttrString(o,n)   PyObject_DelAttrString((o),((char *)(n)))
+#else
+  #define __Pyx_GetAttrString(o,n)   PyObject_GetAttrString((o),(n))
+  #define __Pyx_SetAttrString(o,n,a) PyObject_SetAttrString((o),(n),(a))
+  #define __Pyx_DelAttrString(o,n)   PyObject_DelAttrString((o),(n))
+#endif
+
+#if PY_VERSION_HEX < 0x02050000
+  #define __Pyx_NAMESTR(n) ((char *)(n))
+  #define __Pyx_DOCSTR(n)  ((char *)(n))
+#else
+  #define __Pyx_NAMESTR(n) (n)
+  #define __Pyx_DOCSTR(n)  (n)
+#endif
+
+#ifndef __PYX_EXTERN_C
+  #ifdef __cplusplus
+    #define __PYX_EXTERN_C extern "C"
+  #else
+    #define __PYX_EXTERN_C extern
+  #endif
+#endif
+
+#if defined(WIN32) || defined(MS_WINDOWS)
+#define _USE_MATH_DEFINES
+#endif
+#include <math.h>
+#define __PYX_HAVE__Metropolis
+#define __PYX_HAVE_API__Metropolis
+#ifdef _OPENMP
+#include <omp.h>
+#endif /* _OPENMP */
+
+#ifdef PYREX_WITHOUT_ASSERTIONS
+#define CYTHON_WITHOUT_ASSERTIONS
+#endif
+
+
+/* inline attribute */
+#ifndef CYTHON_INLINE
+  #if defined(__GNUC__)
+    #define CYTHON_INLINE __inline__
+  #elif defined(_MSC_VER)
+    #define CYTHON_INLINE __inline
+  #elif defined (__STDC_VERSION__) && __STDC_VERSION__ >= 199901L
+    #define CYTHON_INLINE inline
+  #else
+    #define CYTHON_INLINE
+  #endif
+#endif
+
+/* unused attribute */
+#ifndef CYTHON_UNUSED
+# if defined(__GNUC__)
+#   if !(defined(__cplusplus)) || (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4))
+#     define CYTHON_UNUSED __attribute__ ((__unused__))
+#   else
+#     define CYTHON_UNUSED
+#   endif
+# elif defined(__ICC) || defined(__INTEL_COMPILER)
+#   define CYTHON_UNUSED __attribute__ ((__unused__))
+# else
+#   define CYTHON_UNUSED
+# endif
+#endif
+
+typedef struct {PyObject **p; char *s; const long n; const char* encoding; const char is_unicode; const char is_str; const char intern; } __Pyx_StringTabEntry; /*proto*/
+
+
+/* Type Conversion Predeclarations */
+
+#define __Pyx_PyBytes_FromUString(s) PyBytes_FromString((char*)s)
+#define __Pyx_PyBytes_AsUString(s)   ((unsigned char*) PyBytes_AsString(s))
+
+#define __Pyx_Owned_Py_None(b) (Py_INCREF(Py_None), Py_None)
+#define __Pyx_PyBool_FromLong(b) ((b) ? (Py_INCREF(Py_True), Py_True) : (Py_INCREF(Py_False), Py_False))
+static CYTHON_INLINE int __Pyx_PyObject_IsTrue(PyObject*);
+static CYTHON_INLINE PyObject* __Pyx_PyNumber_Int(PyObject* x);
+
+static CYTHON_INLINE Py_ssize_t __Pyx_PyIndex_AsSsize_t(PyObject*);
+static CYTHON_INLINE PyObject * __Pyx_PyInt_FromSize_t(size_t);
+static CYTHON_INLINE size_t __Pyx_PyInt_AsSize_t(PyObject*);
+
+#define __pyx_PyFloat_AsDouble(x) (PyFloat_CheckExact(x) ? PyFloat_AS_DOUBLE(x) : PyFloat_AsDouble(x))
+
+
+#ifdef __GNUC__
+  /* Test for GCC > 2.95 */
+  #if __GNUC__ > 2 || (__GNUC__ == 2 && (__GNUC_MINOR__ > 95))
+    #define likely(x)   __builtin_expect(!!(x), 1)
+    #define unlikely(x) __builtin_expect(!!(x), 0)
+  #else /* __GNUC__ > 2 ... */
+    #define likely(x)   (x)
+    #define unlikely(x) (x)
+  #endif /* __GNUC__ > 2 ... */
+#else /* __GNUC__ */
+  #define likely(x)   (x)
+  #define unlikely(x) (x)
+#endif /* __GNUC__ */
+    
+static PyObject *__pyx_m;
+static PyObject *__pyx_b;
+static PyObject *__pyx_empty_tuple;
+static PyObject *__pyx_empty_bytes;
+static int __pyx_lineno;
+static int __pyx_clineno = 0;
+static const char * __pyx_cfilenm= __FILE__;
+static const char *__pyx_filename;
+
+
+static const char *__pyx_f[] = {
+  "Metropolis.pyx",
+};
+
+/*--- Type declarations ---*/
+
+#ifndef CYTHON_REFNANNY
+  #define CYTHON_REFNANNY 0
+#endif
+
+#if CYTHON_REFNANNY
+  typedef struct {
+    void (*INCREF)(void*, PyObject*, int);
+    void (*DECREF)(void*, PyObject*, int);
+    void (*GOTREF)(void*, PyObject*, int);
+    void (*GIVEREF)(void*, PyObject*, int);
+    void* (*SetupContext)(const char*, int, const char*);
+    void (*FinishContext)(void**);
+  } __Pyx_RefNannyAPIStruct;
+  static __Pyx_RefNannyAPIStruct *__Pyx_RefNanny = NULL;
+  static __Pyx_RefNannyAPIStruct *__Pyx_RefNannyImportAPI(const char *modname); /*proto*/
+  #define __Pyx_RefNannyDeclarations void *__pyx_refnanny = NULL;
+  #define __Pyx_RefNannySetupContext(name)           __pyx_refnanny = __Pyx_RefNanny->SetupContext((name), __LINE__, __FILE__)
+  #define __Pyx_RefNannyFinishContext()           __Pyx_RefNanny->FinishContext(&__pyx_refnanny)
+  #define __Pyx_INCREF(r)  __Pyx_RefNanny->INCREF(__pyx_refnanny, (PyObject *)(r), __LINE__)
+  #define __Pyx_DECREF(r)  __Pyx_RefNanny->DECREF(__pyx_refnanny, (PyObject *)(r), __LINE__)
+  #define __Pyx_GOTREF(r)  __Pyx_RefNanny->GOTREF(__pyx_refnanny, (PyObject *)(r), __LINE__)
+  #define __Pyx_GIVEREF(r) __Pyx_RefNanny->GIVEREF(__pyx_refnanny, (PyObject *)(r), __LINE__)
+  #define __Pyx_XINCREF(r)  do { if((r) != NULL) {__Pyx_INCREF(r); }} while(0)
+  #define __Pyx_XDECREF(r)  do { if((r) != NULL) {__Pyx_DECREF(r); }} while(0)
+  #define __Pyx_XGOTREF(r)  do { if((r) != NULL) {__Pyx_GOTREF(r); }} while(0)
+  #define __Pyx_XGIVEREF(r) do { if((r) != NULL) {__Pyx_GIVEREF(r);}} while(0)
+#else
+  #define __Pyx_RefNannyDeclarations
+  #define __Pyx_RefNannySetupContext(name)
+  #define __Pyx_RefNannyFinishContext()
+  #define __Pyx_INCREF(r) Py_INCREF(r)
+  #define __Pyx_DECREF(r) Py_DECREF(r)
+  #define __Pyx_GOTREF(r)
+  #define __Pyx_GIVEREF(r)
+  #define __Pyx_XINCREF(r) Py_XINCREF(r)
+  #define __Pyx_XDECREF(r) Py_XDECREF(r)
+  #define __Pyx_XGOTREF(r)
+  #define __Pyx_XGIVEREF(r)
+#endif /* CYTHON_REFNANNY */
+
+static PyObject *__Pyx_GetName(PyObject *dict, PyObject *name); /*proto*/
+
+static void __Pyx_RaiseArgtupleInvalid(const char* func_name, int exact,
+    Py_ssize_t num_min, Py_ssize_t num_max, Py_ssize_t num_found); /*proto*/
+
+static void __Pyx_RaiseDoubleKeywordsError(
+    const char* func_name, PyObject* kw_name); /*proto*/
+
+static int __Pyx_ParseOptionalKeywords(PyObject *kwds, PyObject **argnames[],     PyObject *kwds2, PyObject *values[], Py_ssize_t num_pos_args,     const char* function_name); /*proto*/
+
+static CYTHON_INLINE void __Pyx_RaiseNeedMoreValuesError(Py_ssize_t index);
+
+static CYTHON_INLINE void __Pyx_RaiseTooManyValuesError(Py_ssize_t expected);
+
+static int __Pyx_IternextUnpackEndCheck(PyObject *retval, Py_ssize_t expected); /*proto*/
+
+static PyObject *__Pyx_Import(PyObject *name, PyObject *from_list, long level); /*proto*/
+
+static CYTHON_INLINE unsigned char __Pyx_PyInt_AsUnsignedChar(PyObject *);
+
+static CYTHON_INLINE unsigned short __Pyx_PyInt_AsUnsignedShort(PyObject *);
+
+static CYTHON_INLINE unsigned int __Pyx_PyInt_AsUnsignedInt(PyObject *);
+
+static CYTHON_INLINE char __Pyx_PyInt_AsChar(PyObject *);
+
+static CYTHON_INLINE short __Pyx_PyInt_AsShort(PyObject *);
+
+static CYTHON_INLINE int __Pyx_PyInt_AsInt(PyObject *);
+
+static CYTHON_INLINE signed char __Pyx_PyInt_AsSignedChar(PyObject *);
+
+static CYTHON_INLINE signed short __Pyx_PyInt_AsSignedShort(PyObject *);
+
+static CYTHON_INLINE signed int __Pyx_PyInt_AsSignedInt(PyObject *);
+
+static CYTHON_INLINE int __Pyx_PyInt_AsLongDouble(PyObject *);
+
+static CYTHON_INLINE unsigned long __Pyx_PyInt_AsUnsignedLong(PyObject *);
+
+static CYTHON_INLINE unsigned PY_LONG_LONG __Pyx_PyInt_AsUnsignedLongLong(PyObject *);
+
+static CYTHON_INLINE long __Pyx_PyInt_AsLong(PyObject *);
+
+static CYTHON_INLINE PY_LONG_LONG __Pyx_PyInt_AsLongLong(PyObject *);
+
+static CYTHON_INLINE signed long __Pyx_PyInt_AsSignedLong(PyObject *);
+
+static CYTHON_INLINE signed PY_LONG_LONG __Pyx_PyInt_AsSignedLongLong(PyObject *);
+
+static int __Pyx_check_binary_version(void);
+
+static void __Pyx_AddTraceback(const char *funcname, int __pyx_clineno,
+                               int __pyx_lineno, const char *__pyx_filename); /*proto*/
+
+static int __Pyx_InitStrings(__Pyx_StringTabEntry *t); /*proto*/
+
+/* Module declarations from 'Metropolis' */
+#define __Pyx_MODULE_NAME "Metropolis"
+int __pyx_module_is_main_Metropolis = 0;
+
+/* Implementation of 'Metropolis' */
+static PyObject *__pyx_builtin_xrange;
+static char __pyx_k_1[] = "matplotlib.pyplot";
+static char __pyx_k_2[] = "*";
+static char __pyx_k__B[] = "B";
+static char __pyx_k__J[] = "J";
+static char __pyx_k__T[] = "T";
+static char __pyx_k__PIL[] = "PIL";
+static char __pyx_k__exp[] = "exp";
+static char __pyx_k__plt[] = "plt";
+static char __pyx_k__sys[] = "sys";
+static char __pyx_k__math[] = "math";
+static char __pyx_k__time[] = "time";
+static char __pyx_k__Image[] = "Image";
+static char __pyx_k__numpy[] = "numpy";
+static char __pyx_k__range[] = "range";
+static char __pyx_k__shape[] = "shape";
+static char __pyx_k__sigma[] = "sigma";
+static char __pyx_k__getopt[] = "getopt";
+static char __pyx_k__random[] = "random";
+static char __pyx_k__xrange[] = "xrange";
+static char __pyx_k__randint[] = "randint";
+static char __pyx_k____main__[] = "__main__";
+static char __pyx_k____test__[] = "__test__";
+static char __pyx_k__Metropolis[] = "Metropolis";
+static char __pyx_k__iterations[] = "iterations";
+static PyObject *__pyx_n_s_1;
+static PyObject *__pyx_n_s_2;
+static PyObject *__pyx_n_s__B;
+static PyObject *__pyx_n_s__Image;
+static PyObject *__pyx_n_s__J;
+static PyObject *__pyx_n_s__Metropolis;
+static PyObject *__pyx_n_s__PIL;
+static PyObject *__pyx_n_s__T;
+static PyObject *__pyx_n_s____main__;
+static PyObject *__pyx_n_s____test__;
+static PyObject *__pyx_n_s__exp;
+static PyObject *__pyx_n_s__getopt;
+static PyObject *__pyx_n_s__iterations;
+static PyObject *__pyx_n_s__math;
+static PyObject *__pyx_n_s__numpy;
+static PyObject *__pyx_n_s__plt;
+static PyObject *__pyx_n_s__randint;
+static PyObject *__pyx_n_s__random;
+static PyObject *__pyx_n_s__range;
+static PyObject *__pyx_n_s__shape;
+static PyObject *__pyx_n_s__sigma;
+static PyObject *__pyx_n_s__sys;
+static PyObject *__pyx_n_s__time;
+static PyObject *__pyx_n_s__xrange;
+static PyObject *__pyx_int_0;
+static PyObject *__pyx_int_1;
+static PyObject *__pyx_int_2;
+
+/* "Metropolis.pyx":17
+ * import Metropolis
+ * 
+ * def Metropolis(sigma,J,B,T,iterations):             # <<<<<<<<<<<<<<
+ *     start=time.time()
+ * 
+ */
+
+static PyObject *__pyx_pf_10Metropolis_Metropolis(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
+static PyMethodDef __pyx_mdef_10Metropolis_Metropolis = {__Pyx_NAMESTR("Metropolis"), (PyCFunction)__pyx_pf_10Metropolis_Metropolis, METH_VARARGS|METH_KEYWORDS, __Pyx_DOCSTR(0)};
+static PyObject *__pyx_pf_10Metropolis_Metropolis(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
+  PyObject *__pyx_v_sigma = 0;
+  PyObject *__pyx_v_J = 0;
+  PyObject *__pyx_v_B = 0;
+  PyObject *__pyx_v_T = 0;
+  PyObject *__pyx_v_iterations = 0;
+  PyObject *__pyx_v_start = NULL;
+  PyObject *__pyx_v_SizeX = NULL;
+  PyObject *__pyx_v_SizeY = NULL;
+  PyObject *__pyx_v_p = NULL;
+  PyObject *__pyx_v_X = NULL;
+  PyObject *__pyx_v_Y = NULL;
+  PyObject *__pyx_v_DeltaE = NULL;
+  PyObject *__pyx_v_duration = NULL;
+  PyObject *__pyx_r = NULL;
+  __Pyx_RefNannyDeclarations
+  PyObject *__pyx_t_1 = NULL;
+  PyObject *__pyx_t_2 = NULL;
+  PyObject *__pyx_t_3 = NULL;
+  PyObject *__pyx_t_4 = NULL;
+  PyObject *(*__pyx_t_5)(PyObject *);
+  Py_ssize_t __pyx_t_6;
+  PyObject *(*__pyx_t_7)(PyObject *);
+  PyObject *__pyx_t_8 = NULL;
+  int __pyx_t_9;
+  int __pyx_t_10;
+  int __pyx_t_11;
+  int __pyx_lineno = 0;
+  const char *__pyx_filename = NULL;
+  int __pyx_clineno = 0;
+  static PyObject **__pyx_pyargnames[] = {&__pyx_n_s__sigma,&__pyx_n_s__J,&__pyx_n_s__B,&__pyx_n_s__T,&__pyx_n_s__iterations,0};
+  __Pyx_RefNannySetupContext("Metropolis");
+  __pyx_self = __pyx_self;
+  {
+    PyObject* values[5] = {0,0,0,0,0};
+    if (unlikely(__pyx_kwds)) {
+      Py_ssize_t kw_args;
+      switch (PyTuple_GET_SIZE(__pyx_args)) {
+        case  5: values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
+        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
+        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
+        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
+        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
+        case  0: break;
+        default: goto __pyx_L5_argtuple_error;
+      }
+      kw_args = PyDict_Size(__pyx_kwds);
+      switch (PyTuple_GET_SIZE(__pyx_args)) {
+        case  0:
+        values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__sigma);
+        if (likely(values[0])) kw_args--;
+        else goto __pyx_L5_argtuple_error;
+        case  1:
+        values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__J);
+        if (likely(values[1])) kw_args--;
+        else {
+          __Pyx_RaiseArgtupleInvalid("Metropolis", 1, 5, 5, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+        }
+        case  2:
+        values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__B);
+        if (likely(values[2])) kw_args--;
+        else {
+          __Pyx_RaiseArgtupleInvalid("Metropolis", 1, 5, 5, 2); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+        }
+        case  3:
+        values[3] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__T);
+        if (likely(values[3])) kw_args--;
+        else {
+          __Pyx_RaiseArgtupleInvalid("Metropolis", 1, 5, 5, 3); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+        }
+        case  4:
+        values[4] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__iterations);
+        if (likely(values[4])) kw_args--;
+        else {
+          __Pyx_RaiseArgtupleInvalid("Metropolis", 1, 5, 5, 4); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+        }
+      }
+      if (unlikely(kw_args > 0)) {
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "Metropolis") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+      }
+    } else if (PyTuple_GET_SIZE(__pyx_args) != 5) {
+      goto __pyx_L5_argtuple_error;
+    } else {
+      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
+      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
+      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
+      values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
+      values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
+    }
+    __pyx_v_sigma = values[0];
+    __pyx_v_J = values[1];
+    __pyx_v_B = values[2];
+    __pyx_v_T = values[3];
+    __pyx_v_iterations = values[4];
+  }
+  goto __pyx_L4_argument_unpacking_done;
+  __pyx_L5_argtuple_error:;
+  __Pyx_RaiseArgtupleInvalid("Metropolis", 1, 5, 5, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+  __pyx_L3_error:;
+  __Pyx_AddTraceback("Metropolis.Metropolis", __pyx_clineno, __pyx_lineno, __pyx_filename);
+  __Pyx_RefNannyFinishContext();
+  return NULL;
+  __pyx_L4_argument_unpacking_done:;
+
+  /* "Metropolis.pyx":18
+ * 
+ * def Metropolis(sigma,J,B,T,iterations):
+ *     start=time.time()             # <<<<<<<<<<<<<<
+ * 
+ *     SizeX,SizeY=sigma.shape
+ */
+  __pyx_t_1 = __Pyx_GetName(__pyx_m, __pyx_n_s__time); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 18; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__time); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 18; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_2);
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+  __pyx_t_1 = PyObject_Call(__pyx_t_2, ((PyObject *)__pyx_empty_tuple), NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 18; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+  __pyx_v_start = __pyx_t_1;
+  __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":20
+ *     start=time.time()
+ * 
+ *     SizeX,SizeY=sigma.shape             # <<<<<<<<<<<<<<
+ * 
+ *     for p in xrange(0,iterations):
+ */
+  __pyx_t_1 = PyObject_GetAttr(__pyx_v_sigma, __pyx_n_s__shape); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 20; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  if ((likely(PyTuple_CheckExact(__pyx_t_1))) || (PyList_CheckExact(__pyx_t_1))) {
+    PyObject* sequence = __pyx_t_1;
+    if (likely(PyTuple_CheckExact(sequence))) {
+      if (unlikely(PyTuple_GET_SIZE(sequence) != 2)) {
+        if (PyTuple_GET_SIZE(sequence) > 2) __Pyx_RaiseTooManyValuesError(2);
+        else __Pyx_RaiseNeedMoreValuesError(PyTuple_GET_SIZE(sequence));
+        {__pyx_filename = __pyx_f[0]; __pyx_lineno = 20; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      }
+      __pyx_t_2 = PyTuple_GET_ITEM(sequence, 0); 
+      __pyx_t_3 = PyTuple_GET_ITEM(sequence, 1); 
+    } else {
+      if (unlikely(PyList_GET_SIZE(sequence) != 2)) {
+        if (PyList_GET_SIZE(sequence) > 2) __Pyx_RaiseTooManyValuesError(2);
+        else __Pyx_RaiseNeedMoreValuesError(PyList_GET_SIZE(sequence));
+        {__pyx_filename = __pyx_f[0]; __pyx_lineno = 20; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      }
+      __pyx_t_2 = PyList_GET_ITEM(sequence, 0); 
+      __pyx_t_3 = PyList_GET_ITEM(sequence, 1); 
+    }
+    __Pyx_INCREF(__pyx_t_2);
+    __Pyx_INCREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+  } else {
+    Py_ssize_t index = -1;
+    __pyx_t_4 = PyObject_GetIter(__pyx_t_1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 20; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+    __pyx_t_5 = Py_TYPE(__pyx_t_4)->tp_iternext;
+    index = 0; __pyx_t_2 = __pyx_t_5(__pyx_t_4); if (unlikely(!__pyx_t_2)) goto __pyx_L6_unpacking_failed;
+    __Pyx_GOTREF(__pyx_t_2);
+    index = 1; __pyx_t_3 = __pyx_t_5(__pyx_t_4); if (unlikely(!__pyx_t_3)) goto __pyx_L6_unpacking_failed;
+    __Pyx_GOTREF(__pyx_t_3);
+    if (__Pyx_IternextUnpackEndCheck(__pyx_t_5(__pyx_t_4), 2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 20; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    goto __pyx_L7_unpacking_done;
+    __pyx_L6_unpacking_failed:;
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    if (PyErr_Occurred() && PyErr_ExceptionMatches(PyExc_StopIteration)) PyErr_Clear();
+    if (!PyErr_Occurred()) __Pyx_RaiseNeedMoreValuesError(index);
+    {__pyx_filename = __pyx_f[0]; __pyx_lineno = 20; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __pyx_L7_unpacking_done:;
+  }
+  __pyx_v_SizeX = __pyx_t_2;
+  __pyx_t_2 = 0;
+  __pyx_v_SizeY = __pyx_t_3;
+  __pyx_t_3 = 0;
+
+  /* "Metropolis.pyx":22
+ *     SizeX,SizeY=sigma.shape
+ * 
+ *     for p in xrange(0,iterations):             # <<<<<<<<<<<<<<
+ *         # Random access coordonate
+ *         X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+ */
+  __pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(((PyObject *)__pyx_t_1));
+  __Pyx_INCREF(__pyx_int_0);
+  PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_int_0);
+  __Pyx_GIVEREF(__pyx_int_0);
+  __Pyx_INCREF(__pyx_v_iterations);
+  PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_v_iterations);
+  __Pyx_GIVEREF(__pyx_v_iterations);
+  __pyx_t_3 = PyObject_Call(__pyx_builtin_xrange, ((PyObject *)__pyx_t_1), NULL); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_3);
+  __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
+  if (PyList_CheckExact(__pyx_t_3) || PyTuple_CheckExact(__pyx_t_3)) {
+    __pyx_t_1 = __pyx_t_3; __Pyx_INCREF(__pyx_t_1); __pyx_t_6 = 0;
+    __pyx_t_7 = NULL;
+  } else {
+    __pyx_t_6 = -1; __pyx_t_1 = PyObject_GetIter(__pyx_t_3); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_1);
+    __pyx_t_7 = Py_TYPE(__pyx_t_1)->tp_iternext;
+  }
+  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+  for (;;) {
+    if (PyList_CheckExact(__pyx_t_1)) {
+      if (__pyx_t_6 >= PyList_GET_SIZE(__pyx_t_1)) break;
+      __pyx_t_3 = PyList_GET_ITEM(__pyx_t_1, __pyx_t_6); __Pyx_INCREF(__pyx_t_3); __pyx_t_6++;
+    } else if (PyTuple_CheckExact(__pyx_t_1)) {
+      if (__pyx_t_6 >= PyTuple_GET_SIZE(__pyx_t_1)) break;
+      __pyx_t_3 = PyTuple_GET_ITEM(__pyx_t_1, __pyx_t_6); __Pyx_INCREF(__pyx_t_3); __pyx_t_6++;
+    } else {
+      __pyx_t_3 = __pyx_t_7(__pyx_t_1);
+      if (unlikely(!__pyx_t_3)) {
+        if (PyErr_Occurred()) {
+          if (likely(PyErr_ExceptionMatches(PyExc_StopIteration))) PyErr_Clear();
+          else {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+        }
+        break;
+      }
+      __Pyx_GOTREF(__pyx_t_3);
+    }
+    __Pyx_XDECREF(__pyx_v_p);
+    __pyx_v_p = __pyx_t_3;
+    __pyx_t_3 = 0;
+
+    /* "Metropolis.pyx":24
+ *     for p in xrange(0,iterations):
+ *         # Random access coordonate
+ *         X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)             # <<<<<<<<<<<<<<
+ * 
+ *         DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+ */
+    __pyx_t_3 = __Pyx_GetName(__pyx_m, __pyx_n_s__numpy); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __pyx_t_2 = PyObject_GetAttr(__pyx_t_3, __pyx_n_s__random); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_2);
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __pyx_t_3 = PyObject_GetAttr(__pyx_t_2, __pyx_n_s__randint); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+    __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(((PyObject *)__pyx_t_2));
+    __Pyx_INCREF(__pyx_v_SizeX);
+    PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_v_SizeX);
+    __Pyx_GIVEREF(__pyx_v_SizeX);
+    __pyx_t_4 = PyObject_Call(__pyx_t_3, ((PyObject *)__pyx_t_2), NULL); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+    __pyx_t_2 = __Pyx_GetName(__pyx_m, __pyx_n_s__numpy); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_2);
+    __pyx_t_3 = PyObject_GetAttr(__pyx_t_2, __pyx_n_s__random); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+    __pyx_t_2 = PyObject_GetAttr(__pyx_t_3, __pyx_n_s__randint); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_2);
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(((PyObject *)__pyx_t_3));
+    __Pyx_INCREF(__pyx_v_SizeY);
+    PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_v_SizeY);
+    __Pyx_GIVEREF(__pyx_v_SizeY);
+    __pyx_t_8 = PyObject_Call(__pyx_t_2, ((PyObject *)__pyx_t_3), NULL); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_8);
+    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+    __Pyx_DECREF(((PyObject *)__pyx_t_3)); __pyx_t_3 = 0;
+    __Pyx_XDECREF(__pyx_v_X);
+    __pyx_v_X = __pyx_t_4;
+    __pyx_t_4 = 0;
+    __Pyx_XDECREF(__pyx_v_Y);
+    __pyx_v_Y = __pyx_t_8;
+    __pyx_t_8 = 0;
+
+    /* "Metropolis.pyx":26
+ *         X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+ * 
+ *         DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+             # <<<<<<<<<<<<<<
+ *                                 sigma[X,(Y-1)%SizeY]+
+ *                                 sigma[(X-1)%SizeX,Y]+
+ */
+    __pyx_t_8 = PyTuple_New(2); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(((PyObject *)__pyx_t_8));
+    __Pyx_INCREF(__pyx_v_X);
+    PyTuple_SET_ITEM(__pyx_t_8, 0, __pyx_v_X);
+    __Pyx_GIVEREF(__pyx_v_X);
+    __Pyx_INCREF(__pyx_v_Y);
+    PyTuple_SET_ITEM(__pyx_t_8, 1, __pyx_v_Y);
+    __Pyx_GIVEREF(__pyx_v_Y);
+    __pyx_t_4 = PyObject_GetItem(__pyx_v_sigma, ((PyObject *)__pyx_t_8)); if (!__pyx_t_4) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(((PyObject *)__pyx_t_8)); __pyx_t_8 = 0;
+    __pyx_t_8 = PyNumber_Multiply(__pyx_v_J, __pyx_t_4); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_8);
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+
+    /* "Metropolis.pyx":27
+ * 
+ *         DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+ *                                 sigma[X,(Y-1)%SizeY]+             # <<<<<<<<<<<<<<
+ *                                 sigma[(X-1)%SizeX,Y]+
+ *                                 sigma[(X+1)%SizeX,Y])+B)
+ */
+    __pyx_t_4 = PyNumber_Add(__pyx_v_Y, __pyx_int_1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+
+    /* "Metropolis.pyx":26
+ *         X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+ * 
+ *         DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+             # <<<<<<<<<<<<<<
+ *                                 sigma[X,(Y-1)%SizeY]+
+ *                                 sigma[(X-1)%SizeX,Y]+
+ */
+    __pyx_t_3 = PyNumber_Remainder(__pyx_t_4, __pyx_v_SizeY); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    __pyx_t_4 = PyTuple_New(2); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(((PyObject *)__pyx_t_4));
+    __Pyx_INCREF(__pyx_v_X);
+    PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_v_X);
+    __Pyx_GIVEREF(__pyx_v_X);
+    PyTuple_SET_ITEM(__pyx_t_4, 1, __pyx_t_3);
+    __Pyx_GIVEREF(__pyx_t_3);
+    __pyx_t_3 = 0;
+    __pyx_t_3 = PyObject_GetItem(__pyx_v_sigma, ((PyObject *)__pyx_t_4)); if (!__pyx_t_3) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(((PyObject *)__pyx_t_4)); __pyx_t_4 = 0;
+
+    /* "Metropolis.pyx":27
+ * 
+ *         DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+ *                                 sigma[X,(Y-1)%SizeY]+             # <<<<<<<<<<<<<<
+ *                                 sigma[(X-1)%SizeX,Y]+
+ *                                 sigma[(X+1)%SizeX,Y])+B)
+ */
+    __pyx_t_4 = PyNumber_Subtract(__pyx_v_Y, __pyx_int_1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __pyx_t_2 = PyNumber_Remainder(__pyx_t_4, __pyx_v_SizeY); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_2);
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    __pyx_t_4 = PyTuple_New(2); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(((PyObject *)__pyx_t_4));
+    __Pyx_INCREF(__pyx_v_X);
+    PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_v_X);
+    __Pyx_GIVEREF(__pyx_v_X);
+    PyTuple_SET_ITEM(__pyx_t_4, 1, __pyx_t_2);
+    __Pyx_GIVEREF(__pyx_t_2);
+    __pyx_t_2 = 0;
+    __pyx_t_2 = PyObject_GetItem(__pyx_v_sigma, ((PyObject *)__pyx_t_4)); if (!__pyx_t_2) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_2);
+    __Pyx_DECREF(((PyObject *)__pyx_t_4)); __pyx_t_4 = 0;
+    __pyx_t_4 = PyNumber_Add(__pyx_t_3, __pyx_t_2); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+
+    /* "Metropolis.pyx":28
+ *         DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+ *                                 sigma[X,(Y-1)%SizeY]+
+ *                                 sigma[(X-1)%SizeX,Y]+             # <<<<<<<<<<<<<<
+ *                                 sigma[(X+1)%SizeX,Y])+B)
+ * 
+ */
+    __pyx_t_2 = PyNumber_Subtract(__pyx_v_X, __pyx_int_1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 28; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_2);
+    __pyx_t_3 = PyNumber_Remainder(__pyx_t_2, __pyx_v_SizeX); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 28; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+    __pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 28; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(((PyObject *)__pyx_t_2));
+    PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_3);
+    __Pyx_GIVEREF(__pyx_t_3);
+    __Pyx_INCREF(__pyx_v_Y);
+    PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_v_Y);
+    __Pyx_GIVEREF(__pyx_v_Y);
+    __pyx_t_3 = 0;
+    __pyx_t_3 = PyObject_GetItem(__pyx_v_sigma, ((PyObject *)__pyx_t_2)); if (!__pyx_t_3) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 28; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+    __pyx_t_2 = PyNumber_Add(__pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_2);
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+
+    /* "Metropolis.pyx":29
+ *                                 sigma[X,(Y-1)%SizeY]+
+ *                                 sigma[(X-1)%SizeX,Y]+
+ *                                 sigma[(X+1)%SizeX,Y])+B)             # <<<<<<<<<<<<<<
+ * 
+ *         if DeltaE < 0. or random() < exp(-DeltaE/T):
+ */
+    __pyx_t_3 = PyNumber_Add(__pyx_v_X, __pyx_int_1); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 29; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __pyx_t_4 = PyNumber_Remainder(__pyx_t_3, __pyx_v_SizeX); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 29; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 29; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(((PyObject *)__pyx_t_3));
+    PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_4);
+    __Pyx_GIVEREF(__pyx_t_4);
+    __Pyx_INCREF(__pyx_v_Y);
+    PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_v_Y);
+    __Pyx_GIVEREF(__pyx_v_Y);
+    __pyx_t_4 = 0;
+    __pyx_t_4 = PyObject_GetItem(__pyx_v_sigma, ((PyObject *)__pyx_t_3)); if (!__pyx_t_4) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 29; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(((PyObject *)__pyx_t_3)); __pyx_t_3 = 0;
+    __pyx_t_3 = PyNumber_Add(__pyx_t_2, __pyx_t_4); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 28; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    __pyx_t_4 = PyNumber_Multiply(__pyx_int_2, __pyx_t_3); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __pyx_t_3 = PyNumber_Add(__pyx_t_4, __pyx_v_B); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 29; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    __pyx_t_4 = PyNumber_Multiply(__pyx_t_8, __pyx_t_3); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    __Pyx_XDECREF(__pyx_v_DeltaE);
+    __pyx_v_DeltaE = __pyx_t_4;
+    __pyx_t_4 = 0;
+
+    /* "Metropolis.pyx":31
+ *                                 sigma[(X+1)%SizeX,Y])+B)
+ * 
+ *         if DeltaE < 0. or random() < exp(-DeltaE/T):             # <<<<<<<<<<<<<<
+ *             sigma[X,Y]=-sigma[X,Y]
+ *     duration=time.time()-start
+ */
+    __pyx_t_4 = PyFloat_FromDouble(0.); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_4);
+    __pyx_t_3 = PyObject_RichCompare(__pyx_v_DeltaE, __pyx_t_4, Py_LT); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_GOTREF(__pyx_t_3);
+    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+    __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_t_3); if (unlikely(__pyx_t_9 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+    if (!__pyx_t_9) {
+      __pyx_t_3 = __Pyx_GetName(__pyx_m, __pyx_n_s__random); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_3);
+      __pyx_t_4 = PyObject_Call(__pyx_t_3, ((PyObject *)__pyx_empty_tuple), NULL); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_4);
+      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+      __pyx_t_3 = __Pyx_GetName(__pyx_m, __pyx_n_s__exp); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_3);
+      __pyx_t_8 = PyNumber_Negative(__pyx_v_DeltaE); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_8);
+      __pyx_t_2 = __Pyx_PyNumber_Divide(__pyx_t_8, __pyx_v_T); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_2);
+      __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
+      __pyx_t_8 = PyTuple_New(1); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(((PyObject *)__pyx_t_8));
+      PyTuple_SET_ITEM(__pyx_t_8, 0, __pyx_t_2);
+      __Pyx_GIVEREF(__pyx_t_2);
+      __pyx_t_2 = 0;
+      __pyx_t_2 = PyObject_Call(__pyx_t_3, ((PyObject *)__pyx_t_8), NULL); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_2);
+      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
+      __Pyx_DECREF(((PyObject *)__pyx_t_8)); __pyx_t_8 = 0;
+      __pyx_t_8 = PyObject_RichCompare(__pyx_t_4, __pyx_t_2, Py_LT); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_8);
+      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
+      __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+      __pyx_t_10 = __Pyx_PyObject_IsTrue(__pyx_t_8); if (unlikely(__pyx_t_10 < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
+      __pyx_t_11 = __pyx_t_10;
+    } else {
+      __pyx_t_11 = __pyx_t_9;
+    }
+    if (__pyx_t_11) {
+
+      /* "Metropolis.pyx":32
+ * 
+ *         if DeltaE < 0. or random() < exp(-DeltaE/T):
+ *             sigma[X,Y]=-sigma[X,Y]             # <<<<<<<<<<<<<<
+ *     duration=time.time()-start
+ *     return(duration)
+ */
+      __pyx_t_8 = PyTuple_New(2); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 32; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(((PyObject *)__pyx_t_8));
+      __Pyx_INCREF(__pyx_v_X);
+      PyTuple_SET_ITEM(__pyx_t_8, 0, __pyx_v_X);
+      __Pyx_GIVEREF(__pyx_v_X);
+      __Pyx_INCREF(__pyx_v_Y);
+      PyTuple_SET_ITEM(__pyx_t_8, 1, __pyx_v_Y);
+      __Pyx_GIVEREF(__pyx_v_Y);
+      __pyx_t_2 = PyObject_GetItem(__pyx_v_sigma, ((PyObject *)__pyx_t_8)); if (!__pyx_t_2) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 32; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_2);
+      __Pyx_DECREF(((PyObject *)__pyx_t_8)); __pyx_t_8 = 0;
+      __pyx_t_8 = PyNumber_Negative(__pyx_t_2); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 32; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(__pyx_t_8);
+      __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+      __pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 32; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_GOTREF(((PyObject *)__pyx_t_2));
+      __Pyx_INCREF(__pyx_v_X);
+      PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_v_X);
+      __Pyx_GIVEREF(__pyx_v_X);
+      __Pyx_INCREF(__pyx_v_Y);
+      PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_v_Y);
+      __Pyx_GIVEREF(__pyx_v_Y);
+      if (PyObject_SetItem(__pyx_v_sigma, ((PyObject *)__pyx_t_2), __pyx_t_8) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 32; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+      __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+      __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
+      goto __pyx_L10;
+    }
+    __pyx_L10:;
+  }
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":33
+ *         if DeltaE < 0. or random() < exp(-DeltaE/T):
+ *             sigma[X,Y]=-sigma[X,Y]
+ *     duration=time.time()-start             # <<<<<<<<<<<<<<
+ *     return(duration)
+ */
+  __pyx_t_1 = __Pyx_GetName(__pyx_m, __pyx_n_s__time); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 33; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  __pyx_t_8 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__time); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 33; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_8);
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+  __pyx_t_1 = PyObject_Call(__pyx_t_8, ((PyObject *)__pyx_empty_tuple), NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 33; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
+  __pyx_t_8 = PyNumber_Subtract(__pyx_t_1, __pyx_v_start); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 33; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_8);
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+  __pyx_v_duration = __pyx_t_8;
+  __pyx_t_8 = 0;
+
+  /* "Metropolis.pyx":34
+ *             sigma[X,Y]=-sigma[X,Y]
+ *     duration=time.time()-start
+ *     return(duration)             # <<<<<<<<<<<<<<
+ */
+  __Pyx_XDECREF(__pyx_r);
+  __Pyx_INCREF(__pyx_v_duration);
+  __pyx_r = __pyx_v_duration;
+  goto __pyx_L0;
+
+  __pyx_r = Py_None; __Pyx_INCREF(Py_None);
+  goto __pyx_L0;
+  __pyx_L1_error:;
+  __Pyx_XDECREF(__pyx_t_1);
+  __Pyx_XDECREF(__pyx_t_2);
+  __Pyx_XDECREF(__pyx_t_3);
+  __Pyx_XDECREF(__pyx_t_4);
+  __Pyx_XDECREF(__pyx_t_8);
+  __Pyx_AddTraceback("Metropolis.Metropolis", __pyx_clineno, __pyx_lineno, __pyx_filename);
+  __pyx_r = NULL;
+  __pyx_L0:;
+  __Pyx_XDECREF(__pyx_v_start);
+  __Pyx_XDECREF(__pyx_v_SizeX);
+  __Pyx_XDECREF(__pyx_v_SizeY);
+  __Pyx_XDECREF(__pyx_v_p);
+  __Pyx_XDECREF(__pyx_v_X);
+  __Pyx_XDECREF(__pyx_v_Y);
+  __Pyx_XDECREF(__pyx_v_DeltaE);
+  __Pyx_XDECREF(__pyx_v_duration);
+  __Pyx_XGIVEREF(__pyx_r);
+  __Pyx_RefNannyFinishContext();
+  return __pyx_r;
+}
+
+static PyMethodDef __pyx_methods[] = {
+  {0, 0, 0, 0}
+};
+
+#if PY_MAJOR_VERSION >= 3
+static struct PyModuleDef __pyx_moduledef = {
+    PyModuleDef_HEAD_INIT,
+    __Pyx_NAMESTR("Metropolis"),
+    0, /* m_doc */
+    -1, /* m_size */
+    __pyx_methods /* m_methods */,
+    NULL, /* m_reload */
+    NULL, /* m_traverse */
+    NULL, /* m_clear */
+    NULL /* m_free */
+};
+#endif
+
+static __Pyx_StringTabEntry __pyx_string_tab[] = {
+  {&__pyx_n_s_1, __pyx_k_1, sizeof(__pyx_k_1), 0, 0, 1, 1},
+  {&__pyx_n_s_2, __pyx_k_2, sizeof(__pyx_k_2), 0, 0, 1, 1},
+  {&__pyx_n_s__B, __pyx_k__B, sizeof(__pyx_k__B), 0, 0, 1, 1},
+  {&__pyx_n_s__Image, __pyx_k__Image, sizeof(__pyx_k__Image), 0, 0, 1, 1},
+  {&__pyx_n_s__J, __pyx_k__J, sizeof(__pyx_k__J), 0, 0, 1, 1},
+  {&__pyx_n_s__Metropolis, __pyx_k__Metropolis, sizeof(__pyx_k__Metropolis), 0, 0, 1, 1},
+  {&__pyx_n_s__PIL, __pyx_k__PIL, sizeof(__pyx_k__PIL), 0, 0, 1, 1},
+  {&__pyx_n_s__T, __pyx_k__T, sizeof(__pyx_k__T), 0, 0, 1, 1},
+  {&__pyx_n_s____main__, __pyx_k____main__, sizeof(__pyx_k____main__), 0, 0, 1, 1},
+  {&__pyx_n_s____test__, __pyx_k____test__, sizeof(__pyx_k____test__), 0, 0, 1, 1},
+  {&__pyx_n_s__exp, __pyx_k__exp, sizeof(__pyx_k__exp), 0, 0, 1, 1},
+  {&__pyx_n_s__getopt, __pyx_k__getopt, sizeof(__pyx_k__getopt), 0, 0, 1, 1},
+  {&__pyx_n_s__iterations, __pyx_k__iterations, sizeof(__pyx_k__iterations), 0, 0, 1, 1},
+  {&__pyx_n_s__math, __pyx_k__math, sizeof(__pyx_k__math), 0, 0, 1, 1},
+  {&__pyx_n_s__numpy, __pyx_k__numpy, sizeof(__pyx_k__numpy), 0, 0, 1, 1},
+  {&__pyx_n_s__plt, __pyx_k__plt, sizeof(__pyx_k__plt), 0, 0, 1, 1},
+  {&__pyx_n_s__randint, __pyx_k__randint, sizeof(__pyx_k__randint), 0, 0, 1, 1},
+  {&__pyx_n_s__random, __pyx_k__random, sizeof(__pyx_k__random), 0, 0, 1, 1},
+  {&__pyx_n_s__range, __pyx_k__range, sizeof(__pyx_k__range), 0, 0, 1, 1},
+  {&__pyx_n_s__shape, __pyx_k__shape, sizeof(__pyx_k__shape), 0, 0, 1, 1},
+  {&__pyx_n_s__sigma, __pyx_k__sigma, sizeof(__pyx_k__sigma), 0, 0, 1, 1},
+  {&__pyx_n_s__sys, __pyx_k__sys, sizeof(__pyx_k__sys), 0, 0, 1, 1},
+  {&__pyx_n_s__time, __pyx_k__time, sizeof(__pyx_k__time), 0, 0, 1, 1},
+  {&__pyx_n_s__xrange, __pyx_k__xrange, sizeof(__pyx_k__xrange), 0, 0, 1, 1},
+  {0, 0, 0, 0, 0, 0, 0}
+};
+static int __Pyx_InitCachedBuiltins(void) {
+  #if PY_MAJOR_VERSION >= 3
+  __pyx_builtin_xrange = __Pyx_GetName(__pyx_b, __pyx_n_s__range); if (!__pyx_builtin_xrange) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  #else
+  __pyx_builtin_xrange = __Pyx_GetName(__pyx_b, __pyx_n_s__xrange); if (!__pyx_builtin_xrange) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  #endif
+  return 0;
+  __pyx_L1_error:;
+  return -1;
+}
+
+static int __Pyx_InitCachedConstants(void) {
+  __Pyx_RefNannyDeclarations
+  __Pyx_RefNannySetupContext("__Pyx_InitCachedConstants");
+  __Pyx_RefNannyFinishContext();
+  return 0;
+}
+
+static int __Pyx_InitGlobals(void) {
+  if (__Pyx_InitStrings(__pyx_string_tab) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  __pyx_int_0 = PyInt_FromLong(0); if (unlikely(!__pyx_int_0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  __pyx_int_1 = PyInt_FromLong(1); if (unlikely(!__pyx_int_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  __pyx_int_2 = PyInt_FromLong(2); if (unlikely(!__pyx_int_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  return 0;
+  __pyx_L1_error:;
+  return -1;
+}
+
+#if PY_MAJOR_VERSION < 3
+PyMODINIT_FUNC initMetropolis(void); /*proto*/
+PyMODINIT_FUNC initMetropolis(void)
+#else
+PyMODINIT_FUNC PyInit_Metropolis(void); /*proto*/
+PyMODINIT_FUNC PyInit_Metropolis(void)
+#endif
+{
+  PyObject *__pyx_t_1 = NULL;
+  PyObject *__pyx_t_2 = NULL;
+  __Pyx_RefNannyDeclarations
+  #if CYTHON_REFNANNY
+  __Pyx_RefNanny = __Pyx_RefNannyImportAPI("refnanny");
+  if (!__Pyx_RefNanny) {
+      PyErr_Clear();
+      __Pyx_RefNanny = __Pyx_RefNannyImportAPI("Cython.Runtime.refnanny");
+      if (!__Pyx_RefNanny)
+          Py_FatalError("failed to import 'refnanny' module");
+  }
+  #endif
+  __Pyx_RefNannySetupContext("PyMODINIT_FUNC PyInit_Metropolis(void)");
+  if ( __Pyx_check_binary_version() < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __pyx_empty_tuple = PyTuple_New(0); if (unlikely(!__pyx_empty_tuple)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __pyx_empty_bytes = PyBytes_FromStringAndSize("", 0); if (unlikely(!__pyx_empty_bytes)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  #ifdef __pyx_binding_PyCFunctionType_USED
+  if (__pyx_binding_PyCFunctionType_init() < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  #endif
+  /*--- Library function declarations ---*/
+  /*--- Threads initialization code ---*/
+  #if defined(__PYX_FORCE_INIT_THREADS) && __PYX_FORCE_INIT_THREADS
+  #ifdef WITH_THREAD /* Python build with threading support? */
+  PyEval_InitThreads();
+  #endif
+  #endif
+  /*--- Module creation code ---*/
+  #if PY_MAJOR_VERSION < 3
+  __pyx_m = Py_InitModule4(__Pyx_NAMESTR("Metropolis"), __pyx_methods, 0, 0, PYTHON_API_VERSION);
+  #else
+  __pyx_m = PyModule_Create(&__pyx_moduledef);
+  #endif
+  if (!__pyx_m) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  #if PY_MAJOR_VERSION < 3
+  Py_INCREF(__pyx_m);
+  #endif
+  __pyx_b = PyImport_AddModule(__Pyx_NAMESTR(__Pyx_BUILTIN_MODULE_NAME));
+  if (!__pyx_b) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  if (__Pyx_SetAttrString(__pyx_m, "__builtins__", __pyx_b) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  /*--- Initialize various global constants etc. ---*/
+  if (unlikely(__Pyx_InitGlobals() < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  if (__pyx_module_is_main_Metropolis) {
+    if (__Pyx_SetAttrString(__pyx_m, "__name__", __pyx_n_s____main__) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+  }
+  /*--- Builtin init code ---*/
+  if (unlikely(__Pyx_InitCachedBuiltins() < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  /*--- Constants init code ---*/
+  if (unlikely(__Pyx_InitCachedConstants() < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  /*--- Global init code ---*/
+  /*--- Variable export code ---*/
+  /*--- Function export code ---*/
+  /*--- Type init code ---*/
+  /*--- Type import code ---*/
+  /*--- Variable import code ---*/
+  /*--- Function import code ---*/
+  /*--- Execution code ---*/
+
+  /* "Metropolis.pyx":7
+ * # CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr>
+ * 
+ * import sys             # <<<<<<<<<<<<<<
+ * import numpy
+ * from PIL import Image
+ */
+  __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s__sys), 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__sys, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":8
+ * 
+ * import sys
+ * import numpy             # <<<<<<<<<<<<<<
+ * from PIL import Image
+ * from math import exp
+ */
+  __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s__numpy), 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__numpy, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":9
+ * import sys
+ * import numpy
+ * from PIL import Image             # <<<<<<<<<<<<<<
+ * from math import exp
+ * from random import random
+ */
+  __pyx_t_1 = PyList_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(((PyObject *)__pyx_t_1));
+  __Pyx_INCREF(((PyObject *)__pyx_n_s__Image));
+  PyList_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_n_s__Image));
+  __Pyx_GIVEREF(((PyObject *)__pyx_n_s__Image));
+  __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__PIL), ((PyObject *)__pyx_t_1), -1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_2);
+  __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
+  __pyx_t_1 = PyObject_GetAttr(__pyx_t_2, __pyx_n_s__Image); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__Image, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+
+  /* "Metropolis.pyx":10
+ * import numpy
+ * from PIL import Image
+ * from math import exp             # <<<<<<<<<<<<<<
+ * from random import random
+ * import time
+ */
+  __pyx_t_2 = PyList_New(1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 10; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(((PyObject *)__pyx_t_2));
+  __Pyx_INCREF(((PyObject *)__pyx_n_s__exp));
+  PyList_SET_ITEM(__pyx_t_2, 0, ((PyObject *)__pyx_n_s__exp));
+  __Pyx_GIVEREF(((PyObject *)__pyx_n_s__exp));
+  __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s__math), ((PyObject *)__pyx_t_2), -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 10; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+  __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__exp); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 10; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_2);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__exp, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 10; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":11
+ * from PIL import Image
+ * from math import exp
+ * from random import random             # <<<<<<<<<<<<<<
+ * import time
+ * import getopt
+ */
+  __pyx_t_1 = PyList_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(((PyObject *)__pyx_t_1));
+  __Pyx_INCREF(((PyObject *)__pyx_n_s__random));
+  PyList_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_n_s__random));
+  __Pyx_GIVEREF(((PyObject *)__pyx_n_s__random));
+  __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__random), ((PyObject *)__pyx_t_1), -1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_2);
+  __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
+  __pyx_t_1 = PyObject_GetAttr(__pyx_t_2, __pyx_n_s__random); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__random, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+
+  /* "Metropolis.pyx":12
+ * from math import exp
+ * from random import random
+ * import time             # <<<<<<<<<<<<<<
+ * import getopt
+ * import matplotlib.pyplot as plt
+ */
+  __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__time), 0, -1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_2);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__time, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+
+  /* "Metropolis.pyx":13
+ * from random import random
+ * import time
+ * import getopt             # <<<<<<<<<<<<<<
+ * import matplotlib.pyplot as plt
+ * import Metropolis
+ */
+  __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__getopt), 0, -1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_2);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__getopt, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+
+  /* "Metropolis.pyx":14
+ * import time
+ * import getopt
+ * import matplotlib.pyplot as plt             # <<<<<<<<<<<<<<
+ * import Metropolis
+ * 
+ */
+  __pyx_t_2 = PyList_New(1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(((PyObject *)__pyx_t_2));
+  __Pyx_INCREF(((PyObject *)__pyx_n_s_2));
+  PyList_SET_ITEM(__pyx_t_2, 0, ((PyObject *)__pyx_n_s_2));
+  __Pyx_GIVEREF(((PyObject *)__pyx_n_s_2));
+  __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s_1), ((PyObject *)__pyx_t_2), -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__plt, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":15
+ * import getopt
+ * import matplotlib.pyplot as plt
+ * import Metropolis             # <<<<<<<<<<<<<<
+ * 
+ * def Metropolis(sigma,J,B,T,iterations):
+ */
+  __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s__Metropolis), 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__Metropolis, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":17
+ * import Metropolis
+ * 
+ * def Metropolis(sigma,J,B,T,iterations):             # <<<<<<<<<<<<<<
+ *     start=time.time()
+ * 
+ */
+  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_10Metropolis_Metropolis, NULL, __pyx_n_s__Metropolis); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(__pyx_t_1);
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s__Metropolis, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
+  /* "Metropolis.pyx":1
+ * #!/usr/bin/env python             # <<<<<<<<<<<<<<
+ * #
+ * # Ising2D model in serial mode
+ */
+  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_GOTREF(((PyObject *)__pyx_t_1));
+  if (PyObject_SetAttr(__pyx_m, __pyx_n_s____test__, ((PyObject *)__pyx_t_1)) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+  __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
+  goto __pyx_L0;
+  __pyx_L1_error:;
+  __Pyx_XDECREF(__pyx_t_1);
+  __Pyx_XDECREF(__pyx_t_2);
+  if (__pyx_m) {
+    __Pyx_AddTraceback("init Metropolis", __pyx_clineno, __pyx_lineno, __pyx_filename);
+    Py_DECREF(__pyx_m); __pyx_m = 0;
+  } else if (!PyErr_Occurred()) {
+    PyErr_SetString(PyExc_ImportError, "init Metropolis");
+  }
+  __pyx_L0:;
+  __Pyx_RefNannyFinishContext();
+  #if PY_MAJOR_VERSION < 3
+  return;
+  #else
+  return __pyx_m;
+  #endif
+}
+
+/* Runtime support code */
+
+#if CYTHON_REFNANNY
+static __Pyx_RefNannyAPIStruct *__Pyx_RefNannyImportAPI(const char *modname) {
+    PyObject *m = NULL, *p = NULL;
+    void *r = NULL;
+    m = PyImport_ImportModule((char *)modname);
+    if (!m) goto end;
+    p = PyObject_GetAttrString(m, (char *)"RefNannyAPI");
+    if (!p) goto end;
+    r = PyLong_AsVoidPtr(p);
+end:
+    Py_XDECREF(p);
+    Py_XDECREF(m);
+    return (__Pyx_RefNannyAPIStruct *)r;
+}
+#endif /* CYTHON_REFNANNY */
+
+static PyObject *__Pyx_GetName(PyObject *dict, PyObject *name) {
+    PyObject *result;
+    result = PyObject_GetAttr(dict, name);
+    if (!result) {
+        if (dict != __pyx_b) {
+            PyErr_Clear();
+            result = PyObject_GetAttr(__pyx_b, name);
+        }
+        if (!result) {
+            PyErr_SetObject(PyExc_NameError, name);
+        }
+    }
+    return result;
+}
+
+static void __Pyx_RaiseArgtupleInvalid(
+    const char* func_name,
+    int exact,
+    Py_ssize_t num_min,
+    Py_ssize_t num_max,
+    Py_ssize_t num_found)
+{
+    Py_ssize_t num_expected;
+    const char *more_or_less;
+
+    if (num_found < num_min) {
+        num_expected = num_min;
+        more_or_less = "at least";
+    } else {
+        num_expected = num_max;
+        more_or_less = "at most";
+    }
+    if (exact) {
+        more_or_less = "exactly";
+    }
+    PyErr_Format(PyExc_TypeError,
+                 "%s() takes %s %"PY_FORMAT_SIZE_T"d positional argument%s (%"PY_FORMAT_SIZE_T"d given)",
+                 func_name, more_or_less, num_expected,
+                 (num_expected == 1) ? "" : "s", num_found);
+}
+
+static void __Pyx_RaiseDoubleKeywordsError(
+    const char* func_name,
+    PyObject* kw_name)
+{
+    PyErr_Format(PyExc_TypeError,
+        #if PY_MAJOR_VERSION >= 3
+        "%s() got multiple values for keyword argument '%U'", func_name, kw_name);
+        #else
+        "%s() got multiple values for keyword argument '%s'", func_name,
+        PyString_AS_STRING(kw_name));
+        #endif
+}
+
+static int __Pyx_ParseOptionalKeywords(
+    PyObject *kwds,
+    PyObject **argnames[],
+    PyObject *kwds2,
+    PyObject *values[],
+    Py_ssize_t num_pos_args,
+    const char* function_name)
+{
+    PyObject *key = 0, *value = 0;
+    Py_ssize_t pos = 0;
+    PyObject*** name;
+    PyObject*** first_kw_arg = argnames + num_pos_args;
+
+    while (PyDict_Next(kwds, &pos, &key, &value)) {
+        name = first_kw_arg;
+        while (*name && (**name != key)) name++;
+        if (*name) {
+            values[name-argnames] = value;
+        } else {
+            #if PY_MAJOR_VERSION < 3
+            if (unlikely(!PyString_CheckExact(key)) && unlikely(!PyString_Check(key))) {
+            #else
+            if (unlikely(!PyUnicode_CheckExact(key)) && unlikely(!PyUnicode_Check(key))) {
+            #endif
+                goto invalid_keyword_type;
+            } else {
+                for (name = first_kw_arg; *name; name++) {
+                    #if PY_MAJOR_VERSION >= 3
+                    if (PyUnicode_GET_SIZE(**name) == PyUnicode_GET_SIZE(key) &&
+                        PyUnicode_Compare(**name, key) == 0) break;
+                    #else
+                    if (PyString_GET_SIZE(**name) == PyString_GET_SIZE(key) &&
+                        _PyString_Eq(**name, key)) break;
+                    #endif
+                }
+                if (*name) {
+                    values[name-argnames] = value;
+                } else {
+                    /* unexpected keyword found */
+                    for (name=argnames; name != first_kw_arg; name++) {
+                        if (**name == key) goto arg_passed_twice;
+                        #if PY_MAJOR_VERSION >= 3
+                        if (PyUnicode_GET_SIZE(**name) == PyUnicode_GET_SIZE(key) &&
+                            PyUnicode_Compare(**name, key) == 0) goto arg_passed_twice;
+                        #else
+                        if (PyString_GET_SIZE(**name) == PyString_GET_SIZE(key) &&
+                            _PyString_Eq(**name, key)) goto arg_passed_twice;
+                        #endif
+                    }
+                    if (kwds2) {
+                        if (unlikely(PyDict_SetItem(kwds2, key, value))) goto bad;
+                    } else {
+                        goto invalid_keyword;
+                    }
+                }
+            }
+        }
+    }
+    return 0;
+arg_passed_twice:
+    __Pyx_RaiseDoubleKeywordsError(function_name, **name);
+    goto bad;
+invalid_keyword_type:
+    PyErr_Format(PyExc_TypeError,
+        "%s() keywords must be strings", function_name);
+    goto bad;
+invalid_keyword:
+    PyErr_Format(PyExc_TypeError,
+    #if PY_MAJOR_VERSION < 3
+        "%s() got an unexpected keyword argument '%s'",
+        function_name, PyString_AsString(key));
+    #else
+        "%s() got an unexpected keyword argument '%U'",
+        function_name, key);
+    #endif
+bad:
+    return -1;
+}
+
+static CYTHON_INLINE void __Pyx_RaiseNeedMoreValuesError(Py_ssize_t index) {
+    PyErr_Format(PyExc_ValueError,
+                 "need more than %"PY_FORMAT_SIZE_T"d value%s to unpack",
+                 index, (index == 1) ? "" : "s");
+}
+
+static CYTHON_INLINE void __Pyx_RaiseTooManyValuesError(Py_ssize_t expected) {
+    PyErr_Format(PyExc_ValueError,
+                 "too many values to unpack (expected %"PY_FORMAT_SIZE_T"d)", expected);
+}
+
+static int __Pyx_IternextUnpackEndCheck(PyObject *retval, Py_ssize_t expected) {
+    if (unlikely(retval)) {
+        Py_DECREF(retval);
+        __Pyx_RaiseTooManyValuesError(expected);
+        return -1;
+    } else if (PyErr_Occurred()) {
+        if (likely(PyErr_ExceptionMatches(PyExc_StopIteration))) {
+            PyErr_Clear();
+            return 0;
+        } else {
+            return -1;
+        }
+    }
+    return 0;
+}
+
+static PyObject *__Pyx_Import(PyObject *name, PyObject *from_list, long level) {
+    PyObject *py_import = 0;
+    PyObject *empty_list = 0;
+    PyObject *module = 0;
+    PyObject *global_dict = 0;
+    PyObject *empty_dict = 0;
+    PyObject *list;
+    py_import = __Pyx_GetAttrString(__pyx_b, "__import__");
+    if (!py_import)
+        goto bad;
+    if (from_list)
+        list = from_list;
+    else {
+        empty_list = PyList_New(0);
+        if (!empty_list)
+            goto bad;
+        list = empty_list;
+    }
+    global_dict = PyModule_GetDict(__pyx_m);
+    if (!global_dict)
+        goto bad;
+    empty_dict = PyDict_New();
+    if (!empty_dict)
+        goto bad;
+    #if PY_VERSION_HEX >= 0x02050000
+    {
+        PyObject *py_level = PyInt_FromLong(level);
+        if (!py_level)
+            goto bad;
+        module = PyObject_CallFunctionObjArgs(py_import,
+            name, global_dict, empty_dict, list, py_level, NULL);
+        Py_DECREF(py_level);
+    }
+    #else
+    if (level>0) {
+        PyErr_SetString(PyExc_RuntimeError, "Relative import is not supported for Python <=2.4.");
+        goto bad;
+    }
+    module = PyObject_CallFunctionObjArgs(py_import,
+        name, global_dict, empty_dict, list, NULL);
+    #endif
+bad:
+    Py_XDECREF(empty_list);
+    Py_XDECREF(py_import);
+    Py_XDECREF(empty_dict);
+    return module;
+}
+
+static CYTHON_INLINE unsigned char __Pyx_PyInt_AsUnsignedChar(PyObject* x) {
+    const unsigned char neg_one = (unsigned char)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(unsigned char) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(unsigned char)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to unsigned char" :
+                    "value too large to convert to unsigned char");
+            }
+            return (unsigned char)-1;
+        }
+        return (unsigned char)val;
+    }
+    return (unsigned char)__Pyx_PyInt_AsUnsignedLong(x);
+}
+
+static CYTHON_INLINE unsigned short __Pyx_PyInt_AsUnsignedShort(PyObject* x) {
+    const unsigned short neg_one = (unsigned short)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(unsigned short) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(unsigned short)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to unsigned short" :
+                    "value too large to convert to unsigned short");
+            }
+            return (unsigned short)-1;
+        }
+        return (unsigned short)val;
+    }
+    return (unsigned short)__Pyx_PyInt_AsUnsignedLong(x);
+}
+
+static CYTHON_INLINE unsigned int __Pyx_PyInt_AsUnsignedInt(PyObject* x) {
+    const unsigned int neg_one = (unsigned int)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(unsigned int) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(unsigned int)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to unsigned int" :
+                    "value too large to convert to unsigned int");
+            }
+            return (unsigned int)-1;
+        }
+        return (unsigned int)val;
+    }
+    return (unsigned int)__Pyx_PyInt_AsUnsignedLong(x);
+}
+
+static CYTHON_INLINE char __Pyx_PyInt_AsChar(PyObject* x) {
+    const char neg_one = (char)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(char) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(char)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to char" :
+                    "value too large to convert to char");
+            }
+            return (char)-1;
+        }
+        return (char)val;
+    }
+    return (char)__Pyx_PyInt_AsLong(x);
+}
+
+static CYTHON_INLINE short __Pyx_PyInt_AsShort(PyObject* x) {
+    const short neg_one = (short)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(short) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(short)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to short" :
+                    "value too large to convert to short");
+            }
+            return (short)-1;
+        }
+        return (short)val;
+    }
+    return (short)__Pyx_PyInt_AsLong(x);
+}
+
+static CYTHON_INLINE int __Pyx_PyInt_AsInt(PyObject* x) {
+    const int neg_one = (int)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(int) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(int)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to int" :
+                    "value too large to convert to int");
+            }
+            return (int)-1;
+        }
+        return (int)val;
+    }
+    return (int)__Pyx_PyInt_AsLong(x);
+}
+
+static CYTHON_INLINE signed char __Pyx_PyInt_AsSignedChar(PyObject* x) {
+    const signed char neg_one = (signed char)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(signed char) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(signed char)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to signed char" :
+                    "value too large to convert to signed char");
+            }
+            return (signed char)-1;
+        }
+        return (signed char)val;
+    }
+    return (signed char)__Pyx_PyInt_AsSignedLong(x);
+}
+
+static CYTHON_INLINE signed short __Pyx_PyInt_AsSignedShort(PyObject* x) {
+    const signed short neg_one = (signed short)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(signed short) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(signed short)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to signed short" :
+                    "value too large to convert to signed short");
+            }
+            return (signed short)-1;
+        }
+        return (signed short)val;
+    }
+    return (signed short)__Pyx_PyInt_AsSignedLong(x);
+}
+
+static CYTHON_INLINE signed int __Pyx_PyInt_AsSignedInt(PyObject* x) {
+    const signed int neg_one = (signed int)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(signed int) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(signed int)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to signed int" :
+                    "value too large to convert to signed int");
+            }
+            return (signed int)-1;
+        }
+        return (signed int)val;
+    }
+    return (signed int)__Pyx_PyInt_AsSignedLong(x);
+}
+
+static CYTHON_INLINE int __Pyx_PyInt_AsLongDouble(PyObject* x) {
+    const int neg_one = (int)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+    if (sizeof(int) < sizeof(long)) {
+        long val = __Pyx_PyInt_AsLong(x);
+        if (unlikely(val != (long)(int)val)) {
+            if (!unlikely(val == -1 && PyErr_Occurred())) {
+                PyErr_SetString(PyExc_OverflowError,
+                    (is_unsigned && unlikely(val < 0)) ?
+                    "can't convert negative value to int" :
+                    "value too large to convert to int");
+            }
+            return (int)-1;
+        }
+        return (int)val;
+    }
+    return (int)__Pyx_PyInt_AsLong(x);
+}
+
+static CYTHON_INLINE unsigned long __Pyx_PyInt_AsUnsignedLong(PyObject* x) {
+    const unsigned long neg_one = (unsigned long)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+#if PY_VERSION_HEX < 0x03000000
+    if (likely(PyInt_Check(x))) {
+        long val = PyInt_AS_LONG(x);
+        if (is_unsigned && unlikely(val < 0)) {
+            PyErr_SetString(PyExc_OverflowError,
+                            "can't convert negative value to unsigned long");
+            return (unsigned long)-1;
+        }
+        return (unsigned long)val;
+    } else
+#endif
+    if (likely(PyLong_Check(x))) {
+        if (is_unsigned) {
+            if (unlikely(Py_SIZE(x) < 0)) {
+                PyErr_SetString(PyExc_OverflowError,
+                                "can't convert negative value to unsigned long");
+                return (unsigned long)-1;
+            }
+            return (unsigned long)PyLong_AsUnsignedLong(x);
+        } else {
+            return (unsigned long)PyLong_AsLong(x);
+        }
+    } else {
+        unsigned long val;
+        PyObject *tmp = __Pyx_PyNumber_Int(x);
+        if (!tmp) return (unsigned long)-1;
+        val = __Pyx_PyInt_AsUnsignedLong(tmp);
+        Py_DECREF(tmp);
+        return val;
+    }
+}
+
+static CYTHON_INLINE unsigned PY_LONG_LONG __Pyx_PyInt_AsUnsignedLongLong(PyObject* x) {
+    const unsigned PY_LONG_LONG neg_one = (unsigned PY_LONG_LONG)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+#if PY_VERSION_HEX < 0x03000000
+    if (likely(PyInt_Check(x))) {
+        long val = PyInt_AS_LONG(x);
+        if (is_unsigned && unlikely(val < 0)) {
+            PyErr_SetString(PyExc_OverflowError,
+                            "can't convert negative value to unsigned PY_LONG_LONG");
+            return (unsigned PY_LONG_LONG)-1;
+        }
+        return (unsigned PY_LONG_LONG)val;
+    } else
+#endif
+    if (likely(PyLong_Check(x))) {
+        if (is_unsigned) {
+            if (unlikely(Py_SIZE(x) < 0)) {
+                PyErr_SetString(PyExc_OverflowError,
+                                "can't convert negative value to unsigned PY_LONG_LONG");
+                return (unsigned PY_LONG_LONG)-1;
+            }
+            return (unsigned PY_LONG_LONG)PyLong_AsUnsignedLongLong(x);
+        } else {
+            return (unsigned PY_LONG_LONG)PyLong_AsLongLong(x);
+        }
+    } else {
+        unsigned PY_LONG_LONG val;
+        PyObject *tmp = __Pyx_PyNumber_Int(x);
+        if (!tmp) return (unsigned PY_LONG_LONG)-1;
+        val = __Pyx_PyInt_AsUnsignedLongLong(tmp);
+        Py_DECREF(tmp);
+        return val;
+    }
+}
+
+static CYTHON_INLINE long __Pyx_PyInt_AsLong(PyObject* x) {
+    const long neg_one = (long)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+#if PY_VERSION_HEX < 0x03000000
+    if (likely(PyInt_Check(x))) {
+        long val = PyInt_AS_LONG(x);
+        if (is_unsigned && unlikely(val < 0)) {
+            PyErr_SetString(PyExc_OverflowError,
+                            "can't convert negative value to long");
+            return (long)-1;
+        }
+        return (long)val;
+    } else
+#endif
+    if (likely(PyLong_Check(x))) {
+        if (is_unsigned) {
+            if (unlikely(Py_SIZE(x) < 0)) {
+                PyErr_SetString(PyExc_OverflowError,
+                                "can't convert negative value to long");
+                return (long)-1;
+            }
+            return (long)PyLong_AsUnsignedLong(x);
+        } else {
+            return (long)PyLong_AsLong(x);
+        }
+    } else {
+        long val;
+        PyObject *tmp = __Pyx_PyNumber_Int(x);
+        if (!tmp) return (long)-1;
+        val = __Pyx_PyInt_AsLong(tmp);
+        Py_DECREF(tmp);
+        return val;
+    }
+}
+
+static CYTHON_INLINE PY_LONG_LONG __Pyx_PyInt_AsLongLong(PyObject* x) {
+    const PY_LONG_LONG neg_one = (PY_LONG_LONG)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+#if PY_VERSION_HEX < 0x03000000
+    if (likely(PyInt_Check(x))) {
+        long val = PyInt_AS_LONG(x);
+        if (is_unsigned && unlikely(val < 0)) {
+            PyErr_SetString(PyExc_OverflowError,
+                            "can't convert negative value to PY_LONG_LONG");
+            return (PY_LONG_LONG)-1;
+        }
+        return (PY_LONG_LONG)val;
+    } else
+#endif
+    if (likely(PyLong_Check(x))) {
+        if (is_unsigned) {
+            if (unlikely(Py_SIZE(x) < 0)) {
+                PyErr_SetString(PyExc_OverflowError,
+                                "can't convert negative value to PY_LONG_LONG");
+                return (PY_LONG_LONG)-1;
+            }
+            return (PY_LONG_LONG)PyLong_AsUnsignedLongLong(x);
+        } else {
+            return (PY_LONG_LONG)PyLong_AsLongLong(x);
+        }
+    } else {
+        PY_LONG_LONG val;
+        PyObject *tmp = __Pyx_PyNumber_Int(x);
+        if (!tmp) return (PY_LONG_LONG)-1;
+        val = __Pyx_PyInt_AsLongLong(tmp);
+        Py_DECREF(tmp);
+        return val;
+    }
+}
+
+static CYTHON_INLINE signed long __Pyx_PyInt_AsSignedLong(PyObject* x) {
+    const signed long neg_one = (signed long)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+#if PY_VERSION_HEX < 0x03000000
+    if (likely(PyInt_Check(x))) {
+        long val = PyInt_AS_LONG(x);
+        if (is_unsigned && unlikely(val < 0)) {
+            PyErr_SetString(PyExc_OverflowError,
+                            "can't convert negative value to signed long");
+            return (signed long)-1;
+        }
+        return (signed long)val;
+    } else
+#endif
+    if (likely(PyLong_Check(x))) {
+        if (is_unsigned) {
+            if (unlikely(Py_SIZE(x) < 0)) {
+                PyErr_SetString(PyExc_OverflowError,
+                                "can't convert negative value to signed long");
+                return (signed long)-1;
+            }
+            return (signed long)PyLong_AsUnsignedLong(x);
+        } else {
+            return (signed long)PyLong_AsLong(x);
+        }
+    } else {
+        signed long val;
+        PyObject *tmp = __Pyx_PyNumber_Int(x);
+        if (!tmp) return (signed long)-1;
+        val = __Pyx_PyInt_AsSignedLong(tmp);
+        Py_DECREF(tmp);
+        return val;
+    }
+}
+
+static CYTHON_INLINE signed PY_LONG_LONG __Pyx_PyInt_AsSignedLongLong(PyObject* x) {
+    const signed PY_LONG_LONG neg_one = (signed PY_LONG_LONG)-1, const_zero = 0;
+    const int is_unsigned = neg_one > const_zero;
+#if PY_VERSION_HEX < 0x03000000
+    if (likely(PyInt_Check(x))) {
+        long val = PyInt_AS_LONG(x);
+        if (is_unsigned && unlikely(val < 0)) {
+            PyErr_SetString(PyExc_OverflowError,
+                            "can't convert negative value to signed PY_LONG_LONG");
+            return (signed PY_LONG_LONG)-1;
+        }
+        return (signed PY_LONG_LONG)val;
+    } else
+#endif
+    if (likely(PyLong_Check(x))) {
+        if (is_unsigned) {
+            if (unlikely(Py_SIZE(x) < 0)) {
+                PyErr_SetString(PyExc_OverflowError,
+                                "can't convert negative value to signed PY_LONG_LONG");
+                return (signed PY_LONG_LONG)-1;
+            }
+            return (signed PY_LONG_LONG)PyLong_AsUnsignedLongLong(x);
+        } else {
+            return (signed PY_LONG_LONG)PyLong_AsLongLong(x);
+        }
+    } else {
+        signed PY_LONG_LONG val;
+        PyObject *tmp = __Pyx_PyNumber_Int(x);
+        if (!tmp) return (signed PY_LONG_LONG)-1;
+        val = __Pyx_PyInt_AsSignedLongLong(tmp);
+        Py_DECREF(tmp);
+        return val;
+    }
+}
+
+static int __Pyx_check_binary_version(void) {
+    char ctversion[4], rtversion[4];
+    PyOS_snprintf(ctversion, 4, "%d.%d", PY_MAJOR_VERSION, PY_MINOR_VERSION);
+    PyOS_snprintf(rtversion, 4, "%s", Py_GetVersion());
+    if (ctversion[0] != rtversion[0] || ctversion[2] != rtversion[2]) {
+        char message[200];
+        PyOS_snprintf(message, sizeof(message),
+                      "compiletime version %s of module '%.100s' "
+                      "does not match runtime version %s",
+                      ctversion, __Pyx_MODULE_NAME, rtversion);
+        #if PY_VERSION_HEX < 0x02050000
+        return PyErr_Warn(NULL, message);
+        #else
+        return PyErr_WarnEx(NULL, message, 1);
+        #endif
+    }
+    return 0;
+}
+
+#include "compile.h"
+#include "frameobject.h"
+#include "traceback.h"
+
+static void __Pyx_AddTraceback(const char *funcname, int __pyx_clineno,
+                               int __pyx_lineno, const char *__pyx_filename) {
+    PyObject *py_srcfile = 0;
+    PyObject *py_funcname = 0;
+    PyObject *py_globals = 0;
+    PyCodeObject *py_code = 0;
+    PyFrameObject *py_frame = 0;
+
+    #if PY_MAJOR_VERSION < 3
+    py_srcfile = PyString_FromString(__pyx_filename);
+    #else
+    py_srcfile = PyUnicode_FromString(__pyx_filename);
+    #endif
+    if (!py_srcfile) goto bad;
+    if (__pyx_clineno) {
+        #if PY_MAJOR_VERSION < 3
+        py_funcname = PyString_FromFormat( "%s (%s:%d)", funcname, __pyx_cfilenm, __pyx_clineno);
+        #else
+        py_funcname = PyUnicode_FromFormat( "%s (%s:%d)", funcname, __pyx_cfilenm, __pyx_clineno);
+        #endif
+    }
+    else {
+        #if PY_MAJOR_VERSION < 3
+        py_funcname = PyString_FromString(funcname);
+        #else
+        py_funcname = PyUnicode_FromString(funcname);
+        #endif
+    }
+    if (!py_funcname) goto bad;
+    py_globals = PyModule_GetDict(__pyx_m);
+    if (!py_globals) goto bad;
+    py_code = PyCode_New(
+        0,            /*int argcount,*/
+        #if PY_MAJOR_VERSION >= 3
+        0,            /*int kwonlyargcount,*/
+        #endif
+        0,            /*int nlocals,*/
+        0,            /*int stacksize,*/
+        0,            /*int flags,*/
+        __pyx_empty_bytes, /*PyObject *code,*/
+        __pyx_empty_tuple,  /*PyObject *consts,*/
+        __pyx_empty_tuple,  /*PyObject *names,*/
+        __pyx_empty_tuple,  /*PyObject *varnames,*/
+        __pyx_empty_tuple,  /*PyObject *freevars,*/
+        __pyx_empty_tuple,  /*PyObject *cellvars,*/
+        py_srcfile,   /*PyObject *filename,*/
+        py_funcname,  /*PyObject *name,*/
+        __pyx_lineno,   /*int firstlineno,*/
+        __pyx_empty_bytes  /*PyObject *lnotab*/
+    );
+    if (!py_code) goto bad;
+    py_frame = PyFrame_New(
+        PyThreadState_GET(), /*PyThreadState *tstate,*/
+        py_code,             /*PyCodeObject *code,*/
+        py_globals,          /*PyObject *globals,*/
+        0                    /*PyObject *locals*/
+    );
+    if (!py_frame) goto bad;
+    py_frame->f_lineno = __pyx_lineno;
+    PyTraceBack_Here(py_frame);
+bad:
+    Py_XDECREF(py_srcfile);
+    Py_XDECREF(py_funcname);
+    Py_XDECREF(py_code);
+    Py_XDECREF(py_frame);
+}
+
+static int __Pyx_InitStrings(__Pyx_StringTabEntry *t) {
+    while (t->p) {
+        #if PY_MAJOR_VERSION < 3
+        if (t->is_unicode) {
+            *t->p = PyUnicode_DecodeUTF8(t->s, t->n - 1, NULL);
+        } else if (t->intern) {
+            *t->p = PyString_InternFromString(t->s);
+        } else {
+            *t->p = PyString_FromStringAndSize(t->s, t->n - 1);
+        }
+        #else  /* Python 3+ has unicode identifiers */
+        if (t->is_unicode | t->is_str) {
+            if (t->intern) {
+                *t->p = PyUnicode_InternFromString(t->s);
+            } else if (t->encoding) {
+                *t->p = PyUnicode_Decode(t->s, t->n - 1, t->encoding, NULL);
+            } else {
+                *t->p = PyUnicode_FromStringAndSize(t->s, t->n - 1);
+            }
+        } else {
+            *t->p = PyBytes_FromStringAndSize(t->s, t->n - 1);
+        }
+        #endif
+        if (!*t->p)
+            return -1;
+        ++t;
+    }
+    return 0;
+}
+
+/* Type Conversion Functions */
+
+static CYTHON_INLINE int __Pyx_PyObject_IsTrue(PyObject* x) {
+   int is_true = x == Py_True;
+   if (is_true | (x == Py_False) | (x == Py_None)) return is_true;
+   else return PyObject_IsTrue(x);
+}
+
+static CYTHON_INLINE PyObject* __Pyx_PyNumber_Int(PyObject* x) {
+  PyNumberMethods *m;
+  const char *name = NULL;
+  PyObject *res = NULL;
+#if PY_VERSION_HEX < 0x03000000
+  if (PyInt_Check(x) || PyLong_Check(x))
+#else
+  if (PyLong_Check(x))
+#endif
+    return Py_INCREF(x), x;
+  m = Py_TYPE(x)->tp_as_number;
+#if PY_VERSION_HEX < 0x03000000
+  if (m && m->nb_int) {
+    name = "int";
+    res = PyNumber_Int(x);
+  }
+  else if (m && m->nb_long) {
+    name = "long";
+    res = PyNumber_Long(x);
+  }
+#else
+  if (m && m->nb_int) {
+    name = "int";
+    res = PyNumber_Long(x);
+  }
+#endif
+  if (res) {
+#if PY_VERSION_HEX < 0x03000000
+    if (!PyInt_Check(res) && !PyLong_Check(res)) {
+#else
+    if (!PyLong_Check(res)) {
+#endif
+      PyErr_Format(PyExc_TypeError,
+                   "__%s__ returned non-%s (type %.200s)",
+                   name, name, Py_TYPE(res)->tp_name);
+      Py_DECREF(res);
+      return NULL;
+    }
+  }
+  else if (!PyErr_Occurred()) {
+    PyErr_SetString(PyExc_TypeError,
+                    "an integer is required");
+  }
+  return res;
+}
+
+static CYTHON_INLINE Py_ssize_t __Pyx_PyIndex_AsSsize_t(PyObject* b) {
+  Py_ssize_t ival;
+  PyObject* x = PyNumber_Index(b);
+  if (!x) return -1;
+  ival = PyInt_AsSsize_t(x);
+  Py_DECREF(x);
+  return ival;
+}
+
+static CYTHON_INLINE PyObject * __Pyx_PyInt_FromSize_t(size_t ival) {
+#if PY_VERSION_HEX < 0x02050000
+   if (ival <= LONG_MAX)
+       return PyInt_FromLong((long)ival);
+   else {
+       unsigned char *bytes = (unsigned char *) &ival;
+       int one = 1; int little = (int)*(unsigned char*)&one;
+       return _PyLong_FromByteArray(bytes, sizeof(size_t), little, 0);
+   }
+#else
+   return PyInt_FromSize_t(ival);
+#endif
+}
+
+static CYTHON_INLINE size_t __Pyx_PyInt_AsSize_t(PyObject* x) {
+   unsigned PY_LONG_LONG val = __Pyx_PyInt_AsUnsignedLongLong(x);
+   if (unlikely(val == (unsigned PY_LONG_LONG)-1 && PyErr_Occurred())) {
+       return (size_t)-1;
+   } else if (unlikely(val != (unsigned PY_LONG_LONG)(size_t)val)) {
+       PyErr_SetString(PyExc_OverflowError,
+                       "value too large to convert to size_t");
+       return (size_t)-1;
+   }
+   return (size_t)val;
+}
+
+
+#endif /* Py_PYTHON_H */
Index: Ising/Cython/Ising2D.py
===================================================================
--- Ising/Cython/Ising2D.py	(revision 0)
+++ Ising/Cython/Ising2D.py	(revision 18)
@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+import Metropolis
+from Metropolis import Metropolis
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep="])
+    except getopt.GetoptError:
+        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]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = arg
+        elif opt in ("-p", "--tempstep"):
+            Tstep = numpy.uint32(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+      
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
+
+    E=[]
+    M=[]
+
+    for T in Trange:
+        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
+        sigma=numpy.copy(sigmaIn)
+        duration=Metropolis(sigma,J,B,T,Iterations)
+        E=numpy.append(E,Energy(sigma,J))
+        M=numpy.append(M,Magnetization(sigma,B))
+        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+    
+        print "CPU Time : %f" % (duration)
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    # Save output
+    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/Cython/Ising2D.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Cython/setup.py
===================================================================
--- Ising/Cython/setup.py	(revision 0)
+++ Ising/Cython/setup.py	(revision 18)
@@ -0,0 +1,14 @@
+# http://docs.cython.org/src/userguide/tutorial.html
+
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Distutils import build_ext
+
+setup(
+    cmdclass = {'build_ext': build_ext},
+    ext_modules = [Extension("Metropolis", ["Metropolis.pyx"],
+                   extra_compile_args=["-O3",])]
+)
+
+# To compile
+# python setup.py build_ext --inplace
Index: Ising/Cython/Metropolis.pyx
===================================================================
--- Ising/Cython/Metropolis.pyx	(revision 0)
+++ Ising/Cython/Metropolis.pyx	(revision 18)
@@ -0,0 +1,34 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+import Metropolis
+
+def Metropolis(sigma,J,B,T,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)

Property changes on: Ising/Cython/Metropolis.pyx
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Numpy-C/array_module_np.c
===================================================================
--- Ising/Numpy-C/array_module_np.c	(revision 0)
+++ Ising/Numpy-C/array_module_np.c	(revision 18)
@@ -0,0 +1,177 @@
+/*  Example of change values in array in the Numpy-C-API. */
+/* http://scipy-lectures.github.io/advanced/interfacing_with_c/interfacing_with_c.html */
+/* http://dsnra.jpl.nasa.gov/software/Python/numpydoc/numpy-13.html */
+/* http://enzo.googlecode.com/.../Grid_ConvertToNumpy.C */
+/* http://docs.scipy.org/doc/numpy/reference/c-api.array.html#data-access */
+
+/* 
+   What is absolutely necessary to know:
+   - 
+*/
+
+#include <Python.h>
+#include <numpy/arrayobject.h>
+#include <math.h>
+#include <stdio.h>
+#include <sys/time.h>
+
+/* All RNG numbers of Marsaglia are Unsigned int32 ! */
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+#define LFIB4 (t[c]=t[c]+t[c+58]+t[c+119]+t[++c+178])
+#define SWB   (t[c+237]=(x=t[c+15])-(y=t[++c]+(x<y)))
+#define UNI   (KISS*2.328306e-10)
+#define VNI   ((long) KISS)*4.656613e-10
+
+#define MWCfp (float)MWC * 2.328306435454494e-10f
+#define KISSfp (float)KISS * 2.328306435454494e-10f
+#define SHR3fp (float)SHR3 * 2.328306435454494e-10f
+
+// static PyObject* array_metropolis_np(PyObject* self, PyObject *args)
+static PyObject* array_metropolis_np(PyObject* self, PyObject *args)
+{
+
+  PyObject *in_array;
+  PyArrayObject *local;
+  unsigned int z,w;
+  long iterations=1;
+  int seed_w=521288629,seed_z=362436069,jcong=380116160,jsr=123456789;
+  int factor;
+  float B=0,J=0,T=0;
+  double elapsed=0;
+  unsigned long i;
+  struct timeval tv1,tv2;
+  struct timezone tz;
+
+  /*  parse single numpy array argument */
+  /* http://docs.python.org/release/1.5.2p2/ext/parseTuple.html */
+  /* O for in_array */
+  /* fff  for 3 float values (Coupling, Magnetic Field, Temperature)*/
+  /* l for 1 long int value (iterations) */
+  /* ii for 2 int values (2 seeds) */
+  if (!PyArg_ParseTuple(args, "Offflii", &in_array,
+			&J,&B,&T,&iterations,&seed_w,&seed_z))
+    {
+      printf("Not null argument provided !\n");
+      return NULL;
+    }
+  
+  /* display the extra parameters */
+  printf("Parameters for physics: J=%f\nB=%f\nT=%f\n",J,B,T);
+  printf("Parameters for simulation: iterations=%ld\nseed_w=%i\nseed_z=%i\n",
+	 iterations,seed_w,seed_z);
+  
+  local = (PyArrayObject *) PyArray_ContiguousFromAny(in_array, NPY_INT32, 2, 2);
+  
+
+  if (PyArray_NDIM(local)!=2) {
+    printf("This simulation is only for 2D arrays !\n");
+    return NULL;
+  }
+  
+  npy_intp size_x=*(npy_intp *)(PyArray_DIMS(local));
+  npy_intp size_y=*(npy_intp *)(PyArray_DIMS(local)+1);
+  /* Cast necessary, PyArg_ParseTuple() does not support any unsigned type ! */
+  w=(unsigned int)seed_w;
+  z=(unsigned int)seed_z;
+ 
+  gettimeofday(&tv1, &tz);
+  printf("Simulation started: ...");
+
+  for (i=0;i<(unsigned long)iterations;i++) {
+    
+    npy_intp x=MWC%size_x ;
+    npy_intp y=MWC%size_y ;
+
+    int p=*(int *)PyArray_GETPTR2(local,x,y);
+    
+    int d=*(int *)PyArray_GETPTR2(local,x,(y+1)%size_y);
+    int u=*(int *)PyArray_GETPTR2(local,x,(y-1)%size_y);
+    int l=*(int *)PyArray_GETPTR2(local,(x-1)%size_x,y);
+    int r=*(int *)PyArray_GETPTR2(local,(x+1)%size_x,y);
+    
+    float DeltaE=J*(float)p*(2.*(float)(u+d+l+r)+B);
+    
+    factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+
+    *(int *)PyArray_GETPTR2(local,x,y)=(factor*p);
+  }
+  gettimeofday(&tv2, &tz);    
+  elapsed=(double)((tv2.tv_sec-tv1.tv_sec) * 1000000L +	\
+                   (tv2.tv_usec-tv1.tv_usec))/1000000.; 
+  printf("done in %.2f seconds !\n",elapsed);
+  
+  return PyArray_Return(local);
+  
+}
+
+static PyObject* array_display_np(PyObject* self, PyObject *args)
+{
+  
+  PyObject *in_array;
+  PyArrayObject      *local;
+  unsigned int i,x,y;
+
+  /*  parse single numpy array argument */
+  /* http://docs.python.org/release/1.5.2p2/ext/parseTuple.html */
+  /* O for in_array */
+  if (!PyArg_ParseTuple(args, "O", &in_array)) {
+    printf("Not null argument provided !\n");
+    return NULL;
+  }
+  
+  local = (PyArrayObject *) PyArray_ContiguousFromObject(in_array, PyArray_NOTYPE, 2, 2);
+  
+  /* display the size of array */
+  printf("Array has %i dimension(s) : size of ",PyArray_NDIM(local));
+  if (PyArray_NDIM(local)==1) {
+      printf("%i\n",(int)*(PyArray_DIMS(local)));
+  }
+  else {
+    for (i=0;i<PyArray_NDIM(local)-1;i++) {
+      printf("%ix",(int)*(PyArray_DIMS(local)+i));
+    }
+    printf("%i\n",(int)*(PyArray_DIMS(local)+PyArray_NDIM(local)-1));
+  }
+  
+  /* display the array */
+  printf("[");
+  for (x=0;x<(unsigned int)(int)*(PyArray_DIMS(local));x++) {
+    printf("[");
+    for (y=0;y<(unsigned int)(int)*(PyArray_DIMS(local)+1);y++)
+      printf("%2i ",* (int *)PyArray_GETPTR2(local,x,y));
+    printf("]\n ");
+  }
+  printf("]\n");
+  
+  // Py_INCREF(local);
+
+  /*  construct the output array, like the input array */
+  return PyArray_Return(local);
+
+}
+
+
+/*  define functions in module */
+static PyMethodDef ArrayMethods[] =
+{
+     {"array_metropolis_np", array_metropolis_np, METH_VARARGS,
+         "evaluate a metropolis simulation on a numpy array"},
+     {"array_display_np", array_display_np, METH_VARARGS,
+         "evaluate a metropolis simulation on a numpy array"},
+     {NULL, NULL, 0, NULL}
+};
+
+/* module initialization */
+PyMODINIT_FUNC
+
+initarray_module_np(void)
+{
+     (void) Py_InitModule("array_module_np", ArrayMethods);
+     /* IMPORTANT: this must be called */
+     import_array();
+}
Index: Ising/Numpy-C/launch.py
===================================================================
--- Ising/Numpy-C/launch.py	(revision 0)
+++ Ising/Numpy-C/launch.py	(revision 18)
@@ -0,0 +1,33 @@
+import array_module_np
+import numpy as np
+import pylab
+
+# x = np.arange(0, 2 * np.pi, 0.1)
+# y = array_module_np.array_cos_np(x)
+# z = array_module_np.array_sin_np(x)
+# pylab.plot(x, y,x,z)
+# pylab.show()
+
+# x = np.arange(0, 2 * np.pi, 0.1)
+# print array_module_np.array_cos_arg_np(x,1)
+
+# x = np.arange(0,64,1).reshape(8,8).astype(np.int64)
+
+# print x
+
+# print array_module_np.array_operation_np(x,1,1,1,1000000000,2008,1010)
+
+SIZE=8
+
+x=np.where(np.random.randn(SIZE,SIZE)>0,1,-1).astype('int32')
+
+print x
+
+y=array_module_np.array_metropolis_np(x,1,0,0.1,10000000,128,256)
+
+array_module_np.array_display_np(x)
+
+
+print y,y.dtype
+
+print x,x.dtype
Index: Ising/Numpy-C/Ising2D.py.orig
===================================================================
--- Ising/Numpy-C/Ising2D.py.orig	(revision 0)
+++ Ising/Numpy-C/Ising2D.py.orig	(revision 18)
@@ -0,0 +1,159 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,J,B,T,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep="])
+    except getopt.GetoptError:
+        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]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = arg
+        elif opt in ("-p", "--tempstep"):
+            Tstep = numpy.uint32(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+      
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
+
+    E=[]
+    M=[]
+
+    for T in Trange:
+        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
+        sigma=numpy.copy(sigmaIn)
+        duration=Metropolis(sigma,J,B,T,Iterations)
+        E=numpy.append(E,Energy(sigma,J))
+        M=numpy.append(M,Magnetization(sigma,B))
+        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+    
+        print "CPU Time : %f" % (duration)
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    # Save output
+    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/Numpy-C/Ising2D.py.orig
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Numpy-C/Ising2D.py
===================================================================
--- Ising/Numpy-C/Ising2D.py	(revision 0)
+++ Ising/Numpy-C/Ising2D.py	(revision 18)
@@ -0,0 +1,166 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode using external C library
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+import array_module_np
+from numpy.random import randint as nprnd
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,J,B,T,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep="])
+    except getopt.GetoptError:
+        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]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = float(arg)
+        elif opt in ("-p", "--tempstep"):
+            Tstep = float(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+      
+    print "Coupling Factor J : %s" % J
+    print "Magnetic Field B :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax 
+    print "Temperature step : %s" % Tstep
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
+    print "Temperatures explored : %s" % Trange
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int32)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+
+    E=[]
+    M=[]
+
+    for T in Trange:
+        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
+        sigma=numpy.copy(sigmaIn)
+        # duration=Metropolis(sigma,J,B,T,Iterations)
+        SeedW,SeedZ=numpy.int32(nprnd(2**31-1)),numpy.int32(nprnd(2**31-1))
+        start=time.time()
+        array_module_np.array_metropolis_np(sigma,J,B,T,Iterations,SeedW,SeedZ)
+        duration=time.time()-start
+        E=numpy.append(E,Energy(sigma,J))
+        M=numpy.append(M,Magnetization(sigma,B))
+        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+    
+        print "CPU Time : %f" % (duration)
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    # Save output
+    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/Numpy-C/Ising2D.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Numpy-C/setup.py
===================================================================
--- Ising/Numpy-C/setup.py	(revision 0)
+++ Ising/Numpy-C/setup.py	(revision 18)
@@ -0,0 +1,8 @@
+from distutils.core import setup, Extension
+import numpy
+
+# define the extension module
+array_module_np = Extension('array_module_np', sources=['array_module_np.c'],include_dirs=[numpy.get_include()])
+
+# run the setup
+setup(ext_modules=[array_module_np])
Index: Ising/Numpy-C/build.sh
===================================================================
--- Ising/Numpy-C/build.sh	(revision 0)
+++ Ising/Numpy-C/build.sh	(revision 18)
@@ -0,0 +1,3 @@
+python setup.py build
+[ -e array_module_np.so ] && rm array_module_np.so
+ln -s build/lib.linux-x86_64-2.7/array_module_np.so array_module_np.so
\ No newline at end of file
Index: Ising/GPU/Ising2D-GPU-One.py
===================================================================
--- Ising/GPU/Ising2D-GPU-One.py	(revision 0)
+++ Ising/GPU/Ising2D-GPU-One.py	(revision 18)
@@ -0,0 +1,421 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+from numpy.random import randint as nprnd
+
+KERNEL_CODE_OPENCL="""
+
+// Marsaglia RNG very simple implementation
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+
+#define MWCfp MWC * 2.328306435454494e-10f
+#define KISSfp KISS * 2.328306435454494e-10f
+
+__kernel void MainLoopOne(__global char *s,float T,float J,float B,
+                          uint sizex,uint sizey,
+                          uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z;
+   uint w=seed_w;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*y];
+
+      int d=s[x+sizex*((y+1)%sizey)];
+      int u=s[x+sizex*((y-1)%sizey)];
+      int l=s[((x-1)%sizex)+sizex*y];
+      int r=s[((x+1)%sizex)+sizex*y];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey)] = (char)factor*p;
+   }
+   barrier(CLK_GLOBAL_MEM_FENCE);
+   
+}
+"""
+
+KERNEL_CODE_CUDA="""
+
+// Marsaglia RNG very simple implementation
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+
+#define MWCfp MWC * 2.328306435454494e-10f
+#define KISSfp KISS * 2.328306435454494e-10f
+
+__global__ void MainLoopOne(char *s,float T,float J,float B,
+                            uint sizex,uint sizey,
+                            uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z;
+   uint w=seed_w;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*y];
+
+      int d=s[x+sizex*((y+1)%sizey)];
+      int u=s[x+sizex*((y-1)%sizey)];
+      int l=s[((x-1)%sizex)+sizex*y];
+      int r=s[((x+1)%sizex)+sizex*y];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey)] = (char)factor*p;
+   }
+   __syncthreads();
+   
+}
+"""
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,T,J,B,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def MetropolisCuda(sigma,T,J,B,iterations,ParaStyle,Alu,Device):
+    
+    # Avec PyCUDA autoinit, rien a faire !
+
+    sigmaCU=cuda.InOut(sigma)
+
+    mod = SourceModule(KERNEL_CODE_CUDA)
+    
+    MetropolisCU=mod.get_function("MainLoopOne")
+
+    start = pycuda.driver.Event()
+    stop = pycuda.driver.Event()
+
+    SizeX,SizeY=sigma.shape
+
+    start.record()
+    start.synchronize()
+    MetropolisCU(sigmaCU,
+                 numpy.float32(T), 
+                 numpy.float32(J),
+                 numpy.float32(B),
+                 numpy.uint32(SizeX),
+                 numpy.uint32(SizeY),
+                 numpy.uint32(iterations),
+                 numpy.uint32(nprnd(2**31-1)),
+                 numpy.uint32(nprnd(2**31-1)),
+                 grid=(1,1),
+                 block=(1,1,1))
+    
+    print "%s with %i %s done" % (Alu,1,ParaStyle)
+
+    stop.record()
+    stop.synchronize()
+                
+    #elapsed = stop.time_since(start)*1e-3
+    elapsed = start.time_till(stop)*1e-3
+
+    return(elapsed)
+
+
+def MetropolisOpenCL(sigma,T,J,B,iterations,ParaStyle,Alu,Device):
+
+    # Initialisation des variables en les CASTant correctement
+    
+    # Je detecte un peripherique GPU dans la liste des peripheriques
+    # for platform in cl.get_platforms():
+    #     for device in platform.get_devices():
+    #             if cl.device_type.to_string(device.type)=='GPU':
+    #                     GPU=device
+    #print "GPU detected: ",device.name
+    
+    HasGPU=False
+    Id=1
+    # Primary Device selection based on Device Id
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            deviceType=cl.device_type.to_string(device.type)
+            if Id==Device and not HasGPU:
+                GPU=device
+                print "CPU/GPU selected: ",device.name
+                HasGPU=True
+            Id=Id+1
+    # Default Device selection based on ALU Type
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            deviceType=cl.device_type.to_string(device.type)
+            if deviceType=="GPU" and Alu=="GPU" and not HasGPU:
+                GPU=device
+                print "GPU selected: ",device.name
+                HasGPU=True
+            if deviceType=="CPU" and Alu=="CPU" and not HasGPU:        
+                GPU=device
+                print "CPU selected: ",device.name
+                HasGPU=True
+            
+    # Je cree le contexte et la queue pour son execution
+    # ctx = cl.create_some_context()
+    ctx = cl.Context([GPU])
+    queue = cl.CommandQueue(ctx,
+                            properties=cl.command_queue_properties.PROFILING_ENABLE)
+    
+    # Je recupere les flag possibles pour les buffers
+    mf = cl.mem_flags
+
+    # Attention au CAST ! C'est un int8 soit un char en OpenCL !
+    sigmaCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=sigma)
+   
+    MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build( \
+        options = "-cl-mad-enable -cl-fast-relaxed-math")
+
+    SizeX,SizeY=sigma.shape
+
+    if ParaStyle=='Blocks':
+        # Call OpenCL kernel
+        # (1,) is Global work size (only 1 work size)
+        # (1,) is local work size
+        CLLaunch=MetropolisCL.MainLoopOne(queue,(1,),None,
+                                          sigmaCL,
+                                          numpy.float32(T), 
+                                          numpy.float32(J),
+                                          numpy.float32(B),
+                                          numpy.uint32(SizeX),
+                                          numpy.uint32(SizeY),
+                                          numpy.uint32(iterations),
+                                          numpy.uint32(nprnd(2**31-1)),
+                                          numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,1,ParaStyle)
+    else:
+        # en OpenCL, necessaire de mettre un Global_id identique au local_id
+        CLLaunch=MetropolisCL.MainLoopOne(queue,(1,),(1,),
+                                          sigmaCL,
+                                          numpy.float32(T),
+                                          numpy.float32(J),
+                                          numpy.float32(B),
+                                          numpy.uint32(SizeX),
+                                          numpy.uint32(SizeY),
+                                          numpy.uint32(iterations),
+                                          numpy.uint32(nprnd(2**31-1)),
+                                          numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,1,ParaStyle)
+        
+    CLLaunch.wait()
+    cl.enqueue_copy(queue, sigma, sigmaCL).wait()
+    elapsed = 1e-9*(CLLaunch.profile.end - CLLaunch.profile.start)
+    sigmaCL.release()
+        
+    return(elapsed)
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Alu can be CPU or GPU
+    Alu='CPU'
+    # Id of GPU : 0 will use the first find !
+    Device=0
+    # GPU style can be Cuda (Nvidia implementation) or OpenCL
+    GpuStyle='OpenCL'
+    # Parallel distribution can be on Threads or Blocks
+    ParaStyle='Blocks'
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:a:d:g:t:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep=","alu=","gpustyle=","parastyle="])
+    except getopt.GetoptError:
+        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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Blocks> -t <ParaStyle>' % sys.argv[0]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Blocks> -t <ParaStyle>' % sys.argv[0]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = float(arg)
+        elif opt in ("-p", "--tempstep"):
+            Tstep = float(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+        elif opt in ("-a", "--alu"):
+            Alu = arg
+        elif opt in ("-d", "--device"):
+            Device = int(arg)
+        elif opt in ("-g", "--gpustyle"):
+            GpuStyle = arg
+        elif opt in ("-t", "--parastyle"):
+            ParaStyle = arg
+
+    if Alu=='CPU' and GpuStyle=='CUDA':
+        print "Alu can't be CPU for CUDA, set Alu to GPU"
+        Alu='GPU'
+
+    if ParaStyle not in ('Blocks','Threads','Hybrid'):
+        print "%s not exists, ParaStyle set as Threads !" % ParaStyle
+        ParaStyle='Blocks'
+   
+    print "Compute unit : %s" % Alu
+    print "Device Identification : %s" % Device
+    print "GpuStyle used : %s" % GpuStyle
+    print "Parallel Style used : %s" % ParaStyle
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    if GpuStyle=='CUDA':
+        # For PyCUDA import
+        import pycuda.driver as cuda
+        import pycuda.gpuarray as gpuarray
+        import pycuda.autoinit
+        from pycuda.compiler import SourceModule
+
+    if GpuStyle=='OpenCL':
+        # For PyOpenCL import
+        import pyopencl as cl
+        Id=1
+        for platform in cl.get_platforms():
+            for device in platform.get_devices():
+                deviceType=cl.device_type.to_string(device.type)
+                print "Device #%i of type %s : %s" % (Id,deviceType,device.name)
+                Id=Id+1
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
+
+    E=[]
+    M=[]
+
+    for T in Trange:
+        sigma=numpy.copy(sigmaIn)
+        if GpuStyle=='CUDA':
+            duration=MetropolisCuda(sigma,T,J,B,Iterations,ParaStyle,Alu,Device)
+        else:
+            duration=MetropolisOpenCL(sigma,T,J,B,Iterations,ParaStyle,Alu,Device)
+            
+        E=numpy.append(E,Energy(sigma,J))
+        M=numpy.append(M,Magnetization(sigma,B))
+        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+
+        print "CPU Time : %f" % (duration)
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
+       
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    # Save output
+    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/GPU/Ising2D-GPU-One.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/GPU/Ising2D-GPU.py-130527
===================================================================
--- Ising/GPU/Ising2D-GPU.py-130527	(revision 0)
+++ Ising/GPU/Ising2D-GPU.py-130527	(revision 18)
@@ -0,0 +1,500 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+import math
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+from numpy.random import randint as nprnd
+
+KERNEL_CODE_OPENCL="""
+
+// Marsaglia RNG very simple implementation
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+
+#define MWCfp MWC * 2.328306435454494e-10f
+#define KISSfp KISS * 2.328306435454494e-10f
+
+__kernel void MainLoopOne(__global char *s,float T,float J,
+                          uint sizex,uint sizey,
+                          uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z;
+   uint w=seed_w;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*y];
+
+      int d=s[x+sizex*((y+1)%sizey)];
+      int u=s[x+sizex*((y-1)%sizey)];
+      int l=s[((x-1)%sizex)+sizex*y];
+      int r=s[((x+1)%sizex)+sizex*y];
+
+      float DeltaE=2.0f*J*p*(u+d+l+r);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey)] = (char)factor*p;
+   }
+   barrier(CLK_GLOBAL_MEM_FENCE);      
+}
+
+__kernel void MainLoopGlobal(__global char *s,__global float *T,float J,
+                             uint sizex,uint sizey,
+                             uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z/(get_global_id(0)+1);
+   uint w=seed_w/(get_global_id(0)+1);
+   float t;
+   uint ind=get_global_id(0);
+
+   t=T[get_global_id(0)];
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*J*p*(u+d+l+r);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+      
+   }
+
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+
+__kernel void MainLoopLocal(__global char *s,__global float *T,float J,
+                            uint sizex,uint sizey,
+                            uint iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(get_local_id(0)+1);
+   uint w=seed_w/(get_local_id(0)+1);
+   //float t=T[get_local_id(0)+get_global_id(0)*get_local_size(0)];
+   //uint ind=get_local_id(0)+get_global_id(0)*get_local_size(0);
+   float t=T[get_local_id(0)];
+   uint ind=get_local_id(0);
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*J*p*(u+d+l+r);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+   }
+
+   barrier(CLK_LOCAL_MEM_FENCE);
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+"""
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,J,B,T,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def MetropolisOpenCL(sigma,J,B,T,iterations,jobs,ParaStyle,Alu,Device):
+
+    # Initialisation des variables en les CASTant correctement
+    
+    # Je detecte un peripherique GPU dans la liste des peripheriques
+    # for platform in cl.get_platforms():
+    #     for device in platform.get_devices():
+    #             if cl.device_type.to_string(device.type)=='GPU':
+    #                     GPU=device
+    #print "GPU detected: ",device.name
+    
+    HasGPU=False
+    Id=1
+    # Device selection based on choice (default is GPU)
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            if not HasGPU:
+                deviceType=cl.device_type.to_string(device.type)
+                if deviceType=="GPU" and Alu=="GPU" and Id==Device:
+                    GPU=device
+                    print "GPU selected: ",device.name
+                    HasGPU=True
+                if deviceType=="CPU" and Alu=="CPU":        
+                    GPU=device
+                    print "CPU selected: ",device.name
+                    HasGPU=True
+            Id=Id+1
+                                
+    # Je cree le contexte et la queue pour son execution
+    # ctx = cl.create_some_context()
+    ctx = cl.Context([GPU])
+    queue = cl.CommandQueue(ctx,
+                            properties=cl.command_queue_properties.PROFILING_ENABLE)
+    
+    # Je recupere les flag possibles pour les buffers
+    mf = cl.mem_flags
+
+    print sigma,sigma.shape
+        
+    # Attention au CAST ! C'est un int8 soit un char en OpenCL !
+    sigmaCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=sigma)
+   
+    MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build( \
+        options = "-cl-mad-enable -cl-fast-relaxed-math")
+
+    SizeX,SizeY=sigma.shape
+
+    if ParaStyle=='Blocks':
+        # Call OpenCL kernel
+        # (1,) is Global work size (only 1 work size)
+        # (1,) is local work size
+        CLLaunch=MetropolisCL.MainLoopOne(queue,(jobs,),None,
+                                          sigmaCL,
+                                          numpy.float32(T), 
+                                          numpy.float32(J),
+                                          numpy.uint32(SizeX),
+                                          numpy.uint32(SizeY),
+                                          numpy.uint32(iterations),
+                                          numpy.uint32(nprnd(2**31-1)),
+                                          numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+    else:
+        # en OpenCL, necessaire de mettre un Global_id identique au local_id
+        CLLaunch=MetropolisCL.MainLoopOne(queue,(jobs,),(jobs,),
+                                          sigmaCL,
+                                          numpy.float32(T), 
+                                          numpy.float32(J),
+                                          numpy.uint32(SizeX),
+                                          numpy.uint32(SizeY),
+                                          numpy.uint32(iterations),
+                                          numpy.uint32(nprnd(2**31-1)),
+                                          numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+        
+    CLLaunch.wait()
+    cl.enqueue_copy(queue, sigma, sigmaCL).wait()
+    elapsed = 1e-9*(CLLaunch.profile.end - CLLaunch.profile.start)
+    sigmaCL.release()
+        
+    return(elapsed)
+
+def MetropolisAllOpenCL(sigmaDict,J,B,TList,iterations,jobs,ParaStyle,Alu,Device):
+
+    # sigmaDict & Tlist are NOT respectively array & float
+    # sigmaDict : dict of array for each temperatoire
+    # TList : list of temperatures
+          
+    # Initialisation des variables en les CASTant correctement
+    
+    # Je detecte un peripherique GPU dans la liste des peripheriques
+    # for platform in cl.get_platforms():
+    #     for device in platform.get_devices():
+    #             if cl.device_type.to_string(device.type)=='GPU':
+    #                     GPU=device
+    #print "GPU detected: ",device.name
+    
+    HasGPU=False
+    Id=1
+    # Device selection based on choice (default is GPU)
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            if not HasGPU:
+                deviceType=cl.device_type.to_string(device.type)
+                if deviceType=="GPU" and Alu=="GPU" and Id==Device:
+                    GPU=device
+                    print "GPU selected: ",device.name
+                    HasGPU=True
+                if deviceType=="CPU" and Alu=="CPU":        
+                    GPU=device
+                    print "CPU selected: ",device.name
+                    HasGPU=True
+            Id=Id+1
+                                
+    # Je cree le contexte et la queue pour son execution
+    # ctx = cl.create_some_context()
+    ctx = cl.Context([GPU])
+    queue = cl.CommandQueue(ctx,
+                            properties=cl.command_queue_properties.PROFILING_ENABLE)
+    
+    # Je recupere les flag possibles pour les buffers
+    mf = cl.mem_flags
+
+    # Concatenate all sigma in single array
+    sigma=numpy.copy(sigmaDict[TList[0]])
+    for T in TList[1:]:
+        sigma=numpy.concatenate((sigma,sigmaDict[T]),axis=1)
+
+    print sigma,sigma.shape
+        
+    sigmaCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=sigma)
+    TCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=TList)
+   
+    MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build( \
+        options = "-cl-mad-enable -cl-fast-relaxed-math")
+
+    SizeX,SizeY=sigmaDict[TList[0]].shape
+
+    if ParaStyle=='Blocks':
+        # Call OpenCL kernel
+        # (1,) is Global work size (only 1 work size)
+        # (1,) is local work size
+        # SeedZCL is lattice translated in CL format
+        # SeedWCL is lattice translated in CL format
+        # step is number of iterations
+        CLLaunch=MetropolisCL.MainLoopGlobal(queue,(jobs,),None,
+                                             sigmaCL,
+                                             TCL, 
+                                             numpy.float32(J),
+                                             numpy.uint32(SizeX),
+                                             numpy.uint32(SizeY),
+                                             numpy.uint32(iterations),
+                                             numpy.uint32(nprnd(2**31-1)),
+                                             numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+    else:
+        blocks=int(math.sqrt(jobs))
+        # en OpenCL, necessaire de mettre un Global_id identique au local_id
+        CLLaunch=MetropolisCL.MainLoopLocal(queue,(jobs,),(2,),
+                                            sigmaCL,
+                                            TCL,
+                                            numpy.float32(J),
+                                            numpy.uint32(SizeX),
+                                            numpy.uint32(SizeY),
+                                            numpy.uint32(iterations),
+                                            numpy.uint32(nprnd(2**31-1)),
+                                            numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+        
+    CLLaunch.wait()
+    cl.enqueue_copy(queue, sigma, sigmaCL).wait()
+    elapsed = 1e-9*(CLLaunch.profile.end - CLLaunch.profile.start)
+    sigmaCL.release()
+
+    results=numpy.split(sigma,len(TList),axis=1)
+    for T in TList:
+        sigmaDict[T]=numpy.copy(results[numpy.nonzero(TList == T)[0][0]])
+
+    return(elapsed)
+
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Alu can be CPU or GPU
+    Alu='CPU'
+    # Id of GPU
+    Device=1
+    # GPU style can be Cuda (Nvidia implementation) or OpenCL
+    GpuStyle='OpenCL'
+    # Parallel distribution can be on Threads or Blocks
+    ParaStyle='Blocks'
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:a:d:g:t:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep=","alu=","gpustyle=","parastyle="])
+    except getopt.GetoptError:
+        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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Blocks> -t <ParaStyle>' % sys.argv[0]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Blocks> -t <ParaStyle>' % sys.argv[0]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = float(arg)
+        elif opt in ("-p", "--tempstep"):
+            Tstep = float(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+        elif opt in ("-a", "--alu"):
+            Alu = arg
+        elif opt in ("-d", "--device"):
+            Device = int(arg)
+        elif opt in ("-g", "--gpustyle"):
+            GpuStyle = arg
+        elif opt in ("-t", "--parastyle"):
+            ParaStyle = arg
+
+    if Alu=='CPU' and GpuStyle=='CUDA':
+        print "Alu can't be CPU for CUDA, set Alu to GPU"
+        Alu='GPU'
+
+    if ParaStyle not in ('Blocks','Threads','Hybrid'):
+        print "%s not exists, ParaStyle set as Threads !" % ParaStyle
+        ParaStyle='Blocks'
+   
+    print "Compute unit : %s" % Alu
+    print "Device Identification : %s" % Device
+    print "GpuStyle used : %s" % GpuStyle
+    print "Parallel Style used : %s" % ParaStyle
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    if GpuStyle=='CUDA':
+        # For PyCUDA import
+        import pycuda.driver as cuda
+        import pycuda.gpuarray as gpuarray
+        import pycuda.autoinit
+        from pycuda.compiler import SourceModule
+
+    if GpuStyle=='OpenCL':
+        # For PyOpenCL import
+        import pyopencl as cl
+        Id=1
+        for platform in cl.get_platforms():
+            for device in platform.get_devices():
+                deviceType=cl.device_type.to_string(device.type)
+                print "Device #%i of type %s : %s" % (Id,deviceType,device.name)
+                Id=Id+1
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    # La temperature est passee comme parametre, attention au CAST !
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep).astype(numpy.float32)
+
+    E=[]
+    M=[]
+
+    print Trange,Trange.shape
+
+    sigma={}
+    for T in Trange:
+        sigma[T]=numpy.copy(sigmaIn)
+
+    # For GPU, all process are launched
+    #MetropolisAllOpenCL(sigma,J,B,Trange,Iterations,len(Trange),
+    #                    ParaStyle,Alu,Device)
+    MetropolisAllOpenCL(sigma,J,B,Trange,Iterations,len(Trange),
+                        ParaStyle,Alu,Device)
+    
+    for T in Trange:        
+        ImageOutput(sigma[T],"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+
+    # if Curves:
+    #     DisplayCurves(Trange,E,M,J,B)
+
+    # # Save output
+    # numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/GPU/Ising2D-GPU.py-130527
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/GPU/Ising2D-GPU.py
===================================================================
--- Ising/GPU/Ising2D-GPU.py	(revision 0)
+++ Ising/GPU/Ising2D-GPU.py	(revision 18)
@@ -0,0 +1,726 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+import math
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+from numpy.random import randint as nprnd
+
+KERNEL_CODE_OPENCL="""
+
+// Marsaglia RNG very simple implementation
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+
+#define MWCfp MWC * 2.328306435454494e-10f
+#define KISSfp KISS * 2.328306435454494e-10f
+
+__kernel void MainLoopOne(__global char *s,float T,float J,float B,
+                          uint sizex,uint sizey,
+                          uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z;
+   uint w=seed_w;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*y];
+
+      int d=s[x+sizex*((y+1)%sizey)];
+      int u=s[x+sizex*((y-1)%sizey)];
+      int l=s[((x-1)%sizex)+sizex*y];
+      int r=s[((x+1)%sizex)+sizex*y];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey)] = (char)factor*p;
+   }
+   barrier(CLK_GLOBAL_MEM_FENCE);      
+}
+
+__kernel void MainLoopGlobal(__global char *s,__global float *T,float J,float B,
+                             uint sizex,uint sizey,
+                             uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z/(get_global_id(0)+1);
+   uint w=seed_w/(get_global_id(0)+1);
+   float t=T[get_global_id(0)];
+   uint ind=get_global_id(0);
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+      
+   }
+
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+
+__kernel void MainLoopHybrid(__global char *s,__global float *T,float J,float B,
+                             uint sizex,uint sizey,
+                             uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z/(get_group_id(0)*get_num_groups(0)+get_local_id(0)+1);
+   uint w=seed_w/(get_group_id(0)*get_num_groups(0)+get_local_id(0)+1);
+   float t=T[get_group_id(0)*get_num_groups(0)+get_local_id(0)];
+   uint ind=get_group_id(0)*get_num_groups(0)+get_local_id(0);
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+      
+   }
+
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+
+__kernel void MainLoopLocal(__global char *s,__global float *T,float J,float B,
+                            uint sizex,uint sizey,
+                            uint iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(get_local_id(0)+1);
+   uint w=seed_w/(get_local_id(0)+1);
+   float t=T[get_local_id(0)];
+   uint ind=get_local_id(0);
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+   }
+
+   barrier(CLK_LOCAL_MEM_FENCE);
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+"""
+
+KERNEL_CODE_CUDA="""
+
+// Marsaglia RNG very simple implementation
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+
+#define MWCfp MWC * 2.328306435454494e-10f
+#define KISSfp KISS * 2.328306435454494e-10f
+
+__global__ void MainLoopOne(char *s,float T,float J,float B,
+                            uint sizex,uint sizey,
+                            uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z;
+   uint w=seed_w;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*y];
+
+      int d=s[x+sizex*((y+1)%sizey)];
+      int u=s[x+sizex*((y-1)%sizey)];
+      int l=s[((x-1)%sizex)+sizex*y];
+      int r=s[((x+1)%sizex)+sizex*y];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey)] = (char)factor*p;
+   }
+   __syncthreads();
+
+}
+
+__global__ void MainLoopGlobal(char *s,float *T,float J,float B,
+                               uint sizex,uint sizey,
+                               uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z/(blockIdx.x+1);
+   uint w=seed_w/(blockIdx.x+1);
+   float t=T[blockIdx.x];
+   uint ind=blockIdx.x;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+      
+   }
+   __syncthreads();
+
+}
+
+__global__ void MainLoopHybrid(char *s,float *T,float J,float B,
+                               uint sizex,uint sizey,
+                               uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z/(blockDim.x*blockIdx.x+threadIdx.x+1);
+   uint w=seed_w/(blockDim.x*blockIdx.x+threadIdx.x+1);
+   float t=T[blockDim.x*blockIdx.x+threadIdx.x];
+   uint ind=blockDim.x*blockIdx.x+threadIdx.x;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+      
+   }
+   __syncthreads();
+
+}
+
+__global__ void MainLoopLocal(char *s,float *T,float J,float B,
+                              uint sizex,uint sizey,
+                              uint iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(threadIdx.x+1);
+   uint w=seed_w/(threadIdx.x+1);
+   float t=T[threadIdx.x];
+   uint ind=threadIdx.x;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+   }
+   __syncthreads();
+
+}
+"""
+
+# find prime factors of a number
+# Get for WWW :
+# http://pythonism.wordpress.com/2008/05/17/looking-at-factorisation-in-python/
+def PrimeFactors(x):
+    factorlist=numpy.array([]).astype('uint32')
+    loop=2
+    while loop<=x:
+        if x%loop==0:
+            x/=loop
+            factorlist=numpy.append(factorlist,[loop])
+        else:
+            loop+=1
+    return factorlist
+
+# Try to find the best thread number in Hybrid approach (Blocks&Threads)
+# output is thread number
+def BestThreadsNumber(jobs):
+    factors=PrimeFactors(jobs)
+    matrix=numpy.append([factors],[factors[::-1]],axis=0)
+    threads=1
+    for factor in matrix.transpose().ravel():
+        threads=threads*factor
+        if threads*threads>jobs:
+            break
+    return(long(threads))
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,T,J,B,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def MetropolisAllOpenCL(sigmaDict,TList,J,B,iterations,jobs,ParaStyle,Alu,Device):
+
+    # sigmaDict & Tlist are NOT respectively array & float
+    # sigmaDict : dict of array for each temperatoire
+    # TList : list of temperatures
+          
+    # Initialisation des variables en les CASTant correctement
+    
+    # Je detecte un peripherique GPU dans la liste des peripheriques
+
+    HasGPU=False
+    Id=1
+    # Primary Device selection based on Device Id
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            deviceType=cl.device_type.to_string(device.type)
+            if Id==Device and not HasGPU:
+                GPU=device
+                print "CPU/GPU selected: ",device.name
+                HasGPU=True
+            Id=Id+1
+    # Default Device selection based on ALU Type
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            deviceType=cl.device_type.to_string(device.type)
+            if deviceType=="GPU" and Alu=="GPU" and not HasGPU:
+                GPU=device
+                print "GPU selected: ",device.name
+                HasGPU=True
+            if deviceType=="CPU" and Alu=="CPU" and not HasGPU:        
+                GPU=device
+                print "CPU selected: ",device.name
+                HasGPU=True
+                                
+    # Je cree le contexte et la queue pour son execution
+    # ctx = cl.create_some_context()
+    ctx = cl.Context([GPU])
+    queue = cl.CommandQueue(ctx,
+                            properties=cl.command_queue_properties.PROFILING_ENABLE)
+    
+    # Je recupere les flag possibles pour les buffers
+    mf = cl.mem_flags
+
+    # Concatenate all sigma in single array
+    sigma=numpy.copy(sigmaDict[TList[0]])
+    for T in TList[1:]:
+        sigma=numpy.concatenate((sigma,sigmaDict[T]),axis=1)
+
+    sigmaCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=sigma)
+    TCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=TList)
+   
+    MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build( \
+        options = "-cl-mad-enable -cl-fast-relaxed-math")
+
+    SizeX,SizeY=sigmaDict[TList[0]].shape
+
+    if ParaStyle=='Blocks':
+        # Call OpenCL kernel
+        # (1,) is Global work size (only 1 work size)
+        # (1,) is local work size
+        # SeedZCL is lattice translated in CL format
+        # SeedWCL is lattice translated in CL format
+        # step is number of iterations
+        CLLaunch=MetropolisCL.MainLoopGlobal(queue,(jobs,),None,
+                                             sigmaCL,
+                                             TCL, 
+                                             numpy.float32(J),
+                                             numpy.float32(B),
+                                             numpy.uint32(SizeX),
+                                             numpy.uint32(SizeY),
+                                             numpy.uint32(iterations),
+                                             numpy.uint32(nprnd(2**31-1)),
+                                             numpy.uint32(nprnd(2**31-1)))
+        print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
+              (Alu,jobs,1,ParaStyle)
+    elif ParaStyle=='Threads':
+        # It's necessary to put a Local_ID equal to a Global_ID
+        # Jobs are to be considerated as global number of jobs to do
+        # And to be distributed to entities
+        # For example :
+        # G_ID=10 & L_ID=10 : 10 Threads on 1 UC
+        # G_ID=10 & L_ID=1 : 10 Threads on 1 UC
+        
+        CLLaunch=MetropolisCL.MainLoopLocal(queue,(jobs,),(jobs,),
+                                            sigmaCL,
+                                            TCL,
+                                            numpy.float32(J),
+                                            numpy.float32(B),
+                                            numpy.uint32(SizeX),
+                                            numpy.uint32(SizeY),
+                                            numpy.uint32(iterations),
+                                            numpy.uint32(nprnd(2**31-1)),
+                                            numpy.uint32(nprnd(2**31-1)))
+        print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
+              (Alu,1,jobs,ParaStyle)
+    else:
+        threads=BestThreadsNumber(jobs)
+        # en OpenCL, necessaire de mettre un Global_id identique au local_id
+        CLLaunch=MetropolisCL.MainLoopHybrid(queue,(jobs,),(threads,),
+                                             sigmaCL,
+                                             TCL,
+                                             numpy.float32(J),
+                                             numpy.float32(B),
+                                             numpy.uint32(SizeX),
+                                             numpy.uint32(SizeY),
+                                             numpy.uint32(iterations),
+                                             numpy.uint32(nprnd(2**31-1)),
+                                             numpy.uint32(nprnd(2**31-1)))
+        print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
+              (Alu,jobs/threads,threads,ParaStyle)
+        
+    CLLaunch.wait()
+    cl.enqueue_copy(queue, sigma, sigmaCL).wait()
+    elapsed = 1e-9*(CLLaunch.profile.end - CLLaunch.profile.start)
+    sigmaCL.release()
+
+    results=numpy.split(sigma,len(TList),axis=1)
+    for T in TList:
+        sigmaDict[T]=numpy.copy(results[numpy.nonzero(TList == T)[0][0]])
+
+    return(elapsed)
+
+def MetropolisAllCuda(sigmaDict,TList,J,B,iterations,jobs,ParaStyle,Alu,Device):
+
+    # sigmaDict & Tlist are NOT respectively array & float
+    # sigmaDict : dict of array for each temperatoire
+    # TList : list of temperatures
+          
+    # Avec PyCUDA autoinit, rien a faire !
+
+    mod = SourceModule(KERNEL_CODE_CUDA)
+    
+    MetropolisBlocksCuda=mod.get_function("MainLoopGlobal")
+    MetropolisThreadsCuda=mod.get_function("MainLoopLocal")
+    MetropolisHybridCuda=mod.get_function("MainLoopHybrid")
+
+    # Concatenate all sigma in single array
+    sigma=numpy.copy(sigmaDict[TList[0]])
+    for T in TList[1:]:
+        sigma=numpy.concatenate((sigma,sigmaDict[T]),axis=1)
+
+    sigmaCU=cuda.InOut(sigma)
+    TCU=cuda.InOut(TList)
+
+    SizeX,SizeY=sigmaDict[TList[0]].shape
+
+    start = pycuda.driver.Event()
+    stop = pycuda.driver.Event()
+    
+    start.record()
+    start.synchronize()
+    if ParaStyle=='Blocks':
+        # Call CUDA kernel
+        # (1,) is Global work size (only 1 work size)
+        # (1,) is local work size
+        # SeedZCL is lattice translated in CL format
+        # SeedWCL is lattice translated in CL format
+        # step is number of iterations
+        MetropolisBlocksCuda(sigmaCU,
+                             TCU, 
+                             numpy.float32(J),
+                             numpy.float32(B),
+                             numpy.uint32(SizeX),
+                             numpy.uint32(SizeY),
+                             numpy.uint32(iterations),
+                             numpy.uint32(nprnd(2**31-1)),
+                             numpy.uint32(nprnd(2**31-1)),
+                             grid=(jobs,1),block=(1,1,1))
+        print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
+              (Alu,jobs,1,ParaStyle)
+    elif ParaStyle=='Threads':
+        MetropolisThreadsCuda(sigmaCU,
+                              TCU,
+                              numpy.float32(J),
+                              numpy.float32(B),
+                              numpy.uint32(SizeX),
+                              numpy.uint32(SizeY),
+                              numpy.uint32(iterations),
+                              numpy.uint32(nprnd(2**31-1)),
+                              numpy.uint32(nprnd(2**31-1)),
+                              grid=(1,1),block=(jobs,1,1))
+        print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
+              (Alu,1,jobs,ParaStyle)
+    else:
+        threads=BestThreadsNumber(jobs)        
+        MetropolisHybridCuda(sigmaCU,
+                              TCU,
+                              numpy.float32(J),
+                              numpy.float32(B),
+                              numpy.uint32(SizeX),
+                              numpy.uint32(SizeY),
+                              numpy.uint32(iterations),
+                              numpy.uint32(nprnd(2**31-1)),
+                              numpy.uint32(nprnd(2**31-1)),
+                              grid=(jobs/threads,1),block=(threads,1,1))
+        print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
+              (Alu,jobs/threads,threads,ParaStyle)
+        
+    stop.record()
+    stop.synchronize()
+    elapsed = start.time_till(stop)*1e-3
+
+    results=numpy.split(sigma,len(TList),axis=1)
+    for T in TList:
+        sigmaDict[T]=numpy.copy(results[numpy.nonzero(TList == T)[0][0]])
+
+    return(elapsed)
+
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Alu can be CPU or GPU
+    Alu='CPU'
+    # Id of GPU : 0 
+    Device=0
+    # GPU style can be Cuda (Nvidia implementation) or OpenCL
+    GpuStyle='OpenCL'
+    # Parallel distribution can be on Threads or Blocks
+    ParaStyle='Blocks'
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:a:d:g:t:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep=","alu=","gpustyle=","parastyle="])
+    except getopt.GetoptError:
+        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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -t <Threads/Blocks>' % sys.argv[0]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -t <Threads/Blocks/Hybrid>' % sys.argv[0]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = float(arg)
+        elif opt in ("-p", "--tempstep"):
+            Tstep = float(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+        elif opt in ("-a", "--alu"):
+            Alu = arg
+        elif opt in ("-d", "--device"):
+            Device = int(arg)
+        elif opt in ("-g", "--gpustyle"):
+            GpuStyle = arg
+        elif opt in ("-t", "--parastyle"):
+            ParaStyle = arg
+
+    if Alu=='CPU' and GpuStyle=='CUDA':
+        print "Alu can't be CPU for CUDA, set Alu to GPU"
+        Alu='GPU'
+
+    if ParaStyle not in ('Blocks','Threads','Hybrid'):
+        print "%s not exists, ParaStyle set as Threads !" % ParaStyle
+        ParaStyle='Blocks'
+   
+    print "Compute unit : %s" % Alu
+    print "Device Identification : %s" % Device
+    print "GpuStyle used : %s" % GpuStyle
+    print "Parallel Style used : %s" % ParaStyle
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    if GpuStyle=='CUDA':
+        # For PyCUDA import
+        import pycuda.driver as cuda
+        import pycuda.gpuarray as gpuarray
+        import pycuda.autoinit
+        from pycuda.compiler import SourceModule
+
+    if GpuStyle=='OpenCL':
+        # For PyOpenCL import
+        import pyopencl as cl
+        Id=1
+        for platform in cl.get_platforms():
+            for device in platform.get_devices():
+                deviceType=cl.device_type.to_string(device.type)
+                print "Device #%i of type %s : %s" % (Id,deviceType,device.name)
+                Id=Id+1
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    # La temperature est passee comme parametre, attention au CAST !
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep).astype(numpy.float32)
+
+    E=[]
+    M=[]
+
+    sigma={}
+    for T in Trange:
+        sigma[T]=numpy.copy(sigmaIn)
+
+    jobs=len(Trange)
+    
+    # For GPU, all process are launched
+    if GpuStyle=='CUDA':
+        duration=MetropolisAllCuda(sigma,Trange,J,B,Iterations,jobs,ParaStyle,Alu,Device)
+    else:
+        duration=MetropolisAllOpenCL(sigma,Trange,J,B,Iterations,jobs,ParaStyle,Alu,Device)
+
+    print BestThreadsNumber(len(Trange))
+    
+    for T in Trange:
+        E=numpy.append(E,Energy(sigma[T],J))
+        M=numpy.append(M,Magnetization(sigma[T],B))
+        print "CPU Time for each : %f" % (duration/len(Trange))
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])        
+        ImageOutput(sigma[T],"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    

Property changes on: Ising/GPU/Ising2D-GPU.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/GPU/Ising2D-GPU.py-130528
===================================================================
--- Ising/GPU/Ising2D-GPU.py-130528	(revision 0)
+++ Ising/GPU/Ising2D-GPU.py-130528	(revision 18)
@@ -0,0 +1,601 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+import math
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+from numpy.random import randint as nprnd
+
+KERNEL_CODE_OPENCL="""
+
+// Marsaglia RNG very simple implementation
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+
+#define MWCfp MWC * 2.328306435454494e-10f
+#define KISSfp KISS * 2.328306435454494e-10f
+
+__kernel void MainLoopOne(__global char *s,float T,float J,float B,
+                          uint sizex,uint sizey,
+                          uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z;
+   uint w=seed_w;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*y];
+
+      int d=s[x+sizex*((y+1)%sizey)];
+      int u=s[x+sizex*((y-1)%sizey)];
+      int l=s[((x-1)%sizex)+sizex*y];
+      int r=s[((x+1)%sizex)+sizex*y];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey)] = (char)factor*p;
+   }
+   barrier(CLK_GLOBAL_MEM_FENCE);      
+}
+
+__kernel void MainLoopGlobal(__global char *s,__global float *T,float J,float B,
+                             uint sizex,uint sizey,
+                             uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z/(get_global_id(0)+1);
+   uint w=seed_w/(get_global_id(0)+1);
+   float t;
+   uint ind=get_global_id(0);
+
+   t=T[get_global_id(0)];
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+      
+   }
+
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+
+__kernel void MainLoopLocal(__global char *s,__global float *T,float J,float B,
+                            uint sizex,uint sizey,
+                            uint iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(get_local_id(0)+1);
+   uint w=seed_w/(get_local_id(0)+1);
+   //float t=T[get_local_id(0)+get_global_id(0)*get_local_size(0)];
+   //uint ind=get_local_id(0)+get_global_id(0)*get_local_size(0);
+   float t=T[get_local_id(0)];
+   uint ind=get_local_id(0);
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+   }
+
+   barrier(CLK_LOCAL_MEM_FENCE);
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+"""
+
+KERNEL_CODE_CUDA="""
+
+// Marsaglia RNG very simple implementation
+#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
+#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+#define MWC   (znew+wnew)
+#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
+#define CONG  (jcong=69069*jcong+1234567)
+#define KISS  ((MWC^CONG)+SHR3)
+
+#define MWCfp MWC * 2.328306435454494e-10f
+#define KISSfp KISS * 2.328306435454494e-10f
+
+__global__ void MainLoopOne(char *s,float T,float J,float B,
+                            uint sizex,uint sizey,
+                            uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z;
+   uint w=seed_w;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*y];
+
+      int d=s[x+sizex*((y+1)%sizey)];
+      int u=s[x+sizex*((y-1)%sizey)];
+      int l=s[((x-1)%sizex)+sizex*y];
+      int r=s[((x+1)%sizex)+sizex*y];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/T))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey)] = (char)factor*p;
+   }
+   __syncthreads();
+
+}
+
+__global__ void MainLoopGlobal(char *s,float *T,float J,float B,
+                               uint sizex,uint sizey,
+                               uint iterations,uint seed_w,uint seed_z)
+
+{
+   uint z=seed_z/(blockIdx.x+1);
+   uint w=seed_w/(blockIdx.x+1);
+   float t=T[blockIdx.x];
+   uint ind=blockIdx.x;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+      
+   }
+   __syncthreads();
+
+}
+
+__global__ void MainLoopLocal(char *s,float *T,float J,float B,
+                              uint sizex,uint sizey,
+                              uint iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(threadIdx.x+1);
+   uint w=seed_w/(threadIdx.x+1);
+   float t=T[threadIdx.x];
+   uint ind=threadIdx.x;
+
+   for (uint i=0;i<iterations;i++) {
+
+      uint x=(uint)(MWC%sizex) ;
+      uint y=(uint)(MWC%sizey) ;
+
+      int p=s[x+sizex*(y+sizey*ind)];
+
+      int d=s[x+sizex*((y+1)%sizey+sizey*ind)];
+      int u=s[x+sizex*((y-1)%sizey+sizey*ind)];
+      int l=s[((x-1)%sizex)+sizex*(y+sizey*ind)];
+      int r=s[((x+1)%sizex)+sizex*(y+sizey*ind)];
+
+      float DeltaE=2.0f*p*(J*(u+d+l+r)+B);
+
+      int factor=((DeltaE < 0.0f) || (MWCfp < exp(-DeltaE/t))) ? -1:1;
+      s[x%sizex+sizex*(y%sizey+sizey*ind)] = (char)factor*p;
+   }
+   __syncthreads();
+
+}
+"""
+
+
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,T,J,B,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def MetropolisAllOpenCL(sigmaDict,TList,J,B,iterations,jobs,ParaStyle,Alu,Device):
+
+    # sigmaDict & Tlist are NOT respectively array & float
+    # sigmaDict : dict of array for each temperatoire
+    # TList : list of temperatures
+          
+    # Initialisation des variables en les CASTant correctement
+    
+    # Je detecte un peripherique GPU dans la liste des peripheriques
+
+    HasGPU=False
+    Id=1
+    # Primary Device selection based on Device Id
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            deviceType=cl.device_type.to_string(device.type)
+            if Id==Device and not HasGPU:
+                GPU=device
+                print "CPU/GPU selected: ",device.name
+                HasGPU=True
+            Id=Id+1
+    # Default Device selection based on ALU Type
+    for platform in cl.get_platforms():
+        for device in platform.get_devices():
+            deviceType=cl.device_type.to_string(device.type)
+            if deviceType=="GPU" and Alu=="GPU" and not HasGPU:
+                GPU=device
+                print "GPU selected: ",device.name
+                HasGPU=True
+            if deviceType=="CPU" and Alu=="CPU" and not HasGPU:        
+                GPU=device
+                print "CPU selected: ",device.name
+                HasGPU=True
+                                
+    # Je cree le contexte et la queue pour son execution
+    # ctx = cl.create_some_context()
+    ctx = cl.Context([GPU])
+    queue = cl.CommandQueue(ctx,
+                            properties=cl.command_queue_properties.PROFILING_ENABLE)
+    
+    # Je recupere les flag possibles pour les buffers
+    mf = cl.mem_flags
+
+    # Concatenate all sigma in single array
+    sigma=numpy.copy(sigmaDict[TList[0]])
+    for T in TList[1:]:
+        sigma=numpy.concatenate((sigma,sigmaDict[T]),axis=1)
+
+    sigmaCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=sigma)
+    TCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=TList)
+   
+    MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build( \
+        options = "-cl-mad-enable -cl-fast-relaxed-math")
+
+    SizeX,SizeY=sigmaDict[TList[0]].shape
+
+    if ParaStyle=='Blocks':
+        # Call OpenCL kernel
+        # (1,) is Global work size (only 1 work size)
+        # (1,) is local work size
+        # SeedZCL is lattice translated in CL format
+        # SeedWCL is lattice translated in CL format
+        # step is number of iterations
+        CLLaunch=MetropolisCL.MainLoopGlobal(queue,(jobs,),None,
+                                             sigmaCL,
+                                             TCL, 
+                                             numpy.float32(J),
+                                             numpy.float32(B),
+                                             numpy.uint32(SizeX),
+                                             numpy.uint32(SizeY),
+                                             numpy.uint32(iterations),
+                                             numpy.uint32(nprnd(2**31-1)),
+                                             numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+    else:
+        blocks=int(math.sqrt(jobs))
+        # en OpenCL, necessaire de mettre un Global_id identique au local_id
+        CLLaunch=MetropolisCL.MainLoopLocal(queue,(jobs,),(jobs,),
+                                            sigmaCL,
+                                            TCL,
+                                            numpy.float32(J),
+                                            numpy.float32(B),
+                                            numpy.uint32(SizeX),
+                                            numpy.uint32(SizeY),
+                                            numpy.uint32(iterations),
+                                            numpy.uint32(nprnd(2**31-1)),
+                                            numpy.uint32(nprnd(2**31-1)))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+        
+    CLLaunch.wait()
+    cl.enqueue_copy(queue, sigma, sigmaCL).wait()
+    elapsed = 1e-9*(CLLaunch.profile.end - CLLaunch.profile.start)
+    sigmaCL.release()
+
+    results=numpy.split(sigma,len(TList),axis=1)
+    for T in TList:
+        sigmaDict[T]=numpy.copy(results[numpy.nonzero(TList == T)[0][0]])
+
+    return(elapsed)
+
+def MetropolisAllCuda(sigmaDict,TList,J,B,iterations,jobs,ParaStyle,Alu,Device):
+
+    # sigmaDict & Tlist are NOT respectively array & float
+    # sigmaDict : dict of array for each temperatoire
+    # TList : list of temperatures
+          
+    # Avec PyCUDA autoinit, rien a faire !
+
+    mod = SourceModule(KERNEL_CODE_CUDA)
+    
+    MetropolisBlocksCuda=mod.get_function("MainLoopGlobal")
+    MetropolisThreadsCuda=mod.get_function("MainLoopLocal")
+
+    # Concatenate all sigma in single array
+    sigma=numpy.copy(sigmaDict[TList[0]])
+    for T in TList[1:]:
+        sigma=numpy.concatenate((sigma,sigmaDict[T]),axis=1)
+
+    sigmaCU=cuda.InOut(sigma)
+    TCU=cuda.InOut(TList)
+
+    SizeX,SizeY=sigmaDict[TList[0]].shape
+
+    start = pycuda.driver.Event()
+    stop = pycuda.driver.Event()
+    
+    start.record()
+    start.synchronize()
+    if ParaStyle=='Blocks':
+        # Call CUDA kernel
+        # (1,) is Global work size (only 1 work size)
+        # (1,) is local work size
+        # SeedZCL is lattice translated in CL format
+        # SeedWCL is lattice translated in CL format
+        # step is number of iterations
+        MetropolisBlocksCuda(sigmaCU,
+                             TCU, 
+                             numpy.float32(J),
+                             numpy.float32(B),
+                             numpy.uint32(SizeX),
+                             numpy.uint32(SizeY),
+                             numpy.uint32(iterations),
+                             numpy.uint32(nprnd(2**31-1)),
+                             numpy.uint32(nprnd(2**31-1)),
+                             grid=(jobs,1),block=(1,1,1))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+    else:
+        blocks=int(math.sqrt(jobs))
+        MetropolisThreadsCuda(sigmaCU,
+                              TCU,
+                              numpy.float32(J),
+                              numpy.float32(B),
+                              numpy.uint32(SizeX),
+                              numpy.uint32(SizeY),
+                              numpy.uint32(iterations),
+                              numpy.uint32(nprnd(2**31-1)),
+                              numpy.uint32(nprnd(2**31-1)),
+                              grid=(1,1),block=(jobs,1,1))
+        print "%s with %i %s done" % (Alu,jobs,ParaStyle)
+        
+    stop.record()
+    stop.synchronize()
+    elapsed = start.time_till(stop)*1e-3
+
+    results=numpy.split(sigma,len(TList),axis=1)
+    for T in TList:
+        sigmaDict[T]=numpy.copy(results[numpy.nonzero(TList == T)[0][0]])
+
+    return(elapsed)
+
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Alu can be CPU or GPU
+    Alu='CPU'
+    # Id of GPU : 0 
+    Device=0
+    # GPU style can be Cuda (Nvidia implementation) or OpenCL
+    GpuStyle='OpenCL'
+    # Parallel distribution can be on Threads or Blocks
+    ParaStyle='Blocks'
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:a:d:g:t:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep=","alu=","gpustyle=","parastyle="])
+    except getopt.GetoptError:
+        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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -t <Threads/Blocks>' % sys.argv[0]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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) -a <CPU/GPU> -d <DeviceId> -g <CUDA/OpenCL> -t <Threads/Blocks>' % sys.argv[0]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = float(arg)
+        elif opt in ("-p", "--tempstep"):
+            Tstep = float(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+        elif opt in ("-a", "--alu"):
+            Alu = arg
+        elif opt in ("-d", "--device"):
+            Device = int(arg)
+        elif opt in ("-g", "--gpustyle"):
+            GpuStyle = arg
+        elif opt in ("-t", "--parastyle"):
+            ParaStyle = arg
+
+    if Alu=='CPU' and GpuStyle=='CUDA':
+        print "Alu can't be CPU for CUDA, set Alu to GPU"
+        Alu='GPU'
+
+    if ParaStyle not in ('Blocks','Threads','Hybrid'):
+        print "%s not exists, ParaStyle set as Threads !" % ParaStyle
+        ParaStyle='Blocks'
+   
+    print "Compute unit : %s" % Alu
+    print "Device Identification : %s" % Device
+    print "GpuStyle used : %s" % GpuStyle
+    print "Parallel Style used : %s" % ParaStyle
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    if GpuStyle=='CUDA':
+        # For PyCUDA import
+        import pycuda.driver as cuda
+        import pycuda.gpuarray as gpuarray
+        import pycuda.autoinit
+        from pycuda.compiler import SourceModule
+
+    if GpuStyle=='OpenCL':
+        # For PyOpenCL import
+        import pyopencl as cl
+        Id=1
+        for platform in cl.get_platforms():
+            for device in platform.get_devices():
+                deviceType=cl.device_type.to_string(device.type)
+                print "Device #%i of type %s : %s" % (Id,deviceType,device.name)
+                Id=Id+1
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    # La temperature est passee comme parametre, attention au CAST !
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep).astype(numpy.float32)
+
+    E=[]
+    M=[]
+
+    sigma={}
+    for T in Trange:
+        sigma[T]=numpy.copy(sigmaIn)
+
+    jobs=len(Trange)
+    
+    # For GPU, all process are launched
+    if GpuStyle=='CUDA':
+        duration=MetropolisAllCuda(sigma,Trange,J,B,Iterations,jobs,ParaStyle,Alu,Device)
+    else:
+        duration=MetropolisAllOpenCL(sigma,Trange,J,B,Iterations,jobs,ParaStyle,Alu,Device)
+        
+    for T in Trange:
+        E=numpy.append(E,Energy(sigma[T],J))
+        M=numpy.append(M,Magnetization(sigma[T],B))
+        print "CPU Time for each : %f" % (duration/len(Trange))
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])        
+        ImageOutput(sigma[T],"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    

Property changes on: Ising/GPU/Ising2D-GPU.py-130528
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/MP/Ising2D-MP.py
===================================================================
--- Ising/MP/Ising2D-MP.py	(revision 0)
+++ Ising/MP/Ising2D-MP.py	(revision 18)
@@ -0,0 +1,170 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+
+from multiprocessing import Pool
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,J,B,T,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+    # Default Number of Procs Used
+    Procs=4
+
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:u:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep=","units"])
+    except getopt.GetoptError:
+        print '%s -j <Coupling Factor> -b <Magnetic Field> -z <Size of Lattice> -i <Iterations> -s <Minimum Temperature> -e <Maximum Temperature> -p <steP Temperature> -u <Unit of process> -c (Print Curves)' % sys.argv[0]
+        sys.exit(2)
+    
+    for opt, arg in opts:
+        if opt == '-h':
+            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]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = arg
+        elif opt in ("-p", "--tempstep"):
+            Tstep = numpy.uint32(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+        elif opt in ("-u", "--unit"):
+            Procs = int(arg)
+      
+    print "Process Units : %s" % Procs
+    print "Coupling Factor J : %s" % J
+    print "Magnetic Field B :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
+
+    E=[]
+    M=[]
+
+    def MetropolisStrip(T):
+        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
+        sigma=numpy.copy(sigmaIn)
+        duration=Metropolis(sigma,J,B,T,Iterations)
+        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+        print "CPU Time : %f" % (duration)
+        E,M=Energy(sigma,J),Magnetization(sigma,B)
+        print "Total Energy at Temperature %f : %f" % (T,E)
+        print "Total Magnetization at Temperature %f : %f" % (T,M)
+        return(T,E,M)
+
+    pool=Pool(processes=Procs)
+    print "Start on %i processors..." % Procs
+    # Apply MAP to POOL
+    # MAP: distribution of usecases T to be applied to MetropolisStrip 
+    Results=pool.map(MetropolisStrip,Trange)
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    # Save output
+    numpy.savez("Ising2D_MP_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/MP/Ising2D-MP.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Serial/Ising2D-All.py
===================================================================
--- Ising/Serial/Ising2D-All.py	(revision 0)
+++ Ising/Serial/Ising2D-All.py	(revision 18)
@@ -0,0 +1,162 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,J,B,T,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep="])
+    except getopt.GetoptError:
+        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]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = float(arg)
+        elif opt in ("-p", "--tempstep"):
+            Tstep = float(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+      
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
+
+    E=[]
+    M=[]
+
+    sigma={}
+    for T in Trange:
+        sigma[T]=numpy.copy(sigmaIn)
+
+    for T in Trange:
+        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
+        duration=Metropolis(sigma[T],J,B,T,Iterations)
+        E=numpy.append(E,Energy(sigma[T],J))
+        M=numpy.append(M,Magnetization(sigma[T],B))
+        ImageOutput(sigma[T],"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+    
+        print "CPU Time : %f" % (duration)
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    # Save output
+    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/Serial/Ising2D-All.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Serial/Ising2D.py
===================================================================
--- Ising/Serial/Ising2D.py	(revision 0)
+++ Ising/Serial/Ising2D.py	(revision 18)
@@ -0,0 +1,159 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+
+import sys
+import numpy
+from PIL import Image
+from math import exp
+from random import random
+import time
+import getopt
+import matplotlib.pyplot as plt
+
+def ImageOutput(sigma,prefix):
+    Max=sigma.max()
+    Min=sigma.min()
+    
+    # Normalize value as 8bits Integer
+    SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+    image = Image.fromarray(SigmaInt)
+    image.save("%s.jpg" % prefix)
+    
+def Metropolis(sigma,J,B,T,iterations): 
+    start=time.time()
+
+    SizeX,SizeY=sigma.shape
+    
+    for p in xrange(0,iterations):
+        # Random access coordonate
+        X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY)
+        
+        DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
+                                sigma[X,(Y-1)%SizeY]+
+                                sigma[(X-1)%SizeX,Y]+
+                                sigma[(X+1)%SizeX,Y])+B)
+        
+        if DeltaE < 0. or random() < exp(-DeltaE/T):
+            sigma[X,Y]=-sigma[X,Y]
+    duration=time.time()-start
+    return(duration)
+
+def Magnetization(sigma,M):
+    return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0))
+
+def Energy(sigma,J):
+    # Copier et caster 
+    E=numpy.copy(sigma).astype(numpy.float32)
+    
+    # Appel par slice
+    E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
+                                  E[1:-1,:-2]+E[1:-1,2:])
+    
+    # Bien nettoyer la peripherie
+    E[:,0]=0
+    E[:,-1]=0
+    E[0,:]=0
+    E[-1,:]=0
+    
+    Energy=numpy.sum(E)
+
+    return(Energy/(E.shape[0]*E.shape[1]*1.0))
+
+def DisplayCurves(T,E,M,J,B):
+
+    plt.xlabel("Temperature")
+    plt.ylabel("Energy")
+
+    Experience,=plt.plot(T,E,label="Energy") 
+
+    plt.legend()
+    plt.show()
+
+
+if __name__=='__main__':
+
+    # Set defaults values
+    # Coupling factor
+    J=1.
+    # Magnetic Field
+    B=0.
+    # Size of Lattice
+    Size=256
+    # Default Temperatures (start, end, step)
+    Tmin=0.1
+    Tmax=5
+    Tstep=0.1
+    # Default Number of Iterations
+    Iterations=Size*Size
+
+    # Curves is True to print the curves
+    Curves=False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep="])
+    except getopt.GetoptError:
+        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]
+        sys.exit(2)
+    
+ 
+    for opt, arg in opts:
+        if opt == '-h':
+            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]
+            sys.exit()
+        elif opt == '-c':
+            Curves=True
+        elif opt in ("-j", "--coupling"):
+            J = float(arg)
+        elif opt in ("-b", "--magneticfield"):
+            B = float(arg)
+        elif opt in ("-s", "--tempmin"):
+            Tmin = float(arg)
+        elif opt in ("-e", "--tempmax"):
+            Tmax = float(arg)
+        elif opt in ("-p", "--tempstep"):
+            Tstep = float(arg)
+        elif opt in ("-i", "--iterations"):
+            Iterations = int(arg)
+        elif opt in ("-z", "--size"):
+            Size = int(arg)
+      
+    print "Coupling Factor : %s" % J
+    print "Magnetic Field :  %s" % B
+    print "Size of lattice : %s" % Size
+    print "Iterations : %s" % Iterations
+    print "Temperature on start : %s" % Tmin
+    print "Temperature on end : %s" % Tmax
+    print "Temperature step : %s" % Tstep
+
+    LAPIMAGE=False
+
+    sigmaIn=numpy.where(numpy.random.randn(Size,Size)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigmaIn,"Ising2D_Serial_%i_Initial" % (Size))
+
+    Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep)
+
+    E=[]
+    M=[]
+
+    for T in Trange:
+        # Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy !
+        sigma=numpy.copy(sigmaIn)
+        duration=Metropolis(sigma,J,B,T,Iterations)
+        E=numpy.append(E,Energy(sigma,J))
+        M=numpy.append(M,Magnetization(sigma,B))
+        ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
+    
+        print "CPU Time : %f" % (duration)
+        print "Total Energy at Temperature %f : %f" % (T,E[-1])
+        print "Total Magnetization at Temperature %f : %f" % (T,M[-1])
+
+    if Curves:
+        DisplayCurves(Trange,E,M,J,B)
+
+    # Save output
+    numpy.savez("Ising2D_Serial_%i_%.8i" % (Size,Iterations),(Trange,E,M))
+    

Property changes on: Ising/Serial/Ising2D.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Serial/Ising2D-MWC.py
===================================================================
--- Ising/Serial/Ising2D-MWC.py	(revision 0)
+++ Ising/Serial/Ising2D-MWC.py	(revision 18)
@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+#
+# Ising2D model in serial mode
+#
+# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 
+#
+# Thanks to George Marsaglia : http://en.wikipedia.org/wiki/Multiply-with-carry
+
+import sys
+import numpy
+#import pylab
+from PIL import Image
+from math import exp
+from random import random
+import time
+
+z=0
+w=0
+
+def ImageOutput(sigma,prefix):
+        Max=sigma.max()
+        Min=sigma.min()
+
+        # Normalize value as 8bits Integer
+        SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+        image = Image.fromarray(SigmaInt)
+        image.save("%s.jpg" % prefix)
+
+SIZE=256
+
+def Metropolis(sigma,J,T,iterations,seed_w,seed_z):
+	Energy=0
+	
+	global z
+	global w
+
+	z=seed_z
+	w=seed_w
+	
+	for p in range(0,iterations):
+		# Random access coordonate
+		#X,Y=numpy.random.randint(SIZE),numpy.random.randint(SIZE)
+		
+		# def znew():
+		# 	global z
+		# 	z=numpy.uint32(36969*(z&65535)+(z>>16))
+		# 	return(z)
+		# def wnew():
+		# 	global w
+		# 	w=numpy.uint32(18000*(w&65535)+(w>>16))
+		# 	return(w)
+		# def MWC(): return(numpy.uint32((znew()<<16)+wnew()))
+		# def MWCfp(): return((MWC() + 1.0) * 2.328306435454494e-10)
+		
+		def MWC(): 
+			global w
+			global z
+			z=numpy.uint32(36969*(z&65535)+(z>>16))
+			w=numpy.uint32(18000*(w&65535)+(w>>16))
+			return(numpy.uint32((z<<16)+w))
+
+		def MWCfp(): 
+			global w
+			global z
+			z=numpy.uint32(36969*(z&65535)+(z>>16))
+			w=numpy.uint32(18000*(w&65535)+(w>>16))
+			return(((numpy.uint32((z<<16)+w))+1.0)*
+			       2.328306435454494e-10)
+
+		X=numpy.uint32(SIZE*MWCfp())
+		Y=numpy.uint32(SIZE*MWCfp())
+
+		#print X,Y
+		
+		DeltaE=2.*J*sigma[X,Y]*(sigma[X,(Y+1)%SIZE]+
+					sigma[X,(Y-1)%SIZE]+
+					sigma[(X-1)%SIZE,Y]+
+					sigma[(X+1)%SIZE,Y])
+		
+		if DeltaE < 0. or random() < exp(-DeltaE/T):
+			sigma[X,Y]=-sigma[X,Y]
+			Energy+=DeltaE
+			
+	return(Energy)
+
+if __name__=='__main__':
+
+    J=1.
+    T=0.5
+
+    iterations=numpy.uint32(SIZE*SIZE*SIZE)
+    
+    sigma=numpy.where(numpy.random.randn(SIZE,SIZE)>0,1,-1).astype(numpy.int8)
+
+    ImageOutput(sigma,"Ising2D_OpenCL_%i_Initial" % (SIZE))
+
+    seed_w=numpy.uint32(37)
+    seed_z=numpy.uint32(91)
+	
+    i=0
+    step=262144
+    while (step*i < iterations):
+        start=time.time()
+        Metropolis(sigma,J,T,step,seed_w*(i+1),seed_z*((i+1)%SIZE))
+        stop=time.time()
+        elapsed = (stop-start)
+        print "Iteration %i in %f: " % (i,elapsed)
+        ImageOutput(sigma,"Ising2D_OpenCL_%i_%.3i_Final" % (SIZE,i))
+        i=i+1
+    

Property changes on: Ising/Serial/Ising2D-MWC.py
___________________________________________________________________
Added: svn:executable
   + *

Index: Ising/Serial/PostProcess.py
===================================================================
--- Ising/Serial/PostProcess.py	(revision 0)
+++ Ising/Serial/PostProcess.py	(revision 18)
@@ -0,0 +1,10 @@
+import numpy
+import sys
+npzfile=numpy.load(sys.argv[1])
+T=npzfile['arr_0'][0]
+E=npzfile['arr_0'][1]
+ED=npzfile['arr_0'][2]
+import matplotlib.pyplot
+matplotlib.pyplot.plot(T,E)
+#matplotlib.pyplot.plot(T,E,T,ED)
+matplotlib.pyplot.show()
