<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Diagram bifurkacyjny dla równania logistycznego $x \to a x (1-x)$
Równanie logistyczne jest niezwykle prostym równaniem iteracyjnym wykazującym zaskakująco złożone zachowanie. Jego własności są od lat siedemdziesiątych przedmiotem poważnych prac matematycznych. Pomimo tego wciąż wiele własności jest niezbadanych i zachowanie się rozwiązań tego równania jest dostępne tylko do analizy numerycznej.
Poniższy przykład wykorzystuje pyCUDA do szybkiego obliczenia tak zwanego diagramu bifurkacyjnego równania logistycznego. Uzyskanie takiego diagramu wymaga jednoczesnej symulacji wielu równań z różnymi warunkami początkowymi i różnymi parametrami. Jest to idealne zadanie dla komputera równoległego.
Sposób implementacji
Pierwszą implementacja naszego algorytmu będzie zastosowanie szablonu jądra zwanego 
    ElementwiseKernel
Jest to prosty sposób na wykonanie tej samej operacji na dużym wektorze danych.
Step1:  Jądro Elementwise
Zdefiniujemy sobie jądro, które dla wektora stanów początkowych, element po  elementcie wykona iteracje rówania logistycznego. Ponieważ będziemy chcieli wykonać powyższe iteracje dla różnych parametrów $a$, zdefiniujemy nasze jądro tak by brało zarówno wektor wartości paramteru $a$ jak i wektor wartości początkowych. Ponieważ będziemy mieli tą samą wartość parametru $a$ dla wielu wartości początkowych to wykorzystamy użyteczną w tym przypadku funkcję numpy
Step3:  Algorytm z pętlą wewnątrz jądra CUDA
Napiszmy teraz algorytm, który będzie iterował równanie Niter razy wewnątrz jednego wywołania jądra CUDA.
Step4:  Porównanie z wersją CPU
Dla porównania napiszemy prosty program, który oblicza iteracje równania logistycznego na CPU. Zatosujemy język cython, który umożliwia automatyczne skompilowanie funkcji do wydajnego kodu, którego wydajność jest porównywalna z kodem napisanym w języku C lub podobnym.
W wyniku działania programu widzimy, że nasze jądro wykonuje obliczenia znacznie szybciej.
Step5:  Wizualizacja wyników
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pycuda.gpuarray as gpuarray
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 ("a tak wogóle to mamy tu:",cuda.Device.count(), " urządzenia")
Explanation: Diagram bifurkacyjny dla równania logistycznego $x \to a x (1-x)$
Równanie logistyczne jest niezwykle prostym równaniem iteracyjnym wykazującym zaskakująco złożone zachowanie. Jego własności są od lat siedemdziesiątych przedmiotem poważnych prac matematycznych. Pomimo tego wciąż wiele własności jest niezbadanych i zachowanie się rozwiązań tego równania jest dostępne tylko do analizy numerycznej.
Poniższy przykład wykorzystuje pyCUDA do szybkiego obliczenia tak zwanego diagramu bifurkacyjnego równania logistycznego. Uzyskanie takiego diagramu wymaga jednoczesnej symulacji wielu równań z różnymi warunkami początkowymi i różnymi parametrami. Jest to idealne zadanie dla komputera równoległego.
Sposób implementacji
Pierwszą implementacja naszego algorytmu będzie zastosowanie szablonu jądra zwanego 
    ElementwiseKernel
Jest to prosty sposób na wykonanie tej samej operacji na dużym wektorze danych.
End of explanation
import numpy as np 
Nx = 1024
Na = 1024
a = np.linspace(3.255,4,Na).astype(np.float32)
a =  np.repeat(a,Nx)
a_gpu = gpuarray.to_gpu(a)
x_gpu = curand((Na*Nx,))
from pycuda.elementwise import ElementwiseKernel
iterate = ElementwiseKernel(
        "float *a, float *x",
        "x[i] = a[i]*x[i]*(1.0f-x[i])",
        "iterate")
%%time
Niter = 1000
for i in range(Niter):
    iterate(a_gpu,x_gpu)
ctx.synchronize()
a,x = a_gpu.get(),x_gpu.get()
plt.figure(num=1, figsize=(10, 6))
every = 10
plt.plot(a[::every],x[::every],'.',markersize=1)
plt.plot([3.83,3.83],[0,1])
Explanation: Jądro Elementwise
Zdefiniujemy sobie jądro, które dla wektora stanów początkowych, element po  elementcie wykona iteracje rówania logistycznego. Ponieważ będziemy chcieli wykonać powyższe iteracje dla różnych parametrów $a$, zdefiniujemy nasze jądro tak by brało zarówno wektor wartości paramteru $a$ jak i wektor wartości początkowych. Ponieważ będziemy mieli tą samą wartość parametru $a$ dla wielu wartości początkowych to wykorzystamy użyteczną w tym przypadku funkcję numpy:
a =  np.repeat(a,Nx)
End of explanation
import pycuda.gpuarray as gpuarray
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()
mod = SourceModule(
    __global__ void logistic_iterations(float *a,float *x,int Niter)
    {
      
      int idx = threadIdx.x + blockDim.x*blockIdx.x;
      float a_ = a[idx];
      float x_ = x[idx];
      int i;
      for (i=0;i<Niter;i++){
        
          x_ = a_*x_*(1-x_);
      }
      
      x[idx] = x_;
    }
    )
logistic_iterations = mod.get_function("logistic_iterations")
block_size=128
Nx = 10240
Na = 1024*2
blocks = Nx*Na//block_size
a = np.linspace(3.255,4,Na).astype(np.float32)
a = np.repeat(a,Nx)
a_gpu = gpuarray.to_gpu(a)
x_gpu = curand((Na*Nx,))
%%time
logistic_iterations(a_gpu,x_gpu, np.int32(10000),block=(block_size,1,1), grid=(blocks,1,1))
ctx.synchronize()
a,x = a_gpu.get(),x_gpu.get()
plt.figure(num=1, figsize=(9, 8))
every = 100
plt.plot(a[::every],x[::every],'.',markersize=1,alpha=1)
plt.plot([3.83,3.83],[0,1])
H, xedges, yedges = np.histogram2d(a,x,bins=(1024,1024))
plt.figure(num=1, figsize=(10,10))
plt.imshow(1-np.log(H.T+5e-1),origin='lower',cmap='gray')
Explanation: Algorytm z pętlą wewnątrz jądra CUDA
Napiszmy teraz algorytm, który będzie iterował równanie Niter razy wewnątrz jednego wywołania jądra CUDA.
End of explanation
%load_ext Cython
%%cython
def logistic_cpu(double a = 3.56994):
    cdef double x
    
    cdef int i
    x = 0.1
    for i in range(1000*1024*1024):
        x = a*x*(1.0-x)
        
    return x
%%time
logistic_cpu(1.235)
print("OK")
Explanation: Porównanie z wersją CPU
Dla porównania napiszemy prosty program, który oblicza iteracje równania logistycznego na CPU. Zatosujemy język cython, który umożliwia automatyczne skompilowanie funkcji do wydajnego kodu, którego wydajność jest porównywalna z kodem napisanym w języku C lub podobnym.
W wyniku działania programu widzimy, że nasze jądro wykonuje obliczenia znacznie szybciej.
End of explanation
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
a1,a2 = 3,3.56994567
Nx = 1024
Na = 1024
a = np.linspace(a1,a2,Na).astype(np.float32)
a =  np.repeat(a,Nx)
a_gpu = gpuarray.to_gpu(a)
x_gpu = curand((Na*Nx,))
x = x_gpu.get()
fig = plt.figure()
every = 1
Niter = 10000
for i in range(Niter):
    if i%every==0:
        
        plt.cla()
        
        plt.xlim(a1,a2)
        plt.ylim(0,1)
        fig.suptitle("iteracja: %05d"%i)
        plt.plot(a,x,'.',markersize=1)
        plt.savefig("/tmp/%05d.png"%i)
    if i>10:
        every=2
    if i>30:
        every=10
    if i>100:
        every=50 
    if i>1000:
        every=500    
    iterate(a_gpu,x_gpu)
    ctx.synchronize()
    a,x = a_gpu.get(),x_gpu.get()
%%sh 
cd /tmp
time convert   -delay 20   -loop 0 *.png anim_double.gif && rm *.png
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
block_size=128
Nx = 1024*5
Na = 1024*3
blocks = Nx*Na//block_size
nframes = 22
for i,(a1,a2) in enumerate(zip(np.linspace(3,3.77,nframes),np.linspace(4,3.83,nframes))):
    a = np.linspace(a1,a2,Na).astype(np.float32)
    a =  np.repeat(a,Nx)
    a_gpu = gpuarray.to_gpu(a)
    x_gpu = curand((Na*Nx,))
    x = x_gpu.get()
    
    logistic_iterations(a_gpu,x_gpu, np.int32(10000),block=(block_size,1,1), grid=(blocks,1,1))
    ctx.synchronize()
    a,x = a_gpu.get(),x_gpu.get()
    H, xedges, yedges = np.histogram2d(a,x,bins=(np.linspace(a1,a2,1024),np.linspace(0,1,1024)))
    fig, ax = plt.subplots(figsize=[10,7])
    
    ax.imshow(1-np.log(H.T+5e-1),origin='lower',cmap='gray',extent=[a1,a2,0,1])
    #plt.xlim(a1,a2)
    #plt.ylim(0,1)
    ax.set_aspect(7/10*(a2-a1))
    #fig.set_size_inches(8, 5)
    fig.savefig("/tmp/zoom%05d.png"%i)
    plt.close(fig)
%%sh 
cd /tmp
time convert   -delay 30   -loop 0 *.png anim_zoom.gif   && rm *.png
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
block_size=128
Nx = 1024*5
Na = 1024*3
blocks = Nx*Na//block_size
a1,a2 = 1,4
x1,x2 = 0., 1
a = np.linspace(a1,a2,Na).astype(np.float32)
a =  np.repeat(a,Nx)
a_gpu = gpuarray.to_gpu(a)
x_gpu = curand((Na*Nx,))
x = x_gpu.get()
logistic_iterations(a_gpu,x_gpu, np.int32(10000),block=(block_size,1,1), grid=(blocks,1,1))
ctx.synchronize()
a,x = a_gpu.get(),x_gpu.get()
H, xedges, yedges = np.histogram2d(a,x,bins=(np.linspace(a1,a2,1024),np.linspace(x1,x2,1024)))
fig, ax = plt.subplots(figsize=[10,7])
ax.imshow(1-np.log(H.T+5e-1),origin='lower',cmap='gray',extent=[a1,a2,x1,x2])
#plt.xlim(a1,a2)
#plt.ylim(0,1)
ax.set_aspect(7/10*(a2-a1)/(x2-x1))
#fig.set_size_inches(8, 5)
fig.savefig("/tmp/zoom.png")
plt.close(fig)
Explanation: Wizualizacja wyników
End of explanation
<END_TASK>