import pycuda.gpuarray as gpuarray

import numpy
from pycuda.curandom import rand as curand
from pycuda.compiler import SourceModule
import pycuda.driver as cuda

try:
    ctx.pop()
    ctx.detach()
except:

    print "No CTX!"

cuda.init()
device = cuda.Device(0)
ctx = device.make_context()
print device.name(), device.compute_capability(),device.total_memory()/1024.**3,"GB"
print "W systemie mamy :",cuda.Device.count(), " urządzenia"


from pycuda.reduction import ReductionKernel

krnl = ReductionKernel(numpy.dtype(numpy.int32), neutral="0",
        reduce_expr="a+b", 
        map_expr="signbit( (powf(x[i]-0.5f,2.0f)+powf(y[i]-0.5f,2.0f))-0.25f )",
        arguments="float *x, float *y")

(free,total)=cuda.mem_get_info()
free,total

free/4

%%time 
print cuda.mem_get_info()[0]/1024**2
N = 800000000
a = curand((N/2,))
b = curand((N/2,))
ctx.synchronize()
print cuda.mem_get_info()[0]/1024**2

%%time 
for i in range(10):
    a1= a.get()

cuda.mem_get_info()[0]/1024.**3

%%time 
pole = float(krnl(a, b).get())/(N/2)
ctx.synchronize()

print "OK"

print 4*pole,"Pi z %d milionów losowań!"%(N/2/1e6)
print np.pi



try:
    ctx.pop()
    ctx.detach()
    print "OK!"
except:
    print "No CTX!"


x,y = a.get()[::100000],b.get()[::100000]
plot(x,y,'.')
c = (x-0.5)**2+(y-0.5)**2<0.25
plot(x[c],y[c],'.')
    



%load_ext cythonmagic

%%cython

cdef extern from "stdlib.h":
    int RAND_MAX

from libc.stdlib cimport rand
from libc.math cimport pow

def cpu_rand():

    cdef double a,b
    cdef long n=0
    for i in range(800000000):
        a = float(rand())/RAND_MAX
        b = float(rand())/RAND_MAX
        if pow(a-0.5,2.0)+pow(b-0.5,2.0)-0.25<0:
            n=n+1
            
    return n

%%time 
cpu_rand()

628317563/800000000.*4