Merge pull request #30 from DeLaVlag/master

Latest cuda model generator + example

dsl_cuda/README.md → dsl/README.md

@@ -1,14 +1,14 @@
 # TVB CUDA model generation using LEMS format
 This readme describes the usage of the code generation for models defined in LEMS based XML to Cuda (C) format.
 The LEMS format PR has been adopted and altered to match TVB model names. 
-In LEMSCUDA.py the function "cuda_templating(Model+'_CUDA')" will start the code generation.
-It expects a [model+'_CUDA'].xml file to be present in tvb/dsl_cuda/NeuroML/XMLmodels. 
-The generated file will be placed in tvb/simulator/models.
+In LEMSCUDA.py the function "cuda_templating(Model+'_CUDA', 'path/to/XMLmodels')" will start the code generation.
+It expects a [model+'_CUDA'].xml file to be present in ['path/to/XMLmodels']. 
+The generated file will be placed in ['installpath']'/tvb-hpc/dsl/dsl_cuda/CUDAmodels/'.
 The produced filename is a lower cased [model].py which contains a class named [model].
-In the directory TVB_testsuite the files to run the models on the GPU can be found.
-Execute './runthings cuda Modelname' to start the parameter sweep based simulation.
 
     .. moduleauthor:: Michiel. A. van der Vlag <[email protected]>
+    .. moduleauthor:: Marmaduke Woodman <[email protected]>
+    .. moduleauthor:: Sandra Diaz <[email protected]>
 
 # The CUDA memory model specification
 ![](GPUmemindex.png)
@@ -29,7 +29,7 @@ Mako templating
 
 # XML LEMS Definitions 
 Based on http://lems.github.io/LEMS/elements.html but attributes are tuned for TVB CUDA models.
-As an example an XML line and its translation to CUDA are given. 
+As an example an XML line and its translation to CUDA are given below.
 
 * Constants\
 If domain = 'none' no domain range will be added.\
@@ -180,9 +180,12 @@ for (unsigned int j_node = 0; j_node < n_node; j_node++)
 ```
 
 
-# Running
-Place model file in directory and execute cuda_templating('modelname') function. Resulting model will be
-placed in the CUDA model directory
-
-# TODO
-Add CUDA model validation tests.
+# Running an example
+Place an xml model file in directory used for your XML model storage and execute "cuda_templating(Model+'_CUDA',
+'path/to/XMLmodels')" function. The resulting model will be placed in the CUDA model directory.
+The directory 'tvb-hpc/dsl/dsl_cuda/example/' holds an example how to run the model generator and the CUDA model
+on a GPU.
+From this directory, execute './runthings cuda [Modelname]' to start model generation corresponding to an xml file
+and a parameters sweep simulation with the produced model file on a CUDA enabled machine.
+The cuda parameter indicates a cuda simulation is to be started  and the [modelname] paramater is the model 
+that is the target of simulation.

dsl_cuda/NeuroML/CUDAmodels/epileptor.c → dsl/dsl_cuda/CUDAmodels/epileptor.c

@@ -27,47 +27,47 @@ __device__ float wrap_it_PI(float x)
 }
 __device__ float wrap_it_x1(float x1)
 {
-    int x1dim[] = {-2.0, 1.0};
+    float x1dim[] = {-2.0, 1.0};
     if (x1 < x1dim[0]) x1 = x1dim[0];
     else if (x1 > x1dim[1]) x1 = x1dim[1];
 
     return x1;
 }
 __device__ float wrap_it_y1(float y1)
 {
-    int y1dim[] = {-20.0, 2.0};
+    float y1dim[] = {-20.0, 2.0};
     if (y1 < y1dim[0]) y1 = y1dim[0];
     else if (y1 > y1dim[1]) y1 = y1dim[1];
 
     return y1;
 }
 __device__ float wrap_it_z(float z)
 {
-    int zdim[] = {-2.0, 5.0};
+    float zdim[] = {-2.0, 5.0};
     if (z < zdim[0]) z = zdim[0];
     else if (z > zdim[1]) z = zdim[1];
 
     return z;
 }
 __device__ float wrap_it_x2(float x2)
 {
-    int x2dim[] = {-2.0, 0.0};
+    float x2dim[] = {-2.0, 0.0};
     if (x2 < x2dim[0]) x2 = x2dim[0];
     else if (x2 > x2dim[1]) x2 = x2dim[1];
 
     return x2;
 }
 __device__ float wrap_it_y2(float y2)
 {
-    int y2dim[] = {0.0, 2.0};
+    float y2dim[] = {0.0, 2.0};
     if (y2 < y2dim[0]) y2 = y2dim[0];
     else if (y2 > y2dim[1]) y2 = y2dim[1];
 
     return y2;
 }
 __device__ float wrap_it_g(float g)
 {
-    int gdim[] = {-1.0, 1.0};
+    float gdim[] = {-1.0, 1.0};
     if (g < gdim[0]) g = gdim[0];
     else if (g > gdim[1]) g = gdim[1];
 
@@ -96,8 +96,8 @@ __global__ void Epileptor(
 
     // unpack params
     // These are the two parameters which are usually explore in fitting in this model
-    const float global_coupling = params(0);
-    const float global_speed = params(1);
+    const float global_speed = params(0);
+    const float global_coupling = params(1);
 
     // regular constants
     const float a = 1.0;
@@ -131,6 +131,7 @@ __global__ void Epileptor(
     const float rec_speed_dt = 1.0f / global_speed / (dt);
     const float nsig = sqrt(dt) * sqrt(2.0 * 1e-5);
 
+
     // conditional_derived variable declaration
     float ydot0 = 0.0;
     float ydot2 = 0.0;
@@ -140,24 +141,30 @@ __global__ void Epileptor(
     curandState crndst;
     curand_init(id * (blockDim.x * gridDim.x * gridDim.y), 0, 0, &crndst);
 
-    double x1 = 0.0;
-    double y1 = 0.0;
-    double z = 0.0;
-    double x2 = 0.0;
-    double y2 = 0.0;
-    double g = 0.0;
+    float x1 = 0.0;
+    float y1 = 0.0;
+    float z = 0.0;
+    float x2 = 0.0;
+    float y2 = 0.0;
+    float g = 0.0;
 
     //***// This is only initialization of the observable
     for (unsigned int i_node = 0; i_node < n_node; i_node++)
+    {
         tavg(i_node) = 0.0f;
+        if (i_step == 0){
+            state(i_step, i_node) = 0.001;
+        }
+    }
 
     //***// This is the loop over time, should stay always the same
     for (unsigned int t = i_step; t < (i_step + n_step); t++)
     {
     //***// This is the loop over nodes, which also should stay the same
         for (unsigned int i_node = threadIdx.y; i_node < n_node; i_node+=blockDim.y)
         {
-            float coupling = 0.0f;
+            c_pop1 = 0.0f;
+            c_pop2 = 0.0f;
 
             x1 = state((t) % nh, i_node + 0 * n_node);
             y1 = state((t) % nh, i_node + 1 * n_node);
@@ -183,55 +190,57 @@ __global__ void Epileptor(
                 float x1_j = state(((t - dij + nh) % nh), j_node + 0 * n_node);
 
                 // Sum it all together using the coupling function. Kuramoto coupling: (postsyn * presyn) == ((a) * (sin(xj - xi))) 
-                coupling += c_a * sin(x1_j - x1);
+                c_pop1 += wij * c_a * sin(x1_j - x1);
 
             } // j_node */
 
             // rec_n is used for the scaling over nodes
-            c_pop1 = global_coupling * coupling;
-            c_pop2 = g;
-
-            // The conditional variables
+            c_pop1 *= global_coupling * coupling;
+            c_pop2 *= g;
 
             if (x1 < 0.0)
+                // The conditional variables
                 ydot0 = -a * powf(x1, 2) + b * x1;
             else
                 ydot0 = slope - x2 + 0.6 * powf((z - 4),2);
             if (z < 0.0)
+                // The conditional variables
                 ydot2 = - 0.1 * (powf(z, 7));
             else
                 ydot2 = 0;
             if (modification)
+                // The conditional variables
                 h = x0 + 3. / (1. + exp(-(x1 + 0.5) / 0.1));
             else
                 h = 4 * (x1 - x0) + ydot2;
             if (x2 < -0.25)
+                // The conditional variables
                 ydot4 = 0.0;
             else
                 ydot4 = aa * (x2 + 0.25);
             // This is dynamics step and the update in the state of the node
-            x1 = tt * (y1 - z + Iext + Kvf * c_pop1 + ydot0 );
-            y1 = tt * (c - d * powf(x1, 2) - y1);
-            z = tt * (r * (h - z + Ks * c_pop1));
-            x2 = tt * (-y2 + x2 - powf(x2, 3) + Iext2 + bb * g - 0.3 * (z - 3.5) + Kf * c_pop2);
-            y2 = tt * (-y2 + ydot4) / tau;
-            g = tt * (-0.01 * (g - 0.1 * x1) );
+            x1 += dt * (tt * (y1 - z + Iext + Kvf * c_pop1 + ydot0 ));
+            y1 += dt * (tt * (c - d * powf(x1, 2) - y1));
+            z += dt * (tt * (r * (h - z + Ks * c_pop1)));
+            x2 += dt * (tt * (-y2 + x2 - powf(x2, 3) + Iext2 + bb * g - 0.3 * (z - 3.5) + Kf * c_pop2));
+            y2 += dt * (tt * (-y2 + ydot4) / tau);
+            g += dt * (tt * (-0.01 * (g - 0.1 * x1) ));
 
             // Add noise (if noise components are present in model), integrate with stochastic forward euler and wrap it up
-            x1 += dt * (nsig * curand_normal2(&crndst).x + tt * (y1 - z + Iext + Kvf * c_pop1 + ydot0 ));
-            y1 += dt * (nsig * curand_normal2(&crndst).x + tt * (c - d * powf(x1, 2) - y1));
-            z += dt * (nsig * curand_normal2(&crndst).x + tt * (r * (h - z + Ks * c_pop1)));
-            x2 += dt * (nsig * curand_normal2(&crndst).x + tt * (-y2 + x2 - powf(x2, 3) + Iext2 + bb * g - 0.3 * (z - 3.5) + Kf * c_pop2));
-            y2 += dt * (nsig * curand_normal2(&crndst).x + tt * (-y2 + ydot4) / tau);
-            g += dt * (nsig * curand_normal2(&crndst).x + tt * (-0.01 * (g - 0.1 * x1) ));
+            x1 += nsig * curand_normal2(&crndst).x;
+            y1 += nsig * curand_normal2(&crndst).x;
+            z += nsig * curand_normal2(&crndst).x;
+            x2 += nsig * curand_normal2(&crndst).x;
+            y2 += nsig * curand_normal2(&crndst).x;
+            g += nsig * curand_normal2(&crndst).x;
 
             // Wrap it within the limits of the model
-            wrap_it_x1(x1);
-            wrap_it_y1(y1);
-            wrap_it_z(z);
-            wrap_it_x2(x2);
-            wrap_it_y2(y2);
-            wrap_it_g(g);
+            x1 = wrap_it_x1(x1);
+            y1 = wrap_it_y1(y1);
+            z = wrap_it_z(z);
+            x2 = wrap_it_x2(x2);
+            y2 = wrap_it_y2(y2);
+            g = wrap_it_g(g);
 
             // Update the state
             state((t + 1) % nh, i_node + 0 * n_node) = x1;

dsl_cuda/NeuroML/CUDAmodels/kuramoto.c → dsl/dsl_cuda/CUDAmodels/kuramoto.c

@@ -48,8 +48,8 @@ __global__ void Kuramoto(
 
     // unpack params
     // These are the two parameters which are usually explore in fitting in this model
-    const float global_coupling = params(0);
-    const float global_speed = params(1);
+    const float global_speed = params(0);
+    const float global_coupling = params(1);
 
     // regular constants
     const float omega = 60.0 * 2.0 * M_PI_F / 1e3;
@@ -70,19 +70,24 @@ __global__ void Kuramoto(
     curandState crndst;
     curand_init(id * (blockDim.x * gridDim.x * gridDim.y), 0, 0, &crndst);
 
-    double V = 0.0;
+    float V = 0.0;
 
     //***// This is only initialization of the observable
     for (unsigned int i_node = 0; i_node < n_node; i_node++)
+    {
         tavg(i_node) = 0.0f;
+        if (i_step == 0){
+            state(i_step, i_node) = 0.001;
+        }
+    }
 
     //***// This is the loop over time, should stay always the same
     for (unsigned int t = i_step; t < (i_step + n_step); t++)
     {
     //***// This is the loop over nodes, which also should stay the same
         for (unsigned int i_node = threadIdx.y; i_node < n_node; i_node+=blockDim.y)
         {
-            float coupling = 0.0f;
+            c_0 = 0.0f;
 
             V = state((t) % nh, i_node + 0 * n_node);
 
@@ -103,12 +108,12 @@ __global__ void Kuramoto(
                 float V_j = state(((t - dij + nh) % nh), j_node + 0 * n_node);
 
                 // Sum it all together using the coupling function. Kuramoto coupling: (postsyn * presyn) == ((a) * (sin(xj - xi))) 
-                coupling += wij * a * sin(V_j - V);
+                c_0 += wij * a * sin(V_j - V);
 
             } // j_node */
 
             // rec_n is used for the scaling over nodes
-            c_0 = global_coupling * rec_n * coupling;
+            c_0 *= global_coupling * rec_n;
 
             // This is dynamics step and the update in the state of the node
             V += dt * (omega + c_0);
@@ -123,7 +128,9 @@ __global__ void Kuramoto(
             state((t + 1) % nh, i_node + 0 * n_node) = V;
 
             // Update the observable only for the last timestep
+            if (t == (i_step + n_step - 1)){
                 tavg(i_node + 0 * n_node) = sin(V);
+            }
 
             // sync across warps executing nodes for single sim, before going on to next time step
             __syncthreads();