 
 struct func_exp {
     __host__ __device__ float operator()(const float &x) const { return exp(x); }
+};
+
+struct func_entropy_kernel {
+    __host__ __device__ float operator()(const float &x) const { float val = x*log2(x); return val != val ? 0 : val; }
+};
+
+struct func_pow2 {
+    __host__ __device__ float operator()(const float &x) const { return pow(x,2); }
 };
 
 thrust::device_vector<float> compute_pij(cublasHandle_t &handle, 

    thrust::device_vector<float> pij_output(N*N);
    cublasSafeCall(cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, N, N, &alpha, thrust::raw_pointer_cast(pij_vals.data()), N, 
                               &beta, thrust::raw_pointer_cast(pij_vals.data()), N, thrust::raw_pointer_cast(pij_output.data()), N));
    return pij_output;
}


    float momentum_weight = 0.9f;


    //printarray(ys, N, 2);
    thrust::device_vector<float> qij(N * N);
    thrust::device_vector<float> dist(N * N);
    float eta = 0.10f;
