/*
void MakeBlurZ(int lengthmz, int numz, int numclose, char *barr, int *closezind, float *mtab, float adductmass, int *nztab, float *dataMZ, int *closeind)
{
	unsigned int i, j, k;
	#pragma omp parallel for private (i,j,k), schedule(auto)
	for (i = 0; i<lengthmz; i++)
	{
		for (j = 0; j<numz; j++)
		{
			if (barr[index2D(numz, i, j)] == 1)
			{
				for (k = 0; k<numclose; k++)
				{
					int indz = (int)(j + closezind[k]);
					int ind;
					if (indz < 0 || indz >= numz || (nztab[j] + closezind[k]) == 0) { closeind[index3D(numz, numclose, i, j, k)] = -1; }
					else
					{
						float point = (float)((mtab[index2D(numz, i, j)] + adductmass*(float)(nztab[j] + closezind[k])) / (float)(nztab[j] + closezind[k]));
						if (point<dataMZ[0] || point>dataMZ[lengthmz - 1]) { closeind[index3D(numz, numclose, i, j, k)] = -1; }
						else {
							ind = nearfast(dataMZ, point, lengthmz);
							closeind[index3D(numz, numclose, i, j, k)] = index2D(numz, ind, indz);
						}
					}
				}
			}
			else
			{
				for (k = 0; k<numclose; k++)
				{
					closeind[index3D(numz, numclose, i, j, k)] = -1;
				}
			}
		}
	}
}


void MakeBlurM(int lengthmz, int numz, int numclose, char *barr, int *closemind, float *mtab, float molig, float adductmass, int *nztab, float *dataMZ, int *closeind)
{
	unsigned int i, j, k;
#pragma omp parallel for private (i,j,k), schedule(auto)
	for (i = 0; i<lengthmz; i++)
	{
		for (j = 0; j<numz; j++)
		{
			if (barr[index2D(numz, i, j)] == 1)
			{
				for (k = 0; k<numclose; k++)
				{
					int ind;
					float point = (float)((mtab[index2D(numz, i, j)] + closemind[k] * molig + adductmass*(float)(nztab[j])) / (float)(nztab[j]));
					if (point<dataMZ[0] || point>dataMZ[lengthmz - 1]) { closeind[index3D(numz, numclose, i, j, k)] = -1; }
					else {
						ind = nearfast(dataMZ, point, lengthmz);
						closeind[index3D(numz, numclose, i, j, k)] = index2D(numz, ind, j);
					}
				}
			}
			else
			{
				for (k = 0; k<numclose; k++)
				{
					closeind[index3D(numz, numclose, i, j, k)] = -1;
				}
			}
		}
	}
}*/


//// "Returns" discrete Fourier transform of input array. Unused.
//// Use fftw library instead
//void discretefouriertransform(double* input, double** output, int length) {
//	double pi = 3.1415926535;
//	for (int i = 0; i < length; i++) {
//		output[i][0] = 0.0; // Real term
//		output[i][1] = 0.0; // Imaginary term
//		for (int j = 0; j < length; j++) {
//			double inner = 2.0 * pi / length * i * j;
//			output[i][0] += input[j] * cos(inner);
//			output[i][1] += input[j] * (-1.0) * sin(inner);
//		}
//	}
//}
//
//// "Returns" discrete inverse Fourier transform of input. Unused.
//// Use fftw library instead
//void inversefouriertransform(double** input, double* output, int length) {
//	double pi = 3.1415926535;
//	for (int i = 0; i < length; i++) {
//		output[i] = 0.0;	// Real term
//		double imag = 0.0;	// Imaginary term
//		for (int j = 0; j < length; j++) {	// (ac - bd) + i(ad + bc)
//			double inner = 2.0 * pi / length * i * j;
//			double c = cos(inner);
//			double d = sin(inner);
//			output[i] += (input[j][0] * c) - (input[j][1] * d);
//			imag += (input[j][0] * d) + (input[j][1] * c);
//		}
//		output[i] = output[i] / ((double)length);
//		imag = imag / ((double)length);
//		printf("Imaginary term is %.2f", imag);
//	}
//}

// Gives convolution of functions a and b. Unused.
// Use cconv2fast instead
/*
void cconv2(double* a, double* b, double* c, int length) {
	double** A = (double **) malloc(length * sizeof(double*));	// (a + bi)
	double** B = (double**) malloc(length * sizeof(double*));	// (c + di)
	double** C = (double**) malloc(length * sizeof(double*));
	for (int i = 0; i < length; i++) {
		A[i] = (double*) calloc(2, sizeof(double));
		B[i] = (double*) calloc(2, sizeof(double));
		C[i] = (double*) calloc(2, sizeof(double));
	}

	discretefouriertransform(a, A, length);
	discretefouriertransform(b, B, length);
	// A * B = (ac - bd) + i(ad + bc)
	for (int j = 0; j < length; j++) {
		C[j][0] = (A[j][0] * B[j][0]) - (A[j][1] * B[j][1]);
		C[j][1] = (A[j][0] * B[j][1]) + (A[j][1] * B[j][0]);
	}
	inversefouriertransform(C, c, length);
	for (int k = 0; k < length; k++) {
		if (c[k] < 0) {
			c[k] = c[k] * -1.0;
		}
		free(A[k]);
		free(B[k]);
		free(C[k]);
	}
	free(A);
	free(B);
	free(C);
}
*/


/*

void softargmax_transposed(float *blur, const int lengthmz, const int numz, const float beta)
{
	float *newblur;
	newblur = calloc(lengthmz*numz, sizeof(float));
	memcpy(newblur, blur, lengthmz*numz * sizeof(float));
#pragma omp parallel for schedule(auto)
	for (int j = 0; j < numz; j++)
	{
		//float sum2 = 0;
		//float sum1 = 0;
		float factor = 0;
		float max1 = 0;
		float max2 = 0;
		float min2 = 1000000000000;
		//float min1 = 1000000000000;
		for (int i = 0; i < lengthmz; i++)
		{
			float d = newblur[index2D(numz, i, j)];
			//sum1 += d;
			//if (d < min1) { min1 = d; }
			if (d > max1) { max1 = d; }
			float e = exp(beta * d);
			//if (beta > 0) { e = exp(beta * d); }
			//else { e = pow(d, -beta); }
			if (e < min2) { min2 = e; }
			if (e > max2) { max2 = e; }
			blur[index2D(numz, i, j)] = e;
			//sum2 += e;
		}
		//float min = min2 - min1;
		//if (beta < 0) { min = min2; }
		//float denom = (sum2 - min2*numz);
		//if (denom != 0) { factor = sum1 / denom; };
		float denom = (max2 - min2);
		if (denom != 0) { factor = max1 / denom; };
		//float sum3 = 0;
		if (factor > 0) {
			for (int i = 0; i < lengthmz; i++)
			{
				blur[index2D(numz, i, j)] -= min2;
				blur[index2D(numz, i, j)] *= factor;
				//sum3 += blur[index2D(numz, i, j)];
			}
		}
		else {
			for (int i = 0; i < lengthmz; i++)
			{
				blur[index2D(numz, i, j)] = 0;
			}
		}
		//printf("%f %f\n", sum1, sum3);
	}
	free(newblur);
	return;
}

void point_smoothing_iso(float *blur, const int lengthmz, const int numz, const int width)
{
	float *newblur;
	newblur = calloc(lengthmz*numz, sizeof(float));
	memcpy(newblur, blur, lengthmz*numz * sizeof(float));
	#pragma omp parallel for schedule(auto)
	for (int i = 0; i < lengthmz; i++)
	{
		for (int j = 0; j < numz; j++)
		{
			int low = i - width;
			if (low < 0) { low = 0; }
			int high = i + width + 1;
			if (high > lengthmz) { high = lengthmz; }

			float sum = 0;
			float val = 0;
			for (int k = low; k < high; k++)
			{
				float e= exp(500*newblur[index2D(numz, k, j)]);
				sum += e;
				if (k == i) { val = e; }
			}
			if(sum!=0){ blur[index2D(numz, i, j)] = val / sum; }
			else { blur[index2D(numz, i, j)] = 0; }
		}
	}
	free(newblur);
	return;
}*/

/*
void point_smoothing_sg(float *blur, const int lengthmz, const int numz, const int width)
{
	float *newblur;
	newblur = calloc(lengthmz*numz, sizeof(float));
	memcpy(newblur, blur, lengthmz*numz * sizeof(float));

	float *c;
	int np = width * 2 + 1;
	c = calloc(np, sizeof(float));
	savgol(c, np, width, width, 0, 2);

	//printf("Savgol: ");
	//for (int i = 0; i < np; i++){printf("%f ", c[i]);}
	//printf("\n");

	#pragma omp parallel for schedule(auto)
	for (int i = 0; i < lengthmz; i++)
	{
		for (int j = 0; j < numz; j++)
		{
			int low = i - width;
			if (low < 0) { low = 0; }
			int high = i + width + 1;
			if (high > lengthmz) { high = lengthmz; }

			float sum = 0;
			for (int k = low; k < high; k++)
			{
				int cind = mod((i - k) , np);
				sum += c[cind] * newblur[index2D(numz, k, j)];
			}

			blur[index2D(numz, i, j)] = sum;
		}
	}
	free(newblur);
	return;
}
*/

/*
void MakeBlur(const int lengthmz, const int numz, const int numclose, char *barr, const int *closezind,
	const int *closemind, const float *mtab, const float molig, const float adductmass, const int *nztab,
	const float *dataMZ,int *closeind, const float threshold, const Config config)
{

	  #pragma omp parallel for schedule(auto)
	  for(int i=0;i<lengthmz;i++)
		{
		  for(int j=0;j<numz;j++)
			{
				if(barr[index2D(numz,i,j)]==1)
				{
					//Reset the threshold if it is zero
					float newthreshold = threshold;
					if (newthreshold == 0) { newthreshold = config.massbins * 3 / (float)nztab[j];}

					int num = 0;
					for(int k=0;k<numclose;k++)
					{
						//Find the z index and test if it is within the charge range
						int indz=(int)(j+closezind[k]);
						if (indz < 0 || indz >= numz || (nztab[j] + closezind[k])==0) { closeind[index3D(numz, numclose, i, j, k)] = -1; }
						else
						{
							//Find the nearest m/z value and test if it is close enough to the predicted one and within appropriate ranges
							float point=(float)((mtab[index2D(numz, i, j)]+closemind[k]*molig+adductmass*(float)(nztab[j]+closezind[k]))/(float)(nztab[j]+closezind[k]));
							if(point<dataMZ[0]- newthreshold || point>dataMZ[lengthmz-1]+ newthreshold){ closeind[index3D(numz, numclose, i, j, k)] = -1; }
							else {
								int ind = nearfast(dataMZ, point, lengthmz);
								float closepoint = dataMZ[ind];
								int newind= index2D(numz, ind, indz);
								if (barr[newind] == 1 && fabs(point- closepoint)< newthreshold) {
									closeind[index3D(numz, numclose, i, j, k)] = newind;
									num += 1;
								}
								else { closeind[index3D(numz, numclose, i, j, k)] = -1; }
								//printf("%d %d %d %f %f %d\n", i, j, k, point, closepoint, closeind[index3D(numz, numclose, i, j, k)]);
							}

						}
					}
					if (num < 2 && config.isotopemode == 0) { barr[index2D(numz, i, j)] = 0; }// printf("%d %d \n", i, j);}
				}
				else
				{
					for (int k=0;k<numclose;k++)
					{
						closeind[index3D(numz, numclose, i, j, k)] = -1;
					}
				}
			}
		}
}*/

/*
//Charge state smooth using a mean filter of the log
void blur_it_geometric_mean(const int lengthmz,
	const int numz,
	const int numclose,
	const int* __restrict closeind,
	float* __restrict newblur,
	const float* __restrict blur,
	const char* __restrict barr,
	const float zerolog)
{
	float mult = 10;
	if (numclose == 1)
	{
		memcpy(newblur, blur, lengthmz * numz * sizeof(float));
	}
	else {
#		pragma omp parallel for schedule(auto)
		for (int i = 0; i < lengthmz * numz; i++)
		{
			float temp = 0;
			if (barr[i] == 1)
			{
				for (int k = 0; k < numclose; k++)
				{
					float temp2 = 0;
					if (closeind[index2D(numclose, i, k)] != -1)
					{
						temp2 = pow(blur[closeind[index2D(numclose, i, k)]], mult);
					}
					temp += temp2;
				}
				temp = pow(temp, 1/mult)/ (float)numclose;
			}
			newblur[i] = temp;
		}
	}
}*/
/*
//Convolution of neighborhood function with gaussian filter.
void blur_it_hybrid1alt(const int lengthmz,
	const int numz,
	const int zlength,
	const int mlength,
	const int * __restrict closeind,
	const int * __restrict closemind,
	const int * __restrict closezind,
	const float * __restrict mdist,
	const float * __restrict zdist,
	float * __restrict newblur,
	const float * __restrict blur,
	const char * __restrict barr,
	const float* __restrict closearray,
	const float zerolog)
{
	int i, j, k, n;
	int numclose = zlength * mlength;
	if (numclose == 1)
	{
		memcpy(newblur, blur, lengthmz*numz * sizeof(float));
	}
	else {
	#pragma omp parallel for private (i,k,j, n), schedule(auto)
		for (i = 0; i < lengthmz; i++)
		{
			for (j = 0; j < numz; j++)
			{
				float temp = 0;
				if (barr[index2D(numz, i, j)] == 1)
				{
					for (k = 0; k < zlength; k++)
					{
						float temp2 = 0;
						for (n = 0; n < mlength; n++)
						{
							int m = index2D(mlength, k, n);
							if (closeind[index3D(numz, numclose, i, j, m)] != -1)
							{
								temp2 += blur[closeind[index3D(numz, numclose, i, j, m)]] * mdist[n] *closearray[index3D(numz, numclose, i, j, m)];
							}
						}
						if (temp2 > 0) { temp += log(temp2); }
						else { temp += zerolog; }
					}
					temp = exp(temp / (float)zlength);
				}
				newblur[index2D(numz, i, j)] = temp;
			}
		}
	}
}*/

/*
void monotopic_to_average(const int lengthmz, const int numz, float *blur, const char *barr, int isolength,
                          const int *__restrict isotopepos, const float *__restrict isotopeval) {
    float *newblur = NULL;
    int l = lengthmz * numz;
    newblur = (float *) calloc(l, sizeof(float));
    if (newblur) {
        for (int i = 0; i < lengthmz; i++) {
            for (int j = 0; j < numz; j++) {
                if (barr[index2D(numz, i, j)] == 1) {
                    float topval = blur[index2D(numz, i, j)];
                    for (int k = 0; k < isolength; k++) {
                        int pos = isotopepos[index3D(numz, isolength, i, j, k)];
                        float val = isotopeval[index3D(numz, isolength, i, j, k)];
                        newblur[index2D(numz, pos, j)] += topval * (float) val;
                    }
                }
            }
        }
        memcpy(blur, newblur, sizeof(float) * lengthmz * numz);
        free(newblur);
    }
}


float isotopemid(float mass, float *isoparams) {
    float a, b, c;
    a = isoparams[4];
    b = isoparams[5];
    c = isoparams[6];
    return a + b * pow(mass, c);
}

float isotopesig(float mass, float *isoparams) {
    float a, b, c;
    a = isoparams[7];
    b = isoparams[8];
    c = isoparams[9];
    return a + b * pow(mass, c);
}

float isotopealpha(float mass, float *isoparams) {
    float a, b;
    a = isoparams[0];
    b = isoparams[1];
    return a * exp(-mass * b);
}

float isotopebeta(float mass, float *isoparams) {
    float a, b;
    a = isoparams[2];
    b = isoparams[3];
    return a * exp(-mass * b);
}

int setup_isotopes(float *isoparams, int *isotopepos, float *isotopeval, float *mtab, int *ztab, char *barr,
                   float *dataMZ, int lengthmz, int numz) {
    float minmass = 100000000;
    float maxmass = 1;
    int i;
    for (i = 0; i < lengthmz * numz; i++) {
        if (barr[i] == 1) {
            float mass = mtab[i];
            if (mass < minmass) { minmass = mass; }
            if (mass > maxmass) { maxmass = mass; }
        }
    }

    float minmid = isotopemid(minmass, isoparams);
    float minsig = isotopesig(minmass, isoparams);
    float maxmid = isotopemid(maxmass, isoparams);
    float maxsig = isotopesig(maxmass, isoparams);

    int isostart = (int) (minmid - 4 * minsig);
    int isoend = (int) (maxmid + 4 * maxsig);
    if (isostart < 0) { isostart = 0; }
    if (isoend < 4) { isoend = 4; }
    int isolength = isoend - isostart;
    return isolength;
}

void make_isotopes(float *isoparams, int *isotopepos, float *isotopeval, float *mtab, int *ztab, char *barr,
                   float *dataMZ, int lengthmz, int numz) {
    float minmass = 100000000;
    float maxmass = 1;
    int i, j, k;
    for (i = 0; i < lengthmz * numz; i++) {
        if (barr[i] == 1) {
            float mass = mtab[i];
            if (mass < minmass) { minmass = mass; }
            if (mass > maxmass) { maxmass = mass; }
        }
    }
    float massdiff = 1.0026;

    float minmid = isotopemid(minmass, isoparams);
    float minsig = isotopesig(minmass, isoparams);
    float maxmid = isotopemid(maxmass, isoparams);
    float maxsig = isotopesig(maxmass, isoparams);

    int isostart = (int) (minmid - 4 * minsig);
    //int isostart = 0;
    int isoend = (int) (maxmid + 4 * maxsig);
    if (isostart < 0) { isostart = 0; }
    if (isoend < 4) { isoend = 4; }
    int isolength = isoend - isostart;
    float *isorange = NULL;
    int *isoindex = NULL;
    isorange = (float *) calloc(isolength, sizeof(float));
    isoindex = (int *) calloc(isolength, sizeof(int));

    if (isorange && isoindex) {
        for (i = 0; i < isolength; i++) {
            isorange[i] = (isostart + i) * massdiff;
            isoindex[i] = (isostart + i);
        }
#pragma omp parallel for private (i,j,k), schedule(auto)
        for (i = 0; i < lengthmz; i++) {
            for (j = 0; j < numz; j++) {
                if (barr[index2D(numz, i, j)] == 1) {
                    float mz = dataMZ[i];
                    int z = ztab[j];
                    for (k = 0; k < isolength; k++) {
                        float newmz = mz + (isorange[k] / ((float) z));
                        int pos = nearfast(dataMZ, newmz, lengthmz);
                        isotopepos[index3D(numz, isolength, i, j, k)] = pos;
                    }
                }
            }
        }

#pragma omp parallel for private (i,j,k), schedule(auto)
        for (i = 0; i < lengthmz; i++) {
            for (j = 0; j < numz; j++) {
                if (barr[index2D(numz, i, j)] == 1) {
                    float mass = mtab[index2D(numz, i, j)];
                    float mid = isotopemid(mass, isoparams);
                    float sig = isotopesig(mass, isoparams);
                    if (sig == 0) {
                        printf("Error: Sigma Isotope Parameter is 0");
                        exit(102);
                    }
                    float alpha = isotopealpha(mass, isoparams);
                    float amp = (1.0 - alpha) / (sig * 2.50662827);
                    float beta = isotopebeta(mass, isoparams);
                    float tot = 0;
                    for (k = 0; k < isolength; k++) {
                        float e = alpha * exp(-isoindex[k] * beta);
                        float g = amp * exp(-pow(isoindex[k] - mid, 2) / (2 * pow(sig, 2)));
                        float temp = e + g;
                        tot += temp;
                        isotopeval[index3D(numz, isolength, i, j, k)] = temp;
                    }
                    for (k = 0; k < isolength; k++) {
                        if (tot > 0) {
                            isotopeval[index3D(numz, isolength, i, j, k)] =
                                    isotopeval[index3D(numz, isolength, i, j, k)] / tot;
                        }
                    }
                }
            }
        }
        free(isorange);
        free(isoindex);
    }
}

void isotope_dist(float mass, int isolength, int *isoindex, float *isovals, float *isoparams) {
    float mid = isotopemid(mass, isoparams);
    float sig = isotopesig(mass, isoparams);
    if (sig == 0) {
        printf("Error: Sigma Isotope Parameter is 0");
        exit(102);
    }
    float alpha = isotopealpha(mass, isoparams);
    float amp = 1.0 - alpha;
    float beta = isotopebeta(mass, isoparams);
    float tot = 0;
    int k;
    for (k = 0; k < isolength; k++) {
        float e = alpha * exp(-isoindex[k] * beta);
        float g = amp / (sig * 2.50662827) * exp(-pow(isoindex[k] - mid, 2) / (2 * pow(sig, 2)));
        tot += e + g;
        isovals[k] = e + g;
    }
    for (k = 0; k < isolength; k++) {
        if (tot > 0) { isovals[k] = isovals[k] / tot; }
    }
}


void test_isotopes(float mass, float *isoparams) {
    int i;
    float maxmid = isotopemid(mass, isoparams);
    float maxsig = isotopesig(mass, isoparams);

    int isostart = 0;
    int isoend = (int) (maxmid + 4 * maxsig);
    if (isoend < 4) { isoend = 4; }
    int isolength = isoend - isostart;
    float *isorange = NULL;
    int *isoindex = NULL;
    isorange = (float *) calloc(isolength, sizeof(float));
    isoindex = (int *) calloc(isolength, sizeof(int));
    if (isorange && isoindex) {
        for (i = 0; i < isolength; i++) {
            isoindex[i] = (isostart + i);
        }
        isotope_dist(mass, isolength, isoindex, isorange, isoparams);
        printf("Distribution: %f", mass);
        for (i = 0; i < isolength; i++) {
            printf("%f ", isorange[i]);
        }
        printf("\n");
        textvectorprint(isorange, isolength);
        free(isorange);
        free(isoindex);
    }
}

void setup_and_make_isotopes(Config *config, Input *inp) {
    printf("Isotope Mode: %d\n", config->isotopemode);

    config->isolength = setup_isotopes(inp->isoparams, inp->isotopepos, inp->isotopeval, inp->mtab, inp->nztab,
                                       inp->barr, inp->dataMZ, config->lengthmz, config->numz);

    int l = config->isolength * config->lengthmz * config->numz;

    inp->isotopepos = (int *) calloc(l, sizeof(int));
    inp->isotopeval = (float *) calloc(l, sizeof(float));

    make_isotopes(inp->isoparams, inp->isotopepos, inp->isotopeval, inp->mtab, inp->nztab, inp->barr, inp->dataMZ,
                  config->lengthmz, config->numz);

    printf("Isotopes set up, Length: %d\n", config->isolength);
}

int nearfast_test(const float *dataMZ, const float point, const int numdat, float cutoff) {
    int index = nearfast(dataMZ, point, numdat);
    float val = dataMZ[index];
    if (fabs(val - point) < cutoff) {
        return index;
    } else { return -1; }
}


//Second Derivative of a Gaussian
float secderndis(float m, float s, float x) {
    if (s == 0) { return 0; }
    return s / 2 * (2 * exp(-pow(m - x, 2) / s)) / s - (4 * exp(-pow(m - x, 2) / s) * pow(m - x, 2)) / pow(s, 2);
}


//Finds nearest factor of two for optimizing FFTs
int twopow(int num) {
    int n, val;
    n = 0;
    val = 1;
    while (n < 100 && val < num) {
        n++;
        val = pow(2, n);
    }
    return val;
}


//Confidence score
float Confidence(float fit, float data) {
    if (data == 0) { return 0; }
    return clip(1 - fabs(fit - data) / data, 0);
}


void textvectorprint(float *arr, int length) {
    int levels = 20;
    float grad = 1.0 / ((float) levels);
    int i, j;
    printf("\n");
    float max = 0;
    for (i = 0; i < length; i++) {
        if (arr[i] > max) { max = arr[i]; }
    }
    if (max != 0) {
        for (i = 0; i < length; i++) {
            arr[i] = arr[i] / max;
        }
    }
    for (i = 0; i < levels; i++) {
        for (j = 0; j < length; j++) {
            if (arr[j] > grad * (levels - i)) { printf("| "); } else { printf("  "); }
        }
        printf("\n");
    }
}


void deepcopy(float *out, const float *in, int len) {
    for (int i = 0; i < len; i++) {
        out[i] = in[i];
    }
}


void blur_noise(float *noise, int lengthmz) {
    float *temp = NULL;
    temp = (float *) calloc(lengthmz, sizeof(float));
    if (temp) {
        memcpy(temp, noise, sizeof(float) * lengthmz);
        int i, j;
        float filter[5] = {-0.1, -0.4, 1, -0.4, -0.1};
#pragma omp parallel for private (i,j), schedule(auto)
        for (i = 0; i < lengthmz; i++) {
            float val = 0;
            for (j = -2; j <= 2; j++) {
                int k = i + j;
                float newval = 0;
                if (k >= 0 && k < lengthmz) {
                    newval = temp[k];
                }
                if (k < 0) {
                    newval = temp[-k];
                }
                if (k >= lengthmz) {
                    newval = temp[2 * lengthmz - k];
                }
                val += newval * filter[j + 2];
            }
            noise[i] = val;
        }
        free(temp);
    }
}

void KillMass(float killmass, int lengthmz, int numz, char *barr, int *nztab, float adductmass, float *dataMZ,
              float psfun, float mzsig) {
    float thresh;
    if (psfun == 0) { thresh = mzsig * 2.35482; } else { thresh = mzsig; }
    int j, i1, i2, k;
    float testmz;
    //#pragma omp parallel for private (j,testmz,i1,i2), schedule(auto)
    for (j = 0; j < numz; j++) {
        testmz = (killmass + adductmass * nztab[j]) / (float) nztab[j];
        i1 = nearfast(dataMZ, testmz - thresh, lengthmz);
        i2 = nearfast(dataMZ, testmz + thresh, lengthmz);
        for (k = i1; k <= i2; k++) { barr[index2D(numz, k, j)] = 0; }
    }
}


void dd_deconv(float *kernel_y, float *data_y, const int length, float *output) {
    // Create flipped point spread function kernel_star
    float *kernel_star = calloc(length, sizeof(float));
    // Create estimate for solution
    float *estimate = calloc(length, sizeof(float));
    // Allocate arrays for convolutions
    float *conv1 = calloc(length, sizeof(float));
    float *conv2 = calloc(length, sizeof(float));

    if (conv1 && conv2 && estimate && kernel_star) {
        for (int i = 0; i < length; i++) {
            kernel_star[i] = kernel_y[length - i - 1];
        }

        for (int i = 0; i < length; i++) {
            estimate[i] = data_y[i];
        }
        // Perform iterations
        int j = 0;
        float diff = 1.0;
        while (j < 50 && diff > 0.0001) {
            // Thresholds same as in Python
            // Find new estimate
            cconv2fast(kernel_y, estimate, conv1, length);
            for (int k = 0; k < length; k++) {
                if (conv1[k] != 0) {
                    conv1[k] = data_y[k] / conv1[k];
                }
            }
            cconv2fast(conv1, kernel_star, conv2, length);
            for (int k = 0; k < length; k++) {
                conv2[k] = conv2[k] * estimate[k]; // Store new estimate in conv2
            }
            // Find how much the estimate changed
            float sum_diff = 0.0;
            float sum_est = 0.0;
            for (int k = 0; k < length; k++) {
                sum_diff += pow((estimate[k] - conv2[k]), 2.0);
                sum_est += estimate[k];
            }
            diff = sum_diff / sum_est;
            // Set new estimate as estimate
            for (int k = 0; k < length; k++) {
                estimate[k] = conv2[k];
            }
            j++;
        }

        // estimate now contains our deconvolution
        // Normalize and "return"
        float estimate_max = 0.0;
        for (int i = 0; i < length; i++) {
            if (estimate_max < estimate[i]) {
                estimate_max = estimate[i];
            }
        }
        for (int i = 0; i < length; i++) {
            output[i] = estimate[i] / estimate_max;
        }

        // Free arrays
        free(kernel_star);
        free(estimate);
        free(conv1);
        free(conv2);
    }
}


void softargmax_transposed(float *blur, const int lengthmz, const int numz, const float beta, const char *barr,
                           const int maxlength, const int speedyflag
                           , int *starttab, int *endtab, float *mzdist, const float mzsig) {
    float *newblur, *deltas, *deltas2, *denom, *denom2;
    int newl = lengthmz * numz;
    newblur = (float *) calloc(newl, sizeof(float));
    deltas = (float *) calloc(lengthmz, sizeof(float));
    deltas2 = (float *) calloc(lengthmz, sizeof(float));
    denom = (float *) calloc(lengthmz, sizeof(float));
    denom2 = (float *) calloc(lengthmz, sizeof(float));
    size_t newl2 = (size_t) lengthmz * numz * sizeof(float);
    memcpy(newblur, blur, newl2);

    //Sum deltas
    sum_deltas(lengthmz, numz, blur, barr, deltas);

    if (mzsig != 0) //Convolve with peak shape
    {
        convolve_simp(lengthmz, maxlength, starttab, endtab, mzdist, deltas, denom, speedyflag);
    } else { memcpy(denom, deltas, sizeof(float) * lengthmz); }

#pragma omp parallel for schedule(auto)
    for (int i = 0; i < lengthmz * numz; i++) {
        blur[i] = exp(beta * newblur[i]) - 1;
    }

    //Sum deltas
    sum_deltas(lengthmz, numz, blur, barr, deltas2);

    if (mzsig != 0) //Convolve with peak shape
    {
        convolve_simp(lengthmz, maxlength, starttab, endtab, mzdist, deltas2, denom2, speedyflag);
    } else { memcpy(denom2, deltas, sizeof(float) * lengthmz); }

#pragma omp parallel for schedule(auto)
    for (int i = 0; i < lengthmz; i++) {
        float factor = 0;
        if (denom2[i] != 0) { factor = denom[i] / denom2[i]; };
        for (int j = 0; j < numz; j++) {
            //blur[index2D(numz, i, j)] -= 1;
            blur[index2D(numz, i, j)] *= factor;
        }
        //printf("%f %f\n", sum1, sum3);
    }
    free(newblur);
    free(deltas);
    free(deltas2);
    free(denom);
    free(denom2);
    return;
}


void softargmax_everything(float *blur, const int lengthmz, const int numz, const float beta) {
    float *newblur;
    int l = lengthmz * numz;
    size_t sizel = (size_t) l * sizeof(float);
    newblur = (float *) calloc(l, sizeof(float));
    if (newblur) {
        memcpy(newblur, blur, sizel);
        //float max1 = 0;
        //float max2 = 0;
        float sum2 = 0;
        float sum1 = 0;
        float min2 = 1000000000000000.0;
        //#pragma omp parallel for schedule(auto), reduction(min:min2), reduction(+:sum1), reduction(+:sum2)
        for (int i = 0; i < lengthmz * numz; i++) {
            float d = newblur[i];
            //if (d > max1) { max1 = d; }
            float e = exp(beta * d);
            //float e = pow(d, -beta);
            //if (e > max2) { max2 = e; }
            blur[i] = e;
            if (e < min2) { min2 = e; }
            sum1 += d;
            sum2 += e;
        }
        float factor = 0;
        //float denom = (max2 - min2);
        //if (denom != 0) { factor = max1 / denom; };
        float denom = (sum2 - min2 * numz * lengthmz);
        if (denom != 0) { factor = sum1 / denom; };
        //if (sum2 != 0) { factor = sum1 / sum2; };
        if (factor > 0) {
#pragma omp parallel for schedule(auto)
            for (int i = 0; i < lengthmz * numz; i++) {
                blur[i] -= min2;
                blur[i] *= factor;
            }
        }
        free(newblur);
    }
    return;
}


void point_smoothing_peak_width(const int lengthmz, const int numz, const int maxlength, const int *starttab,
                                const int *endtab, const float *mzdist, float *blur, const int speedyflag,
                                const char *barr) {
    float *newblur;
    int ln = lengthmz * numz;
    newblur = (float *) calloc(ln, sizeof(float));
    if (newblur) {
        memcpy(newblur, blur, ln * sizeof(float));
        Reconvolve(lengthmz, numz, maxlength, starttab, endtab, mzdist, newblur, blur, speedyflag, barr);
    }
}


*/

	/*
	//Create a list of start and end values to box in arrays based on the above threshold
	starttab = calloc(config.lengthmz, sizeof(int));
	endtab = calloc(config.lengthmz, sizeof(int));
	int maxlength = 1;
	if (config.mzsig != 0) {
		//Gets maxlength and sets start and endtab
		maxlength = SetStartsEnds(config, &inp, starttab, endtab);

		//Changes dimensions of the peak shape function. 1D for speedy and 2D otherwise
		int pslen = config.lengthmz;
		if (config.speedyflag == 0) { pslen = config.lengthmz * maxlength; }
		if (verbose == 1) { printf("Maxlength: %d \t Total Size: %zu GB\n", maxlength, pslen*sizeof(float)/1000000000); }
		mzdist = calloc(pslen, sizeof(float));
		if (pslen * sizeof(float) / 1000000000 > 4) { printf("Danger: Your data may crash the memory. Consider setting the Peak FWHM to 0.\n"); }
		int makereverse = 0;
		if (config.mzsig < 0 || config.beta < 0) { makereverse = 1; rmzdist = calloc(pslen, sizeof(float));
		//memset(rmzdist, 0, pslen * sizeof(float));
		}
		else { rmzdist = calloc(0, sizeof(float)); }

		//Calculates the distance between mz values as a 2D or 3D matrix

		if (config.speedyflag == 0)
		{
			if (verbose == 1) { printf("Making Peak Shape 2D\n"); }
			MakePeakShape2D(config.lengthmz, maxlength, starttab, endtab, inp, fabs(config.mzsig) * config.peakshapeinflate, config.psfun, config.speedyflag, mzdist, rmzdist, makereverse);
		}
		else
		{
			if (verbose == 1) { printf("Making Peak Shape 1D\n"); }
			//Calculates peak shape as a 1D list centered at the first element for circular convolutions
			MakePeakShape1D(inp.dataMZ, config.psmzthresh, config.lengthmz, config.speedyflag, fabs(config.mzsig) * config.peakshapeinflate, config.psfun, mzdist, rmzdist, makereverse);
		}
		if (silent == 0) { printf("mzdist set: %f\t maxlength: %d\n", mzdist[0], maxlength); }

	}
	else
	{
		mzdist = calloc(0, sizeof(float));
		maxlength = 0;
	}*/

	/*
    void mh5readfile2d_grid(hid_t file_id, char* dataname, int length1, int length2, float* data1)
    {
    	//Unfinished
    	const hsize_t length[2] = { length1 ,length2 };
    	if (H5LTpath_valid(file_id, dataname, 1)) {
    		H5Ldelete(file_id, dataname, H5P_DEFAULT);
    	}
    	H5LTmake_dataset_float(file_id, dataname, 2, length, data1);
    }*/

    /*
    void ReadPeaks(const Config config, const char * dataset, float* peakx, float* peaky, float* dscores)
    {
    	//Unfinished
    	char outdat[1024];
    	strjoin(dataset, "/peaks", outdat);
    	float* ptemp = NULL;
    	//ptemp = calloc(decon->plen * 3, sizeof(float));
    }*/


// peak_extracts() extracts for each row in /peakdata and writes to /extracts

/*
int peak_combine(int argc, char* argv[], float* peakx, float* peaky, Config config)
{
	//Unfinished
	char dataset[1024];
	char outdat[1024];
	char strval[1024];

	config.file_id = H5Fopen(argv[1], H5F_ACC_RDWR, H5P_DEFAULT);

	int plen = 1;
	float* peakx = NULL;
	float* peaky = NULL;
	float* dscores = NULL;

	peakx = calloc(plen, sizeof(float));
	peaky = calloc(plen, sizeof(float));
	dscores = calloc(plen, sizeof(float));

	float* numerators = NULL;
	float* denominators = NULL;

	//peakx = realloc(peakx, plen * sizeof(float));
	//peaky = realloc(peaky, plen * sizeof(float));

	//numerators = calloc(plen, sizeof(float));
	//denominators = calloc(plen, sizeof(float));

	//peak_norm(peaky, plen, config.peaknorm);

	int num = 0;
	num = int_attr(config.file_id, "/ms_dataset", "num", num);

	// Get the dscores by scoring the peak at each scan
	for (int i = 0; i < num; i++) {
		//Create a temp array
		float* tempdscores = NULL;
		tempdscores = calloc(plen, sizeof(float));

		//Import everything
		config.metamode = i;
		Decon decon = SetupDecon();
		Input inp = SetupInputs();
		ReadInputs(argc, argv, &config, &inp);
		int status = ReadDecon(&config, inp, &decon);

		//Check to make sure key deconvolution grids are there
		if (status == 1) {
			//Get the scores for each peak
			//score(config, &decon, inp, 0);
			score_from_peaks(plen, peakx, peaky, tempdscores, config, &decon, inp, 0);

			//Average In Dscores
			for (int j = 0; j < plen; j++)
			{
				//Get the peaky value locally
				int index = nearfast(decon.massaxis, peakx[j], decon.mlen);
				float yval = decon.massaxisval[index];
				numerators[j] += yval * tempdscores[j];
				denominators[j] += yval;
			}
		}
		else
		{
			printf("Missing deconvolution outputs. Turn off Fast Profile/Fast Centroid and try deconvolving again.");
		}

		//Free things
		FreeDecon(decon);
		FreeInputs(inp);
		free(tempdscores);
	}

	//Average In Dscores
	for (int j = 0; j < plen; j++)
	{
		if (denominators[j] != 0) { dscores[j] = numerators[j] / denominators[j]; }
	}

	//Write the outputs
	strcpy(dataset, "/peaks");
	makegroup(config.file_id, dataset);
	strjoin(dataset, "/peakdata", outdat);
	printf("\tWriting %d Peaks to: %s\n", plen, outdat);

	float* ptemp = NULL;
	ptemp = calloc(plen * 3, sizeof(float));

	for (int i = 0; i < plen; i++) {
		ptemp[i * 3] = peakx[i];
		ptemp[i * 3 + 1] = peaky[i];
		ptemp[i * 3 + 2] = dscores[i];
	}

	mh5writefile2d_grid(config.file_id, outdat, plen, 3, ptemp);
	free(ptemp);
	//mh5writefile2d(config.file_id, outdat, plen, peakx, peaky);

	peak_extracts(config, peakx, config.file_id, "/mass_data", plen, 0);

	free(peakx);
	free(peaky);
	free(massaxis);
	free(masssum);
	free(dscores);
	free(numerators);
	free(denominators);

	H5Fclose(config.file_id);
}*/