 }
 } likelihooder;
 
+struct supportpenalty
+{
+double factor;
+
+  double __device__ operator ()(thrust::tuple<int, cufftComplex> val)
+{
+  int x = thrust::get<0>(val) % SIDE;
+  int y = thrust::get<0>(val) / SIDE;
+  if (!((x < SHORT_SIDE / 2 || x >= SIDE - SHORT_SIDE / 2 - (SHORT_SIDE % 2)) &&
+      (y < SHORT_SIDE / 2 || y >= SIDE - SHORT_SIDE / 2 - (SHORT_SIDE % 2)))) {
+    cufftComplex input = thrust::get<1>(val);
+    return factor * (input.x * input.x + input.y * input.y);
+  }
+  return 0;
+}
+} supportpenalizer;
+
 struct gradientgetter {
   double __device__ gradient(int count, int mask, cufftComplex val, cufftComplex gTempl)
   {

		t_mask = thrust::device_pointer_cast(d_mask);
		t_truePattern = thrust::device_pointer_cast(d_truePattern);
		t_pattern = thrust::device_pointer_cast(d_pattern);

		printf("Memset: %d\n", cudaMemset(d_data, 0, sizeof(cufftComplex) * SIDE * SIDE));
		cufftPlan2d(&plan, SIDE, SIDE, CUFFT_Z2Z);

		term += (val2 - lim) * factor;
		val2 = lim;
	}
	term += log(val2) * count - val2;
	//term -= (val2 - count) * (val2 - count);	
	term -= fabs(val.y) * 1;

	
