import numpy as np
import pyopencl as cl
import pyopencl.array as cl_array
from numpy.random import randint as nprnd

deviceID = 0
platformID = 0
workGroup=(1,1)

N = 32768
MyData = np.zeros(N, dtype=cl_array.vec.float8)

dev = cl.get_platforms()[platformID].get_devices()[deviceID]

ctx = cl.Context([dev])
queue = cl.CommandQueue(ctx)
mf = cl.mem_flags
clData = cl.Buffer(ctx, mf.READ_WRITE, MyData.nbytes)

MyRoutines = cl.Program(ctx, """
#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
#define SHR3fp SHR3 * 2.328306435454494e-10f
#define CONGfp CONG * 2.328306435454494e-10f

#define LENGTH 1.

#define PI 3.14159265359

#define SMALL_NUM 0.0000001

__kernel void SplutterSpace(__global float8* clData, 
                               uint seed_z,uint seed_w)
{
    int gid = get_global_id(0);
    uint z=seed_z+(uint)gid;
    uint w=seed_w-(uint)gid;

    clData[gid].s01234567 = (float8) (MWCfp,MWCfp,MWCfp,0.,0.,0.,0.,0.);
}

__kernel void ExtendSegment(__global float8* clData, 
                               uint seed_z,uint seed_w)
{
    int gid = get_global_id(0);
    uint z=seed_z+(uint)gid;
    uint w=seed_w-(uint)gid;

    float theta=MWCfp*PI;
    float phi=MWCfp*PI*2.;
    float sinTheta=sin(theta);
    clData[gid].s4=clData[gid].s0+LENGTH*sinTheta*cos(phi);
    clData[gid].s5=clData[gid].s1+LENGTH*sinTheta*sin(phi);
    clData[gid].s6=clData[gid].s2+LENGTH*cos(theta);

}

__kernel void EstimateLength(__global float8* clData,__global float* clSize)
{
    int gid = get_global_id(0);
    
    clSize[gid]=distance(clData[gid].lo,clData[gid].hi);
}

// Get from http://geomalgorithms.com/a07-_distance.html

__kernel void ShortestDistance(__global float8* clData,__global float* clDistance)
{
    int gidx = get_global_id(0);
    int ggsz = get_global_size(0);
    int gidy = get_global_id(1);
    
    float4   u = clData[gidx].hi - clData[gidx].lo;
    float4   v = clData[gidy].hi - clData[gidy].lo;
    float4   w = clData[gidx].lo - clData[gidy].lo;     

    float    a = dot(u,u);         // always >= 0
    float    b = dot(u,v);
    float    c = dot(v,v);         // always >= 0
    float    d = dot(u,w);
    float    e = dot(v,w);
   
    float    D = a*c - b*b;        // always >= 0
    float    sc, sN, sD = D;       // sc = sN / sD, default sD = D >= 0
    float    tc, tN, tD = D;       // tc = tN / tD, default tD = D >= 0

    // compute the line parameters of the two closest points
    if (D < SMALL_NUM) { // the lines are almost parallel
        sN = 0.0;         // force using point P0 on segment S1
        sD = 1.0;         // to prevent possible division by 0.0 later
        tN = e;
        tD = c;
    }
    else {                 // get the closest points on the infinite lines
        sN = (b*e - c*d);
        tN = (a*e - b*d);
        if (sN < 0.0) {        // sc < 0 => the s=0 edge is visible
            sN = 0.0;
            tN = e;
            tD = c;
        }
        else if (sN > sD) {  // sc > 1  => the s=1 edge is visible
            sN = sD;
            tN = e + b;
            tD = c;
        }
    }

    if (tN < 0.0) {            // tc < 0 => the t=0 edge is visible
        tN = 0.0;
        // recompute sc for this edge
        if (-d < 0.0)
            sN = 0.0;
        else if (-d > a)
            sN = sD;
        else {
            sN = -d;
            sD = a;
        }
    }
    else if (tN > tD) {      // tc > 1  => the t=1 edge is visible
        tN = tD;
        // recompute sc for this edge
        if ((-d + b) < 0.0)
            sN = 0;
        else if ((-d + b) > a)
            sN = sD;
        else {
            sN = (-d +  b);
            sD = a;
        }
    }
    // finally do the division to get sc and tc
    sc = (fabs(sN) < SMALL_NUM ? 0.0 : sN / sD);
    tc = (fabs(tN) < SMALL_NUM ? 0.0 : tN / tD);

    // get the difference of the two closest points
    float4   dP = w + (sc * u) - (tc * v);  // =  S1(sc) - S2(tc)

    clDistance[ggsz*gidy+gidx]=length(dP);   // return the closest distance
}

 """).build()

print 'Tous au meme endroit',MyData

MyRoutines.SplutterSpace(queue, (N,1), None, clData, numpy.uint32(nprnd(2**32)),numpy.uint32(nprnd(2**32)))

cl.enqueue_copy(queue, MyData, clData)

print 'Tous distribues',MyData

MyRoutines.ExtendSegment(queue, (N,1), None, clData,numpy.uint32(nprnd(2**32)),numpy.uint32(nprnd(2**32)))

cl.enqueue_copy(queue, MyData, clData)

print 'Tous avec leur extremite',MyData

MySize = np.zeros(len(MyData), dtype=numpy.float32)
clSize = cl.Buffer(ctx, mf.READ_WRITE, MySize.nbytes)

MyRoutines.EstimateLength(queue, (N,1), None, clData, clSize)
cl.enqueue_copy(queue, MySize, clSize)

print 'La distance de chacun avec son extremite',MySize

MyDistance = np.zeros(len(MyData)*len(MyData), dtype=numpy.float32)
clDistance = cl.Buffer(ctx, mf.READ_WRITE, MyDistance.nbytes)

MyRoutines.ShortestDistance(queue, (N,N), None, clData, clDistance)
cl.enqueue_copy(queue, MyDistance, clDistance)

MyDistance=numpy.reshape(MyDistance,(N,N))

print 'La distance de chacun',MyDistance