text:
/*************************************************************************************************************
Calculate shear Gaussian + non-Gaussian covariance using CosmoLike.
Gaussian, super-sample and connected non-Gaussian contributions can each be switched on and off independently.
This loops over all combinations of z bins, saving each to disk separately.

Provide two command line arguments:
  1. Path to ini file with all the configuration (see the example_input/example.ini).
  2. Index of the first power spectrum, spec1_idx. This will then iterate over all spec2_idx <= spec1_idx.
**************************************************************************************************************/

#include <string.h>

#include <gsl/gsl_deriv.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sf_expint.h>

// CosmoLike includes (std & GSL includes must be above here)
#include <cosmolike_core/theory/structs.c>
#include <cosmolike_core/theory/basics.c>
#include <cosmolike_core/theory/recompute.c>
#include <cosmolike_core/theory/cosmo3D.c>
#include <cosmolike_core/theory/redshift_spline.c>
#include <cosmolike_core/theory/cosmo2D_fourier.c>
#include <cosmolike_core/theory/halo.c>
#include <cosmolike_core/theory/covariances_3D.c>
#include <cosmolike_core/theory/covariances_fourier.c>
#include <covs/init.c>

// Cl covariance params
typedef struct {
    int lmin;
    int lmax;
    char ell_str[200];
    char output_dir[200];
    int do_g;
    int do_ss;
    int do_cng;
} clcov_par;
clcov_par clcov_params = {.lmin = 2, .lmax = -1, .ell_str = "NULL", .do_g = 0, .do_ss = 0, .do_cng = 0};

// Set Cl covariance params from ini file
void set_clcov_parameters(char *paramfile, int output)
{
  char line[256];
  int iline=0;

  FILE* input = fopen(paramfile, "r");
  while(fgets(line, 256, input) != NULL) {
    char name[128], val[128];
    iline++;
    if(line[0] == '#') continue;
    sscanf(line, "%128s : %128s", name, val);

    if(strcmp(name, "lmin") == 0) {
      sscanf(val, "%d", &clcov_params.lmin);
      if(output == 1) printf("lmin %d \n", clcov_params.lmin);
    } else if(strcmp(name, "lmax") == 0) {
      sscanf(val, "%d", &clcov_params.lmax);
      if(output == 1) printf("lmax %d \n", clcov_params.lmax);
    } else if(strcmp(name, "output_dir") == 0) {
      sprintf(clcov_params.output_dir, "%s", val);
      if(output == 1) printf("output_dir %s \n", clcov_params.output_dir);
    } else if(strcmp(name, "c_footprint_file") == 0) {
      sprintf(covparams.C_FOOTPRINT_FILE, "%s", val);
      if(output == 1) printf("c_footprint_file %s \n", covparams.C_FOOTPRINT_FILE);
    } else if(strcmp(name, "ell") == 0) {
      sprintf(clcov_params.ell_str, "%s", val);
      if(output == 1) printf("ell_str %s \n", clcov_params.ell_str);
    } else if(strcmp(name, "do_g") == 0) {
      sscanf(val, "%d", &clcov_params.do_g);
      if(output == 1) printf("do_g %d \n", clcov_params.do_g);
    } else if(strcmp(name, "do_ss") == 0) {
      sscanf(val, "%d", &clcov_params.do_ss);
      if(output == 1) printf("do_ss %d \n", clcov_params.do_ss);
    } else if(strcmp(name, "do_cng") == 0) {
      sscanf(val, "%d", &clcov_params.do_cng);
      if(output == 1) printf("do_cng %d \n", clcov_params.do_cng);
    }
  }
}

// 2D version of numpy savetxt. Use header = "0" for no header
int savetxt_2d(char *path, double **array, int n_row, int n_col, char *header) {

    FILE *file = fopen(path, "w");

    if (strncmp(header, "0", 2) != 0) fprintf(file, "# %s\n", header);

    for (int i = 0; i < n_row; i++) {
      printf("Saving %d / %d ... \r", i + 1, n_row);
      for (int j = 0; j < n_col; j++) {
        fprintf(file, "%.10e", array[i][j]);
        if (j < n_col - 1) {
          fprintf(file, " ");
        }
      }
      fprintf(file, "\n");
    }

    fclose(file);
    return 0;
}

int main(int argc, char** argv)
{

  if (argc != 3){
    fprintf(stderr, "Syntax: %s config_file spec1_idx\n", argv[0]);
    exit(1);
  }

  setbuf(stdout, NULL);

  // Set this to 1 to output details about inputs for diagnostics
  int output = 1;

  // Set cosmo, survey and Cl covariance params from ini file
  char *inifile = argv[1];
  set_cosmological_parameters(inifile, output);
  set_survey_parameters(inifile, output);
  covparams.full_tomo = 1;
  set_clcov_parameters(inifile, output);
  init_source_sample(redshift.shear_REDSHIFT_FILE, tomo.shear_Nbin);
  init_lens_sample(redshift.clustering_REDSHIFT_FILE, tomo.clustering_Nbin);

  // Determine which contributions to include
  char cov_str[9] = "";
  printf("Including Gaussian contribution: ");
  if (clcov_params.do_g) {
    printf("yes \n");
    sprintf(cov_str, "%sg_", cov_str);
  } else printf("no \n");
  printf("Including super-sample contribution: ");
  if (clcov_params.do_ss) {
    printf("yes \n");
    sprintf(cov_str, "%sss_", cov_str);
  } else printf("no \n");
  printf("Including connected non-Gaussian contribution: ");
  if (clcov_params.do_cng) {
    printf("yes \n");
    sprintf(cov_str, "%scng_", cov_str);
  } else printf("no \n");
  covparams.ssc = clcov_params.do_ss;
  covparams.cng = clcov_params.do_cng;

  // Obtain the number of z bins from config
  int n_zbin;
  if (tomo.shear_Nbin != tomo.clustering_Nbin) error("tomo.shear_Nbin != tomo.clustering_Nbin");
  else n_zbin = tomo.shear_Nbin;

  // Determine ells, either from the list provided in the config or the provided lmin and lmax
  // Do this in 2 steps: first determine the number of ells, then the ells themselves
  int n_ell;
  if (strcmp(clcov_params.ell_str, "NULL") == 0) {
    if (clcov_params.lmax < clcov_params.lmin) error("lmax < lmin, make sure lmax is provided");
    n_ell = clcov_params.lmax - clcov_params.lmin + 1;
  }
  else { // Count the number of commas + 1
    n_ell = 1;
    const char *remaining_str = clcov_params.ell_str;
    while ((remaining_str = strstr(remaining_str, ",")))
    {
      n_ell++;
      remaining_str++;
    }
  }
  int ell[n_ell];
  if (strcmp(clcov_params.ell_str, "NULL") == 0) {
    for (int l = clcov_params.lmin; l <= clcov_params.lmax; l++) {
      ell[l - clcov_params.lmin] = l;
    }
  }
  else {
    char* rem_str = calloc(strlen(clcov_params.ell_str) + 1, sizeof(char));
    strcpy(rem_str, clcov_params.ell_str);
    int l_idx = 0;
    char *token = strtok(rem_str, ",");
    while (token != NULL) {
      ell[l_idx] = strtol(token, NULL, 10);
      token = strtok(NULL, ",");
      l_idx++;
    }
  }

  // Prepare a string about the ells to go in the header of the output text files
  char ell_header_str[200];
  if (strcmp(clcov_params.ell_str, "NULL") == 0) {
    sprintf(ell_header_str, "lmin %d, lmax %d", clcov_params.lmin, clcov_params.lmax);
  } else {
    sprintf(ell_header_str, "ell: %s", clcov_params.ell_str);
  }
  printf("%s \n", ell_header_str);

  // Generate list of power spectra within the shear data vector
  // Rows are the different spectra, cols are: z1, z2
  int n_spectra_sametype = n_zbin * (n_zbin + 1) / 2;
  int spectra_sametype[n_spectra_sametype][2];
  int spec_idx = 0;
  for (int diag = 0; diag < n_zbin; diag++) {
    for (int row = 0; row < n_zbin - diag; row++) {
      spectra_sametype[spec_idx][0] = row;
      spectra_sametype[spec_idx][1] = row + diag;
      spec_idx++;
    }
  }

  // Read in the provided spec1_idx and check it is allowed
  int spec1_idx = atoi(argv[2]);
  if (spec1_idx < 0 || spec1_idx >= n_spectra_sametype) {
    error("spec1_idx provided is not valid");
  } else {
    printf("Requested spec1_idx is %d \n", spec1_idx);
  }

  // Common variables for the covariance blocks
  int z1, z2, z3, z4, l1, l2;
  double cov;
  char save_path[250];
  char header[250];

  // Shear-shear covariance (symmetric so only do spec2_idx <= spec1_idx)
  for (int spec2_idx = 0; spec2_idx <= spec1_idx; spec2_idx++) {

    // Allocate memory for the block and fill with nans
    double** cov_shear_shear = (double**) malloc(n_ell * sizeof(double*));
    for (int i = 0; i < n_ell; i++) cov_shear_shear[i] = (double*) malloc(n_ell * sizeof(double));
    for (int l1_idx = 0; l1_idx < n_ell; l1_idx++) {
      for (int l2_idx = 0; l2_idx < n_ell; l2_idx++) {
        cov_shear_shear[l1_idx][l2_idx] = NAN;
      }
    }

    // Loop over all combinations of (l1, l2), unless the two spectra are the same in which case only need l2 <= l1
    for (int l1_idx = 0; l1_idx < n_ell; l1_idx++) {
      printf("spec2_idx %d / %d, l1_idx %d / %d \n", spec2_idx, spec1_idx, l1_idx, n_ell);
      int l2_idx_max;
      if (spec2_idx == spec1_idx) {
        l2_idx_max = l1_idx;
      } else {
        l2_idx_max = n_ell - 1;
      }
      for (int l2_idx = 0; l2_idx <= l2_idx_max; l2_idx++) {
        l1 = ell[l1_idx];
        l2 = ell[l2_idx];
        z1 = spectra_sametype[spec1_idx][0];
        z2 = spectra_sametype[spec1_idx][1];
        z3 = spectra_sametype[spec2_idx][0];
        z4 = spectra_sametype[spec2_idx][1];
        cov = 0.0;
        if (clcov_params.do_g) cov += cov_G_shear_shear_tomo(l1, l2, z1, z2, z3, z4);
        if (clcov_params.do_ss || clcov_params.do_cng) cov += cov_NG_shear_shear_tomo(l1, l2, z1, z2, z3, z4);
        cov_shear_shear[l1_idx][l2_idx] = cov;
        if (spec2_idx == spec1_idx) {
          cov_shear_shear[l2_idx][l1_idx] = cov;
        }
      }
    }

    // Save block to disk and free memory
    sprintf(save_path, "%s/cov_%sspec1_%d_spec2_%d.txt", clcov_params.output_dir, cov_str, spec1_idx, spec2_idx);
    sprintf(header, "cov_shear_shear for spec_idx %d and %d with do_g=%d, do_ss=%d, do_cng=%d, %s",
            spec1_idx, spec2_idx, clcov_params.do_g, clcov_params.do_ss, clcov_params.do_cng, ell_header_str);
    if (savetxt_2d(save_path, cov_shear_shear, n_ell, n_ell, header) == 0) printf("\nSaved %s\n", save_path);
    for (int i = 0; i < n_ell; i++) free(cov_shear_shear[i]);
    free(cov_shear_shear);

  } // next spec2_idx for this fixed spec1_idx

  printf("Done \n");
  return 0;
}


text:
#pragma once

#pragma once

#include "core/math/vec3.h"
#include "core/math/vec4.h"
#include "core/container/vector.h"
#include "core/memory/linear_allocator.h"
#include <lapacke.h>

/*
void compute_spline_coeff(uint n, float* coeff, float* y, uint stride, float* h, float constant_h) {
	LinearRegion region(get_temporary_allocator());

	float* DL = TEMPORARY_ARRAY(float, n);
	float* D  = TEMPORARY_ARRAY(float, n+1);
	float* DU = TEMPORARY_ARRAY(float, n);
	float* B = TEMPORARY_ARRAY(float, n + 1);

	for (uint i = 1; i < n; i++) {
		float h0 = h ? h[i - 1] : constant_h;
		float h1 = h ? h[i] : constant_h;

		float a0 = y[(i - 1) * stride];
		float a1 = y[i * stride];
		float a2 = y[(i + 1) * stride];

		DL[i-1] = h0;
		D[i] = 2*(h0 + h1);
		DU[i] = h1;
		B[i] = 3/h1 * (a2 - a1) + 3/h0 * (a1 - a0);
	}

	B[0] = 0;
	B[n] = 0;

	D[0] = 1;
	D[n] = 1;

	DU[0] = 0;
	DL[n - 1] = 1;

	int info = LAPACKE_sgtsv(LAPACK_ROW_MAJOR, n+1, 1, DL, D, DU, B, 1);
	assert(info == 0);

	for (uint i = 0; i < n; i++) {
		float h0 = h ? h[i] : constant_h;

		float a0 = y[i * stride];
		float a1 = y[(i + 1) * stride];

		float c0 = B[i];
		float c1 = B[i + 1];

		coeff[i * 4 * stride + 0] = a0;
		coeff[i * 4 * stride + 1] = 1 / h0 * (a1 - a0) - h0 / 3 * (2 * c0 + c1);
		coeff[i * 4 + stride + 2] = c0;
		coeff[i * 4 + stride + 3] = (c1 - c0) / (3*h0);
	}
}
*/

//Cubic hermitian spline
class Spline {
	struct Cubic {
		vec3 a, b, c, d;
	};

	vector<Cubic> coeff;
	vector<float> lengths;

public:
    float total_length;

	Spline(slice<vec3> points) {
		build(points);
	}

    uint points() {
        return lengths.length;
    }

	void build(slice<vec3> points) {
		uint n = points.length - 1;

		coeff.resize(n+1);
		lengths.resize(n);

		LinearRegion region(get_temporary_allocator());

		vec3* a = TEMPORARY_ARRAY(vec3, n + 1);
		for (int i = 1; i <= n - 1; i++)
			a[i] = 3 * ((points[i + 1] - 2 * points[i] + points[i - 1]));

		float* l = TEMPORARY_ARRAY(float, n+1);
		float* mu = TEMPORARY_ARRAY(float, n + 1);
		vec3* z = TEMPORARY_ARRAY(vec3, n + 1);

		l[0] = l[n] = 1;
		mu[0] = 0;
		z[0] = z[n] = vec3();
		coeff[n].c = vec3();

		for (int i = 1; i <= n - 1; i++)
		{
			l[i] = 4 - mu[i - 1];
			mu[i] = 1 / l[i];
			z[i] = (a[i] - z[i - 1]) / l[i];
		}

		for (int i = 0; i < n+1; i++)
			coeff[i].a = points[i];

		for (int j = n - 1; j >= 0; j--)
		{
			coeff[j].c = z[j] - mu[j] * coeff[j + 1].c;
			coeff[j].d = (1.0f / 3.0f) * (coeff[j + 1].c - coeff[j].c);
			coeff[j].b = points[j + 1] - points[j] - coeff[j].c - coeff[j].d;
		}

		/*
		float* DL = TEMPORARY_ARRAY(float, n);
		float* D = TEMPORARY_ARRAY(float, n + 1);
		float* DU = TEMPORARY_ARRAY(float, n);
		float* B = TEMPORARY_ARRAY(float, n + 1);

		for (uint k = 0; k < 3; k++) {
			for (uint i = 1; i < n; i++) {
				float a0 = points[i - 1][k];
				float a1 = points[i][k];
				float a2 = points[i + 1][k];

				DL[i - 1] = 1.0;
				D[i] = 4.0;
				DU[i] = 1.0;
				B[i] = 3.0*(a2 - a1) + 3.0*(a1 - a0);
			}

			B[0] = 0; B[n] = 0;
			D[0] = 1; D[n] = 1;

			DL[n - 1] = 0;
			DU[0] = 0;

            if (true) {
                char spaces[256];
                for (uint i = 0; i < 256; i++) spaces[i] = ' ';

                for (uint i = 0; i < n+1; i++) {
                    if (i == 0) printf("[ %.2f, %.2f,%.*s][%i] = [%f]\n", D[i], DU[0], (n-1)*6, spaces, i, B[i]);
                    else if (i == n) printf("[%.*s %.2f, %.2f ][%i] = [%f]\n", (i-1)*6, spaces, DL[i-1], D[i], i, B[i]);
                    else printf("[%.*s %.2f, %.2f, %.2f,%.*s][%i] = [%f]\n", (i-1)*6, spaces, DL[i-1], D[i], DU[i], (n-1-i)*6, spaces, i, B[i]);
                }
            }

			int info = LAPACKE_sgtsv(LAPACK_ROW_MAJOR, n + 1, 1, DL, D, DU, B, 1);
			assert(info == 0);

			for (uint i = 0; i < n; i++) {
				float a0 = points[i][k];
				float a1 = points[i + 1][k];

				float c0 = B[i];
				float c1 = B[i + 1];

				coeff[i].a[k] = a0;
				coeff[i].b[k] = (a1 - a0) - 1.0 / 3 * (2 * c0 + c1);
				coeff[i].c[k] = c0;
				coeff[i].d[k] = (c1 - c0) / 3.0;
			}
		}
		*/

        total_length = 0;
		for (uint i = 0; i < n; i++) {
			lengths[i] = integrate(i,1.0);
            total_length += lengths[i];
		}
	}

	vec3 at_time(float t) {
		uint i = min(t, coeff.length-1);
		float x = t - i;
		float xx = x * x;
		float xxx = xx * x;

		vec3 result = coeff[i].a + coeff[i].b * x + coeff[i].c * xx + coeff[i].d * xxx;
		return result;
	}

    vec3 tangent(float t) {
        uint i = min(t, coeff.length - 1);
        float x = t - i;
        float xx = x * x;

        vec3 result = coeff[i].b + 2.0f * coeff[i].c * x + 3.0f * coeff[i].d * xx;
        return normalize(result);
    }

	float arc_length_integrand(uint spline, float t) {
		float tt = t * t;

		vec3 result = coeff[spline].b + 2.0f * coeff[spline].c * t + 3.0f * coeff[spline].d * tt;
		return length(result);
	}

	float integrate(uint spline, float t) {
		uint n = 16;
		float h = t / n;

		float XI0 = arc_length_integrand(spline, t);
		float XI1 = 0;
		float XI2 = 0;

		for (uint i = 0; i < n; i++) {
			float X = i * h;
			if (i % 2 == 0) XI2 += arc_length_integrand(spline, X);
			else XI1 += arc_length_integrand(spline, X);
		}

		float XI = h * (XI0 + 2 * XI2 + 4 * XI1) / 3.0f;
		return XI;
	}

    float const_speed_reparam(float t, float speed) {
        float desired_distance = speed * t;

        float distance = 0.0f;

        uint spline;
        for (spline = 0; spline < lengths.length; spline++) {
            float new_distance = distance + lengths[spline];
            if (new_distance >= desired_distance) break;
			distance = new_distance;
        }

        assert(spline < lengths.length);
        float d = desired_distance - distance;

        float s = 0.5f;
        for (uint i = 0; i < 10; i++) {
            float arc_length = integrate(spline, s);
            float arc_speed = arc_length_integrand(spline, s);
            //printf("S %f, arc length %f, arc speed %f\n", s, arc_length, arc_speed);
            s = s - (arc_length - d) / arc_speed;
        }

        return spline + s;
    }

	vec3 const_speed_at_time(float t, float speed) {
		return at_time(const_speed_reparam(t, speed));
	}
};

ext: h

max_stars_repo_name: CompilerLuke/NextEngine

max_issues_repo_name: CompilerLuke/NextEngine

max_forks_repo_name: CompilerLuke/NextEngine


text:
/**
 *
 * @file core_scasum.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:46 2014
 *
 **/
#include <cblas.h>
#include <math.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_scasum - Computes the sums of the absolute values of elements in a same
 *  row or column.
 *  This function is an auxiliary function to norm computations.
 *
 *******************************************************************************
 *
 * @param[in] storev
 *          Specifies whether the sums are made per column or row.
 *          = PlasmaColumnwise: Computes the sum on each column
 *          = PlasmaRowwise:    Computes the sum on each row
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or lower triangular or general
 *          = PlasmaUpperLower: All matrix A is referenced;
 *          = PlasmaUpper: Upper triangle of A is referenced;
 *          = PlasmaLower: Lower triangle of A is referenced.
 *
 * @param[in] M
 *          M specifies the number of rows of the matrix A. M >= 0.
 *
 * @param[in] N
 *          N specifies the number of columns of the matrix A. N >= 0.
 *
 * @param[in] A
 *          A is a M-by-N matrix.
 *
 * @param[in] lda
 *          The leading dimension of the array A. lda >= max(1,M).
 *
 * @param[out] work
 *          Array of dimension M if storev = PlasmaRowwise; N otherwise.
 *          On exit, contains the sums of the absolute values per column or row.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_scasum = PCORE_scasum
#define CORE_scasum PCORE_scasum
#endif
void CORE_scasum(PLASMA_enum storev, PLASMA_enum uplo, int M, int N,
                 const PLASMA_Complex32_t *A, int lda, float *work)
{
    const PLASMA_Complex32_t *tmpA;
    float *tmpW, sum, abs;
    int i,j;

    switch (uplo) {
    case PlasmaUpper:
        for (j = 0; j < N; j++) {
            tmpA = A+(j*lda);
            sum = 0.0;
            for (i = 0; i < j; i++) {
                abs      = cabsf(*tmpA);
                sum     += abs;
                work[i] += abs;
                tmpA++;
            }
            work[j] += sum + cabsf(*tmpA);
        }
        break;
    case PlasmaLower:
        for (j = 0; j < N; j++) {
            tmpA = A+(j*lda)+j;

            sum = 0.0;
            work[j] += cabsf(*tmpA);

            tmpA++;
            for (i = j+1; i < M; i++) {
                abs      = cabsf(*tmpA);
                sum     += abs;
                work[i] += abs;
                tmpA++;
            }
            work[j] += sum;
        }
        break;
    case PlasmaUpperLower:
    default:
        if (storev == PlasmaColumnwise) {
            for (j = 0; j < N; j++) {
                /* work[j] += cblas_scasum(M, &(A[j*lda]), 1); */
                tmpA = A+(j*lda);
                for (i = 0; i < M; i++) {
                    work[j] +=  cabsf(*tmpA);
                    tmpA++;
                }
            }
        }
        else {
            for (j = 0; j < N; j++) {
                tmpA = A+(j*lda);
                tmpW = work;
                for (i = 0; i < M; i++) {
                    /* work[i] += cabsf( A[j*lda+i] );*/
                    *tmpW += cabsf( *tmpA );
                    tmpA++; tmpW++;
                }
            }
        }
    }
}

ext: c


text:
//
//  main.c
//  ODE_MODEL
//
//  Created by Kari Heine on 28/10/2020.
//  Copyright © 2020 Kari Heine. All rights reserved.
//

#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_eigen.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <Accelerate/Accelerate.h>

#include "particle_filters.h"

#define MAX_TEST_MESH_SIZE 4000 // Mesh size for BPF and MLBPF level 1 (the truth mesh size)
#define N_SETTINGS 500 // Number of different scenarios
#define N_SCALES 1 // Scaling samplesizes to study convergence
#define N_MEASUREMENTS 2 // Number of places where the beam deflection is measured

void generate_signal(double* x, int len, double std_sig, double x0, gsl_rng* rng, double L);
void create_dense_solver_matrix(double* A, int n);
void init_with_zeros(double* array, long length);
void make_lapack_band_matrix(double *A, int n, double *ab, int ku, int kl);
double compute_estimation_error(double* x_hat, double* x_ref, int len);
void particle_allocation_time_matching(int max_mesh_size, int* level_0_mesh_sizes, int n_level_0_mesh_sizes, int bpf_sample_size, double std_sig, double obs_std, int data_length, double x0, double *y, int kl, int ku, double* ab_bpf, int ldab, double w, double d, double E, double I, double h_bpf, double *ab0, double h1, double L, double* meass, int N_meas);
void get_settings(int** settings);
void write_beam_in_file(int data_length, int true_mesh_size, double h, double* beam);

int main(void) {

  gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
  gsl_rng_set(rng,clock());
//  gsl_rng_set(rng,42);

  /* Global parameters */
  /* ----------------- */
  int N_top = 250; // Top leve Monte Carlo iterations
  int kl = 3; // Lower diagonal bandwidth for the ODE solver
  int ku = 3; // upper diagonal bandwidth for the ODE solver
  int ldab = 2 * kl + ku + 1; // parameter for Lapack DGBSV
  int run_time_matching = 0; // This toggle (0 or 1) determines if we want to run time matching or filtering
  int N_ref = 100000; // Sample size for the reference filter

  /* Signal and observations */
  /* ----------------------- */
  int data_length = 50; // number of time steps
  double x0 = 2; // initial state
  double std_sig = 0.02; // signal noise standard deviation
  double obs_std = 0.0002; // observation noise standard deviation

  // Data arrays for the signal and observations
  double* x = (double*) malloc(data_length * sizeof(double));
  double* y = (double*) malloc(N_MEASUREMENTS * data_length * sizeof(double));
  FILE * data_out;

  /* Euler-Bernoulli beam parameters */
  /* ------------------------------- */
  double L = 4; // Length (m)
  double w = 0.3; // Width (m)
  double d = 0.03; // thickness (m)
  double E = 1e10; // Young's modulus
  double I = w * d * d * d / (double) 12; // Inertia
  double measurements[N_MEASUREMENTS] = {1, 1.75}; // deflection measurement points along the beam (m)
  int minds[N_MEASUREMENTS]; // mesh point indices for the deflection measurement points

  /* Beam solver setup for data generation */
  /* ------------------------------------- */
  int true_mesh_size = MAX_TEST_MESH_SIZE;
  double* A = (double*) malloc(MAX_TEST_MESH_SIZE * MAX_TEST_MESH_SIZE * sizeof(double));
  int* ipiv = (int*) malloc(MAX_TEST_MESH_SIZE * sizeof(int));
  double* ab = (double*) malloc(ldab * MAX_TEST_MESH_SIZE * sizeof(double));
  double* ab_copy = (double*) malloc(ldab * MAX_TEST_MESH_SIZE * sizeof(double));
  double* signal_payload = (double*) malloc(MAX_TEST_MESH_SIZE * data_length * sizeof(double));
  double h = L / (double)(true_mesh_size-1);
  int ldb = true_mesh_size, info = 0;

  // Calculate the indices for the deflection measurement points.
  for(int i = 0; i < N_MEASUREMENTS; i++) {
    minds[i] = (int)((measurements[i] / L) * (double) true_mesh_size);
  }

  create_dense_solver_matrix(A, true_mesh_size);

  /* Convert the design matrix to general band representation */
  init_with_zeros(ab, ldab * true_mesh_size);
  make_lapack_band_matrix(A, true_mesh_size, ab, ku, kl);

  /* COMPARISON SETUP */
  /* ---------------- */
  FILE * error_out = fopen("error.txt", "w");
  double bpf_MSE = 0;
  double mlbpf_MSE = 0;
  double scaler = 1;
  double bpf_MSE_error_matched = 0;

  /* BPF SETUP */
  /* --------- */
  double *x_hat_bpf = (double*) malloc(data_length * sizeof(double));
  double *x_ref = (double*) malloc(data_length * sizeof(double));
  double *p_ref = (double*) malloc(data_length * sizeof(double));
  FILE* ref_out;
  double filtering_time[4]; // array for time recording
  int N = 500;
  int mesh_size_bpf = MAX_TEST_MESH_SIZE;
  double* ab_bpf = ab;
  int* ipiv_bpf = ipiv;
  double h_bpf = L / (double) (mesh_size_bpf-1);

  create_dense_solver_matrix(A, mesh_size_bpf);
  init_with_zeros(ab_bpf, ldab*mesh_size_bpf);
  make_lapack_band_matrix(A, mesh_size_bpf, ab_bpf, ku, kl);

  /* MLBPF SETUP */
  /* ----------- */
  double worst_case_sign_ratio[1] = { 100 }; // This is basically sustitute for "Inf"
  double *x_hats_mlbpf = (double*) malloc(data_length * sizeof(double));
  int mesh_sizes[2] = {10, 500}; // these values will be overwritten
  int N0, N1;

  double* ab0 = (double*) malloc(ldab * MAX_TEST_MESH_SIZE * sizeof(double));
  double* ab1 = ab;
  int* ipiv0 = ipiv;
  int* ipiv1 = ipiv;
  double h0;

  init_with_zeros(ab0, ldab * MAX_TEST_MESH_SIZE);
  init_with_zeros(ab1, ldab * MAX_TEST_MESH_SIZE);

  // We can create the solver matrices for level 1 only which has fixed mesh size
  create_dense_solver_matrix(A, MAX_TEST_MESH_SIZE);
  make_lapack_band_matrix(A, MAX_TEST_MESH_SIZE, ab1, ku, kl);
  double h1 = L / (double)(MAX_TEST_MESH_SIZE - 1);

  if(run_time_matching) {
    /*
     * Particle allocation time matching
     *
     * This part of the code should be run when trying to find the correct N0 sample sizes to match the
     * BPF running time
     */
    int max_mesh_size = MAX_TEST_MESH_SIZE;
    int level_0_mesh_sizes[20] = {100, 105, 110, 115, 120, 130, 140, 150, 160, 180, 200, 220, 240, 260, 280, 320, 360, 400, 440, 500};
    int n_level_0_mesh_sizes = 20;

    generate_signal(x, data_length, std_sig, x0, rng, L);
    create_payload(data_length, true_mesh_size, signal_payload, w, d, x, h, E, I, L);

    array_copy(ab, ab_copy, ldab * true_mesh_size); // Restore undamaged band matrix

    // Solve the Euler-Bernoulli beam
    dgbsv_(&true_mesh_size, &kl, &ku, &data_length, ab_copy, &ldab, ipiv, signal_payload, &ldb, &info);

    // Generate measurements
    for(int i = 0; i < data_length; i++)
      for(int j = 0; j < N_MEASUREMENTS; j++)
        y[i + j * data_length] = signal_payload[i * true_mesh_size + minds[j]] + gsl_ran_gaussian(rng, obs_std);

    particle_allocation_time_matching(max_mesh_size, level_0_mesh_sizes, n_level_0_mesh_sizes, N, std_sig, obs_std, data_length, x0, y, kl, ku, ab_bpf, ldab, w, d, E, I, h_bpf, ab0, h1, L, measurements, N_MEASUREMENTS);

    printf("Particle allocation time matching done! Exiting...\n");
    return 0; // the code exits
  }

  /*
   *  THIS IS THE ACTUAL FILTER COMPARISON PART
   */

  // Settings are hard coded here. TODO: This should be changed to reading a file.
  int** settings = (int **)malloc(N_SETTINGS * sizeof(int*));
  for(int i = 0; i < N_SETTINGS; i++)
    settings[i] = (int *) malloc(4*sizeof(int));
  int sample_sizes[2] = {0,0};
  get_settings(settings);

  // Top level iteration
  for(int nmc = 0; nmc < N_top; nmc++) {

    printf("Top level iteration %i out of %i.*****************************\n", nmc + 1, N_top);

    /* To save computational time and also to take into account the possibility that a particular
     data set is more favourable to one or the other algorithm, we create new data every nth top
     level iteration, and compute the reference filter.
     */
    if (nmc % 10 == 0) {

      // Create the signal
      generate_signal(x, data_length, std_sig, x0, rng, L);

      printf("Signal generated\n");
      fflush(stdout);

      create_payload(data_length, true_mesh_size, signal_payload, w, d, x, h, E, I, L);

      printf("Payload created\n");
      fflush(stdout);

      array_copy(ab, ab_copy, ldab * true_mesh_size);

      // Solve the Euler-Bernoulli beam
      dgbsv_(&true_mesh_size, &kl, &ku, &data_length, ab_copy, &ldab, ipiv, signal_payload, &ldb, &info);

      printf("\nBeam solved. Info: %i\n\n",info);
      fflush(stdout);

      // write_beam_in_file(data_length, true_mesh_size, h, signal_payload);

      // Generate the noisy measurements
      data_out = fopen("data.txt","w");
      for(int i = 0; i < data_length; i++) {

        // Noisy measurements for one time step but all measurement points
        for(int j = 0; j < N_MEASUREMENTS; j++) {
          y[i + j*data_length] = signal_payload[i * true_mesh_size + minds[j]] + gsl_ran_gaussian(rng, obs_std);
        }

        fprintf(data_out,"%i %e ",i, x[i]); // write true signal into a file
        for(int j = 0; j < N_MEASUREMENTS; j++) {
          // write noisy and noiseless observation into a file
          fprintf(data_out,"%6.6e %6.6e ", y[i+j*data_length], signal_payload[i * true_mesh_size + minds[j]]);
        }
        fprintf(data_out,"\n");
      }
      fclose(data_out);

      printf("Data generated\n");

      /* Compute a reference solution */
      printf("Reference filter with N_ref = %i\n", N_ref);
      fflush(stdout);
      bootstrapfilter(y, std_sig, obs_std, data_length, x0, N_ref, x_ref, p_ref, filtering_time, mesh_size_bpf, kl, ku, ab_bpf, ldab, ipiv_bpf, mesh_size_bpf, info, w, d, E, I, h_bpf, L, measurements, (int)N_MEASUREMENTS);

      // Write the reference filter output into a file with columns:
      // step index, posterior mean, posterior variance, true signal, observation
      ref_out = fopen("ref_filter_out.txt","w");
      for(int i = 0; i < data_length;i++)
        fprintf(ref_out,"%i %e %e %e %e\n",i,x_ref[i],p_ref[i],x[i],y[i]);
      fclose(ref_out);

      // return 0; // option to exit after reference filter may be useful in debugging

    }

    /* Run the benchmark BPF */
    // NB: This is not the reference, but a comparison filter which is time is matched to MLBPF
    bootstrapfilter(y, std_sig, obs_std, data_length, x0, N, x_hat_bpf, NULL, filtering_time, mesh_size_bpf, kl, ku, ab_bpf, ldab, ipiv_bpf, mesh_size_bpf, info, w, d, E, I, h_bpf, L,  measurements, (int)N_MEASUREMENTS);

    // bpf_MSE = compute_estimation_error(x_hat_bpf,x,data_length); // compute the error w.r.t. the truth
    bpf_MSE = compute_estimation_error(x_hat_bpf,x_ref,data_length);

    // Write the results in a file
    fprintf(error_out, "%e %i %i %f %i %e %e %f %e %f %i\n", bpf_MSE, 0, N, filtering_time[0], 0,         filtering_time[1], filtering_time[2], worst_case_sign_ratio[0], filtering_time[3], scaler, mesh_size_bpf);

    /* Compute the error matched BPF */
    // The ERROR matched filter, i.e. it produces roughly the same level of error, but hopefully takes
    // much longer to achieve it.
    printf("Error matched filter\n");
    bootstrapfilter(y, std_sig, obs_std, data_length, x0, 4*N, x_hat_bpf, NULL, filtering_time, mesh_size_bpf, kl, ku, ab_bpf, ldab, ipiv_bpf, mesh_size_bpf, info, w, d, E, I, h_bpf, L,  measurements, (int)N_MEASUREMENTS);

    //  bpf_MSE = compute_estimation_error(x_hat_bpf,x,data_length); // compute the error w.r.t. the truth
    bpf_MSE_error_matched = compute_estimation_error(x_hat_bpf,x_ref,data_length);

    // Write the results in a file
    fprintf(error_out, "%e %i %i %f %i %e %e %f %e %f %i\n", bpf_MSE_error_matched, 0, N, filtering_time[0], -1,         filtering_time[1], filtering_time[2], worst_case_sign_ratio[0], filtering_time[3], scaler, mesh_size_bpf);

    /* Iterate over different mesh size and particle allocation configurations */
    for(int i_setting = 0; i_setting < N_SETTINGS; i_setting++) {

      // Rename
      sample_sizes[0] = settings[i_setting][0];
      sample_sizes[1] = settings[i_setting][1];
      mesh_sizes[0] = settings[i_setting][2];
      mesh_sizes[1] = settings[i_setting][3];

      // Ensure odd number of particles
      sample_sizes[0] = (sample_sizes[0] + sample_sizes[1]) % 2 == 0 ? sample_sizes[0] + 1 : sample_sizes[0];

      init_with_zeros(ab0, ldab * mesh_sizes[0]);
      create_dense_solver_matrix(A, mesh_sizes[0]);
      make_lapack_band_matrix(A, mesh_sizes[0], ab0, ku, kl);
      h0 = L / (double)(mesh_sizes[0] - 1);

      printf("Setting %i out of %i\n",i_setting+1,N_SETTINGS);
      fflush(stdout);

      MLbootstrapfilter(y, std_sig, obs_std, data_length, x0, sample_sizes, worst_case_sign_ratio, x_hats_mlbpf, filtering_time, mesh_sizes, kl, ku, ab0, ab1, ldab, ipiv0, ipiv1, mesh_sizes[0], mesh_sizes[1], info, w, d, E, I, h0, h1, L, measurements, (int)N_MEASUREMENTS);

      //  mlbpf_MSE = compute_estimation_error(x_hats_mlbpf,x,data_length);
      mlbpf_MSE = compute_estimation_error(x_hats_mlbpf,x_ref,data_length);

      N0 = sample_sizes[0];
      N1 = sample_sizes[1];

      fprintf(error_out, "%e %i %i %f %i %e %e %f %e %f %i\n", mlbpf_MSE, N0, N1, filtering_time[0], i_setting + 1,         filtering_time[1], filtering_time[2], worst_case_sign_ratio[0], filtering_time[3], scaler, mesh_sizes[0]);
      fflush(error_out);

    }
  }

  /* Free all that was allocated */
  free(x);
  free(y);
  free(A);
  free(ipiv);
  free(ab);
  free(ab_copy);
  free(signal_payload);
  free(x_hat_bpf);
  free(x_ref);
  free(p_ref);
  free(x_hats_mlbpf);
  free(ab0);
  for(int i = 0; i < N_SETTINGS; i++)
    free(settings[i]);
  free(settings);

  return 0;
}

void write_beam_in_file(int data_length, int true_mesh_size, double h, double* beam) {

  FILE * payload_out = fopen("payload.txt","w");
  for(int i = 0; i < data_length; i++)
    for(int j = 0; j < true_mesh_size; j++)
      fprintf(payload_out, "%5.5e %5.10e \n", j * h, beam[i * true_mesh_size + j]);
  fclose(payload_out);
  printf("Payload written in file\n");
  fflush(stdout);

}

void particle_allocation_time_matching(int max_mesh_size, int* level_0_mesh_sizes, int n_level_0_mesh_sizes, int bpf_sample_size, double std_sig, double obs_std, int data_length, double x0, double *y, int kl, int ku, double* ab_bpf, int ldab, double w, double d, double E, double I, double h_bpf, double *ab0, double h1, double L, double* meass, int N_meas) {

  int N1_step = 20;
  int N0_step = 10000;
  int N0_init = 10;
  int N_top = 5;
  double time_per_iteration;
  double time_tolerance = 1.00; //
  double filtering_time[4];
  int info = 0;
  int* ipiv = (int *) malloc(max_mesh_size * sizeof(int));
  double worst_case_sign_ratio[1] = { 100 };
  double* ab_copy = (double*) malloc(ldab * max_mesh_size * sizeof(double));
  double* ab0_copy = (double*) malloc(ldab * max_mesh_size * sizeof(double));
  double h0;
  double* A0 = (double*) malloc(MAX_TEST_MESH_SIZE * MAX_TEST_MESH_SIZE * sizeof(double));

  /* Find the reference time per iteration for BPF */
  double bpf_time_per_iteration = 0;
  for(int nmc = 0; nmc < N_top; nmc++) {

    for(int i = 0; i < ldab * max_mesh_size; i++)
      ab_copy[i] = ab_bpf[i];

    bootstrapfilter(y, std_sig, obs_std, data_length, x0, bpf_sample_size, NULL, NULL, filtering_time, max_mesh_size, kl, ku, ab_bpf, ldab, ipiv, max_mesh_size, info, w, d, E, I, h_bpf , L, meass, N_MEASUREMENTS);
    bpf_time_per_iteration += filtering_time[0] / (double) N_top;
  }

  printf("Time per iteration (BPF): %0.5f\n",bpf_time_per_iteration);
  fflush(stdout);

  int level_0_mesh_size;
  int sample_sizes[2] = {0,0};
  int mesh_sizes[2];
  int N0;
  short direction;

  FILE* out = fopen("time_matched_sample_sizes.txt","w");

  /* Iterate over level 0 mesh sizes */
  for(int j = 0; j < n_level_0_mesh_sizes; j++) {

    level_0_mesh_size = level_0_mesh_sizes[j];
    mesh_sizes[0] = level_0_mesh_size;
    mesh_sizes[1] = max_mesh_size;

    // We have to create A0 again for each mesh size
    for(int i = 0; i < ldab * mesh_sizes[0]; i++)
      ab0[i] = 0;
    create_dense_solver_matrix(A0, mesh_sizes[0]);
    make_lapack_band_matrix(A0, mesh_sizes[0], ab0, ku, kl);
    h0 = L / (double)(mesh_sizes[0] - 1);

    /* For given L0 mesh, iterate over N1 */
    for(int N1 = 0; N1 <= bpf_sample_size ; N1 = N1 + N1_step) {

      N0 = N0_init;
      N0_step = 10000;

      /* For given L0 mesh and N1 iterate over increasing N0 until MLBPF time per iteration exceeds
       BPF time per iteration */
      time_per_iteration = 0;

      direction = 1; // up
      while(1) {

        sample_sizes[0] = N0;
        sample_sizes[1] = N1;

        /* Estimate the time per iteration by averaging N_top runs */
        time_per_iteration = 0;
        for(int nmc = 0; nmc < N_top; nmc++) {

          for(int i = 0; i < ldab * max_mesh_size; i++)
            ab_copy[i] = ab_bpf[i];

          for(int i = 0; i < ldab * level_0_mesh_size; i++)
            ab0_copy[i] = ab0[i];

          MLbootstrapfilter(y, std_sig, obs_std, data_length, x0, sample_sizes, worst_case_sign_ratio, NULL, filtering_time, mesh_sizes, kl, ku, ab0_copy, ab_copy, ldab, ipiv, ipiv, mesh_sizes[0], max_mesh_size, info, w, d, E, I, h0, h1, L, meass, N_MEASUREMENTS);
          //          printf("TIMES: %e %e %e %e\n", filtering_time[0], filtering_time[1], filtering_time[2], filtering_time[3]);
          //          fflush(stdout);

          time_per_iteration += filtering_time[0] / (double) N_top;
        }

        if (direction == 1 && time_per_iteration > bpf_time_per_iteration * time_tolerance && N0_step > 10) {
          direction = -1; // start going down
          N0_step /= 2; // halve the step

        } else if (direction == -1 && time_per_iteration < bpf_time_per_iteration * 1.00 && N0_step > 10) {
          direction = 1;
          N0_step /= 2;
        } else if (N0_step < 10) {
          break;
        }
        N0 = N0 + direction * N0_step;

        if(N0 < 1) { // Make sure N0 remains positive and go up if this happens
          N0 = 1;
          direction = 1;
        }

      }
      fprintf(out,"%i %i %i %i\n",N0-N0_step,N1,mesh_sizes[0],mesh_sizes[1]);
      fflush(out);
    }
  }
  fclose(out);

  free(ipiv);
  free(ab_copy);
  free(ab0_copy);
  free(A0);
}

double compute_estimation_error(double* x_hat, double* x_ref, int len) {

  double MSE = 0;

  for(int i = 0; i < len; i++) {
    MSE += (x_hat[i]-x_ref[i])*(x_hat[i]-x_ref[i]) / (double) len;
  }
  return MSE;
}

void init_with_zeros(double *array, long length) {
  for (int i = 0; i < length; i++)
    array[i] = 0;
}

void generate_signal(double* x, int len, double std_sig, double x0, gsl_rng* rng, double L) {

  double lower_bound = 0.25;
  double upper_bound = L - 0.25;
  x[0] = x0;

  for (int i = 1; i < len; i++){
    x[i] = x[i - 1] + gsl_ran_gaussian(rng, std_sig);
    //    x[i] = 0.5 / ((double)1.0 + exp(-(double)10.0 * (x[i - 1] - 1))) + (double) 0.75 + gsl_ran_gaussian(rng, std_sig);
    if(x[i] < lower_bound)
      x[i] = lower_bound + lower_bound-x[i];
    if(x[i] > upper_bound)
      x[i] = upper_bound - (x[i]-upper_bound);
  }

}

void make_lapack_band_matrix(double *A, int n, double *ab, int ku, int kl) {

  int lower, upper, baserow,baserowcand ;
  int n_rows = 2 * kl + ku + 1;

  for (int j = 0; j < n; j++) { // iterate columns

    lower = j - ku > 0 ? j - ku : 0;
    upper = j + kl < n ? j + kl : n-1;

    baserowcand = n_rows - upper - 1;
    baserow = baserowcand < n_rows - (ku+kl+1)? n_rows - (ku+kl+1) : baserowcand;

    for (int i = lower; i <= upper; i++) {
      ab[n_rows * j + baserow + i - lower] = A[j * n + i];
    }
  }

}

void create_dense_solver_matrix(double* A, int n) {

  double pattern[5] = { 1, -4, 6, -4, 1 };
  /* Make sure there are zeros */
  for(int i = 0; i < n*n; i++)
    A[i] = 0;

  /* First two lines */
  for (int i = 0; i < n; i++) {
    A[i * n] = 0;
    A[i * n + 1] = 0;
  }

  A[0] = 16;
  A[1] = -4;
  A[n + 0] = -9;
  A[n + 1] = 6;
  A[2 * n + 0] = (double) 8 / (double) 3;
  A[2 * n + 1] = -4;
  A[3 * n + 0] = -(double) 1 / (double) 4;
  A[3 * n + 1] = (double) 1;

  /* Last two lines */
  for (int i = 0; i < n; i++) {
    A[i * n + n - 2] = 0;
    A[i * n + n - 1] = 0;
  }

  //  A[(n - 4) * n + n - 2] = (double) 16 / (double) 17;
  A[(n - 4) * n + n - 2] = (double) 1;
  //  A[(n - 4) * n + n - 1] = -(double) 12 / (double) 17;
  A[(n - 4) * n + n - 1] = -(double) 1 / (double) 4;
  //  A[(n - 3) * n + n - 2] = -(double) 60 / (double) 17;
  A[(n - 3) * n + n - 2] = -(double) 4;
  //  A[(n - 3) * n + n - 1] = (double) 96 / (double) 17;
  A[(n - 3) * n + n - 1] = (double) 8 / (double) 3;
  //  A[(n - 2) * n + n - 2] = (double) 72 / (double) 17;
  A[(n - 2) * n + n - 2] = (double) 6;
  //  A[(n - 2) * n + n - 1] = -(double) 156 / (double) 17;
  A[(n - 2) * n + n - 1] = -(double) 9;
  //  A[(n - 1) * n + n - 2] = -(double) 28 / (double) 17;
  A[(n - 1) * n + n - 2] = -(double) 4;
  //  A[(n - 1) * n + n - 1] = (double) 72 / (double) 17;
  A[(n - 1) * n + n - 1] = (double) 16;

  /* The rest */
  for (int i = 2; i < n - 2; i++) { // row
    for (int j = 0; j < n; j++) {
      A[j * n + i] = 0;
    }
    for (int j = 0; j < 5; j++)
      A[(j + i - 2) * n + i] = pattern[j];
  }
}

void get_settings(int** settings) {

//  /* THIS IS FOR N = 1000 */
//  int hc_settings[N_SETTINGS][4] = {
//  {38000, 0, 100, 4000},
//{36835, 40, 100, 4000},
//{36054, 80, 100, 4000},
//{36439, 120, 100, 4000},
//{35867, 160, 100, 4000},
//{35658, 200, 100, 4000},
//{34954, 240, 100, 4000},
//{34888, 280, 100, 4000},
//{34570, 320, 100, 4000},
//{34237, 360, 100, 4000},
//{33456, 400, 100, 4000},
//{32811, 440, 100, 4000},
//{32127, 480, 100, 4000},
//{30799, 520, 100, 4000},
//{29146, 560, 100, 4000},
//{27180, 600, 100, 4000},
//{24810, 640, 100, 4000},
//{22713, 680, 100, 4000},
//{20024, 720, 100, 4000},
//{17349, 760, 100, 4000},
//{14049, 800, 100, 4000},
//{10800, 840, 100, 4000},
//{7812, 880, 100, 4000},
//{4960, 920, 100, 4000},
//{2408, 960, 100, 4000},
//{36965, 0, 105, 4000},
//{36354, 40, 105, 4000},
//{35552, 80, 105, 4000},
//{34341, 120, 105, 4000},
//{34073, 160, 105, 4000},
//{33840, 200, 105, 4000},
//{33970, 240, 105, 4000},
//{33444, 280, 105, 4000},
//{33079, 320, 105, 4000},
//{32883, 360, 105, 4000},
//{32615, 400, 105, 4000},
//{32173, 440, 105, 4000},
//{31021, 480, 105, 4000},
//{29393, 520, 105, 4000},
//{27453, 560, 105, 4000},
//{25473, 600, 105, 4000},
//{23605, 640, 105, 4000},
//{21639, 680, 105, 4000},
//{19367, 720, 105, 4000},
//{16693, 760, 105, 4000},
//{13229, 800, 105, 4000},
//{10143, 840, 105, 4000},
//{7134, 880, 105, 4000},
//{4543, 920, 105, 4000},
//{2108, 960, 105, 4000},
//{35266, 0, 110, 4000},
//{35461, 40, 110, 4000},
//{35409, 80, 110, 4000},
//{35058, 120, 110, 4000},
//{34531, 160, 110, 4000},
//{33971, 200, 110, 4000},
//{33437, 240, 110, 4000},
//{33026, 280, 110, 4000},
//{32804, 320, 110, 4000},
//{32765, 360, 110, 4000},
//{32550, 400, 110, 4000},
//{31659, 440, 110, 4000},
//{30162, 480, 110, 4000},
//{28254, 520, 110, 4000},
//{26203, 560, 110, 4000},
//{24055, 600, 110, 4000},
//{21965, 640, 110, 4000},
//{19192, 680, 110, 4000},
//{16087, 720, 110, 4000},
//{12969, 760, 110, 4000},
//{10438, 800, 110, 4000},
//{8172, 840, 110, 4000},
//{5045, 880, 110, 4000},
//{2513, 920, 110, 4000},
//{605, 960, 110, 4000},
//{32765, 0, 115, 4000},
//{33038, 40, 115, 4000},
//{33013, 80, 115, 4000},
//{32440, 120, 115, 4000},
//{31593, 160, 115, 4000},
//{31444, 200, 115, 4000},
//{31444, 240, 115, 4000},
//{31235, 280, 115, 4000},
//{30265, 320, 115, 4000},
//{29243, 360, 115, 4000},
//{29035, 400, 115, 4000},
//{28579, 440, 115, 4000},
//{27928, 480, 115, 4000},
//{26274, 520, 115, 4000},
//{24568, 560, 115, 4000},
//{22960, 600, 115, 4000},
//{21743, 640, 115, 4000},
//{19972, 680, 115, 4000},
//{17707, 720, 115, 4000},
//{14634, 760, 115, 4000},
//{11717, 800, 115, 4000},
//{8586, 840, 115, 4000},
//{6040, 880, 115, 4000},
//{3749, 920, 115, 4000},
//{1867, 960, 115, 4000},
//{34171, 0, 120, 4000},
//{33248, 40, 120, 4000},
//{32916, 80, 120, 4000},
//{33066, 120, 120, 4000},
//{32969, 160, 120, 4000},
//{32323, 200, 120, 4000},
//{31901, 240, 120, 4000},
//{31776, 280, 120, 4000},
//{31932, 320, 120, 4000},
//{31580, 360, 120, 4000},
//{29759, 400, 120, 4000},
//{29077, 440, 120, 4000},
//{27950, 480, 120, 4000},
//{27637, 520, 120, 4000},
//{25709, 560, 120, 4000},
//{23834, 600, 120, 4000},
//{22322, 640, 120, 4000},
//{20388, 680, 120, 4000},
//{18175, 720, 120, 4000},
//{15142, 760, 120, 4000},
//{12278, 800, 120, 4000},
//{9471, 840, 120, 4000},
//{6984, 880, 120, 4000},
//{4640, 920, 120, 4000},
//{2362, 960, 120, 4000},
//{32765, 0, 130, 4000},
//{32765, 40, 130, 4000},
//{32622, 80, 130, 4000},
//{32244, 120, 130, 4000},
//{31763, 160, 130, 4000},
//{31249, 200, 130, 4000},
//{30553, 240, 130, 4000},
//{29629, 280, 130, 4000},
//{29062, 320, 130, 4000},
//{28794, 360, 130, 4000},
//{28585, 400, 130, 4000},
//{27778, 440, 130, 4000},
//{26463, 480, 130, 4000},
//{24966, 520, 130, 4000},
//{23411, 560, 130, 4000},
//{22420, 600, 130, 4000},
//{20878, 640, 130, 4000},
//{18827, 680, 130, 4000},
//{15982, 720, 130, 4000},
//{13273, 760, 130, 4000},
//{10487, 800, 130, 4000},
//{7955, 840, 130, 4000},
//{6021, 880, 130, 4000},
//{4524, 920, 130, 4000},
//{2864, 960, 130, 4000},
//{31926, 0, 140, 4000},
//{31405, 40, 140, 4000},
//{30917, 80, 140, 4000},
//{30312, 120, 140, 4000},
//{29876, 160, 140, 4000},
//{29433, 200, 140, 4000},
//{29055, 240, 140, 4000},
//{28579, 280, 140, 4000},
//{28058, 320, 140, 4000},
//{27524, 360, 140, 4000},
//{26952, 400, 140, 4000},
//{26073, 440, 140, 4000},
//{24790, 480, 140, 4000},
//{23332, 520, 140, 4000},
//{22290, 560, 140, 4000},
//{21386, 600, 140, 4000},
//{20128, 640, 140, 4000},
//{17993, 680, 140, 4000},
//{15493, 720, 140, 4000},
//{12941, 760, 140, 4000},
//{10454, 800, 140, 4000},
//{7915, 840, 140, 4000},
//{4605, 880, 140, 4000},
//{1955, 920, 140, 4000},
//{151, 960, 140, 4000},
//{30538, 0, 150, 4000},
//{29758, 40, 150, 4000},
//{29140, 80, 150, 4000},
//{28560, 120, 150, 4000},
//{27948, 160, 150, 4000},
//{27336, 200, 150, 4000},
//{26984, 240, 150, 4000},
//{26607, 280, 150, 4000},
//{26458, 320, 150, 4000},
//{26295, 360, 150, 4000},
//{26256, 400, 150, 4000},
//{25643, 440, 150, 4000},
//{24537, 480, 150, 4000},
//{23157, 520, 150, 4000},
//{22245, 560, 150, 4000},
//{21144, 600, 150, 4000},
//{19712, 640, 150, 4000},
//{17473, 680, 150, 4000},
//{15058, 720, 150, 4000},
//{12538, 760, 150, 4000},
//{10024, 800, 150, 4000},
//{7544, 840, 150, 4000},
//{5396, 880, 150, 4000},
//{3281, 920, 150, 4000},
//{1581, 960, 150, 4000},
//{22765, 0, 160, 4000},
//{22186, 40, 160, 4000},
//{21873, 80, 160, 4000},
//{21821, 120, 160, 4000},
//{21827, 160, 160, 4000},
//{21827, 200, 160, 4000},
//{21827, 240, 160, 4000},
//{21801, 280, 160, 4000},
//{21217, 320, 160, 4000},
//{20965, 360, 160, 4000},
//{20026, 400, 160, 4000},
//{19770, 440, 160, 4000},
//{18369, 480, 160, 4000},
//{17388, 520, 160, 4000},
//{15949, 560, 160, 4000},
//{15012, 600, 160, 4000},
//{13845, 640, 160, 4000},
//{12595, 680, 160, 4000},
//{11241, 720, 160, 4000},
//{9887, 760, 160, 4000},
//{8306, 800, 160, 4000},
//{6561, 840, 160, 4000},
//{4634, 880, 160, 4000},
//{2800, 920, 160, 4000},
//{1212, 960, 160, 4000},
//{26379, 0, 180, 4000},
//{26815, 40, 180, 4000},
//{26880, 80, 180, 4000},
//{26556, 120, 180, 4000},
//{26235, 160, 180, 4000},
//{25864, 200, 180, 4000},
//{25637, 240, 180, 4000},
//{25246, 280, 180, 4000},
//{24940, 320, 180, 4000},
//{24399, 360, 180, 4000},
//{23846, 400, 180, 4000},
//{23163, 440, 180, 4000},
//{22564, 480, 180, 4000},
//{21874, 520, 180, 4000},
//{20963, 560, 180, 4000},
//{19739, 600, 180, 4000},
//{17562, 640, 180, 4000},
//{14182, 680, 180, 4000},
//{11115, 720, 180, 4000},
//{8364, 760, 180, 4000},
//{6770, 800, 180, 4000},
//{5407, 840, 180, 4000},
//{4314, 880, 180, 4000},
//{3304, 920, 180, 4000},
//{1678, 960, 180, 4000},
//{25989, 0, 200, 4000},
//{25488, 40, 200, 4000},
//{25142, 80, 200, 4000},
//{24654, 120, 200, 4000},
//{24510, 160, 200, 4000},
//{24119, 200, 200, 4000},
//{23787, 240, 200, 4000},
//{23384, 280, 200, 4000},
//{22974, 320, 200, 4000},
//{22558, 360, 200, 4000},
//{22024, 400, 200, 4000},
//{21729, 440, 200, 4000},
//{21020, 480, 200, 4000},
//{19978, 520, 200, 4000},
//{18430, 560, 200, 4000},
//{16893, 600, 200, 4000},
//{15331, 640, 200, 4000},
//{13443, 680, 200, 4000},
//{11502, 720, 200, 4000},
//{9549, 760, 200, 4000},
//{7778, 800, 200, 4000},
//{6007, 840, 200, 4000},
//{4360, 880, 200, 4000},
//{2870, 920, 200, 4000},
//{1486, 960, 200, 4000},
//{23624, 0, 220, 4000},
//{23285, 40, 220, 4000},
//{22934, 80, 220, 4000},
//{22472, 120, 220, 4000},
//{22075, 160, 220, 4000},
//{21827, 200, 220, 4000},
//{21814, 240, 220, 4000},
//{21678, 280, 220, 4000},
//{21327, 320, 220, 4000},
//{20871, 360, 220, 4000},
//{20343, 400, 220, 4000},
//{19595, 440, 220, 4000},
//{18599, 480, 220, 4000},
//{17219, 520, 220, 4000},
//{15890, 560, 220, 4000},
//{14406, 600, 220, 4000},
//{12882, 640, 220, 4000},
//{11105, 680, 220, 4000},
//{9314, 720, 220, 4000},
//{7517, 760, 220, 4000},
//{5818, 800, 220, 4000},
//{4165, 840, 220, 4000},
//{2694, 880, 220, 4000},
//{1382, 920, 220, 4000},
//{497, 960, 220, 4000},
//{21398, 0, 240, 4000},
//{21047, 40, 240, 4000},
//{20754, 80, 240, 4000},
//{20422, 120, 240, 4000},
//{20154, 160, 240, 4000},
//{19861, 200, 240, 4000},
//{19607, 240, 240, 4000},
//{19269, 280, 240, 4000},
//{18892, 320, 240, 4000},
//{18508, 360, 240, 4000},
//{18098, 400, 240, 4000},
//{17518, 440, 240, 4000},
//{16593, 480, 240, 4000},
//{15441, 520, 240, 4000},
//{14230, 560, 240, 4000},
//{12994, 600, 240, 4000},
//{11594, 640, 240, 4000},
//{9999, 680, 240, 4000},
//{8320, 720, 240, 4000},
//{6698, 760, 240, 4000},
//{5188, 800, 240, 4000},
//{3827, 840, 240, 4000},
//{2532, 880, 240, 4000},
//{1271, 920, 240, 4000},
//{432, 960, 240, 4000},
//{21730, 0, 260, 4000},
//{21210, 40, 260, 4000},
//{20839, 80, 260, 4000},
//{20519, 120, 260, 4000},
//{20083, 160, 260, 4000},
//{19627, 200, 260, 4000},
//{19223, 240, 260, 4000},
//{18819, 280, 260, 4000},
//{18383, 320, 260, 4000},
//{18436, 360, 260, 4000},
//{18449, 400, 260, 4000},
//{18221, 440, 260, 4000},
//{17270, 480, 260, 4000},
//{16079, 520, 260, 4000},
//{14790, 560, 260, 4000},
//{13515, 600, 260, 4000},
//{12141, 640, 260, 4000},
//{10780, 680, 260, 4000},
//{9257, 720, 260, 4000},
//{7811, 760, 260, 4000},
//{6313, 800, 260, 4000},
//{4894, 840, 260, 4000},
//{3631, 880, 260, 4000},
//{2545, 920, 260, 4000},
//{1614, 960, 260, 4000},
//{21262, 0, 280, 4000},
//{20890, 40, 280, 4000},
//{20545, 80, 280, 4000},
//{20129, 120, 280, 4000},
//{19771, 160, 280, 4000},
//{19367, 200, 280, 4000},
//{18785, 240, 280, 4000},
//{18355, 280, 280, 4000},
//{17965, 320, 280, 4000},
//{17739, 360, 280, 4000},
//{17310, 400, 280, 4000},
//{16619, 440, 280, 4000},
//{15714, 480, 280, 4000},
//{14588, 520, 280, 4000},
//{13403, 560, 280, 4000},
//{12238, 600, 280, 4000},
//{10988, 640, 280, 4000},
//{9673, 680, 280, 4000},
//{8293, 720, 280, 4000},
//{6984, 760, 280, 4000},
//{5690, 800, 280, 4000},
//{4491, 840, 280, 4000},
//{3359, 880, 280, 4000},
//{2362, 920, 280, 4000},
//{1388, 960, 280, 4000},
//{13585, 0, 320, 4000},
//{13391, 40, 320, 4000},
//{13202, 80, 320, 4000},
//{12936, 120, 320, 4000},
//{12610, 160, 320, 4000},
//{12252, 200, 320, 4000},
//{11985, 240, 320, 4000},
//{11776, 280, 320, 4000},
//{11567, 320, 320, 4000},
//{11320, 360, 320, 4000},
//{10989, 400, 320, 4000},
//{10585, 440, 320, 4000},
//{9993, 480, 320, 4000},
//{9263, 520, 320, 4000},
//{8462, 560, 320, 4000},
//{7785, 600, 320, 4000},
//{7232, 640, 320, 4000},
//{6685, 680, 320, 4000},
//{6047, 720, 320, 4000},
//{5304, 760, 320, 4000},
//{4503, 800, 320, 4000},
//{3598, 840, 320, 4000},
//{2648, 880, 320, 4000},
//{1737, 920, 320, 4000},
//{881, 960, 320, 4000},
//{16086, 0, 360, 4000},
//{15669, 40, 360, 4000},
//{15330, 80, 360, 4000},
//{15096, 120, 360, 4000},
//{14796, 160, 360, 4000},
//{14465, 200, 360, 4000},
//{14075, 240, 360, 4000},
//{13743, 280, 360, 4000},
//{13455, 320, 360, 4000},
//{13150, 360, 360, 4000},
//{12786, 400, 360, 4000},
//{12284, 440, 360, 4000},
//{11554, 480, 360, 4000},
//{10630, 520, 360, 4000},
//{9686, 560, 360, 4000},
//{8827, 600, 360, 4000},
//{7960, 640, 360, 4000},
//{6945, 680, 360, 4000},
//{5904, 720, 360, 4000},
//{4824, 760, 360, 4000},
//{3840, 800, 360, 4000},
//{2922, 840, 360, 4000},
//{2108, 880, 360, 4000},
//{1386, 920, 360, 4000},
//{682, 960, 360, 4000},
//{13802, 0, 400, 4000},
//{13579, 40, 400, 4000},
//{13338, 80, 400, 4000},
//{13098, 120, 400, 4000},
//{12799, 160, 400, 4000},
//{12500, 200, 400, 4000},
//{12180, 240, 400, 4000},
//{11900, 280, 400, 4000},
//{11613, 320, 400, 4000},
//{11327, 360, 400, 4000},
//{11021, 400, 400, 4000},
//{10558, 440, 400, 4000},
//{9920, 480, 400, 4000},
//{9061, 520, 400, 4000},
//{8260, 560, 400, 4000},
//{7551, 600, 400, 4000},
//{6847, 640, 400, 4000},
//{5988, 680, 400, 4000},
//{4992, 720, 400, 4000},
//{4048, 760, 400, 4000},
//{3188, 800, 400, 4000},
//{2406, 840, 400, 4000},
//{1684, 880, 400, 4000},
//{1027, 920, 400, 4000},
//{480, 960, 400, 4000},
//{12200, 0, 440, 4000},
//{11984, 40, 440, 4000},
//{11737, 80, 440, 4000},
//{11483, 120, 440, 4000},
//{11229, 160, 440, 4000},
//{10949, 200, 440, 4000},
//{10631, 240, 440, 4000},
//{10365, 280, 440, 4000},
//{10156, 320, 440, 4000},
//{9946, 360, 440, 4000},
//{9646, 400, 440, 4000},
//{9224, 440, 440, 4000},
//{8657, 480, 440, 4000},
//{7941, 520, 440, 4000},
//{7270, 560, 440, 4000},
//{6593, 600, 440, 4000},
//{5878, 640, 440, 4000},
//{5039, 680, 440, 4000},
//{4200, 720, 440, 4000},
//{3424, 760, 440, 4000},
//{2700, 800, 440, 4000},
//{2017, 840, 440, 4000},
//{1385, 880, 440, 4000},
//{838, 920, 440, 4000},
//{382, 960, 440, 4000},
//{10206, 0, 500, 4000},
//{10018, 40, 500, 4000},
//{9829, 80, 500, 4000},
//{9608, 120, 500, 4000},
//{9386, 160, 500, 4000},
//{9152, 200, 500, 4000},
//{8931, 240, 500, 4000},
//{8697, 280, 500, 4000},
//{8476, 320, 500, 4000},
//{8261, 360, 500, 4000},
//{8027, 400, 500, 4000},
//{7669, 440, 500, 4000},
//{7200, 480, 500, 4000},
//{6608, 520, 500, 4000},
//{5989, 560, 500, 4000},
//{5357, 600, 500, 4000},
//{4699, 640, 500, 4000},
//{4042, 680, 500, 4000},
//{3398, 720, 500, 4000},
//{2778, 760, 500, 4000},
//{2160, 800, 500, 4000},
//{1554, 840, 500, 4000},
//{1036, 880, 500, 4000},
//{538, 920, 500, 4000},
//    {200, 960, 500, 4000}};

  /* THIS IS FOR N = 500 */
  int hc_settings[N_SETTINGS][4] = {
  {19432, 0, 100, 4000},
{17976, 20, 100, 4000},
{17413, 40, 100, 4000},
{17524, 60, 100, 4000},
{17159, 80, 100, 4000},
{16432, 100, 100, 4000},
{15611, 120, 100, 4000},
{14817, 140, 100, 4000},
{13937, 160, 100, 4000},
{13338, 180, 100, 4000},
{12623, 200, 100, 4000},
{12284, 220, 100, 4000},
{11594, 240, 100, 4000},
{11080, 260, 100, 4000},
{10267, 280, 100, 4000},
{9648, 300, 100, 4000},
{8933, 320, 100, 4000},
{8340, 340, 100, 4000},
{7630, 360, 100, 4000},
{6983, 380, 100, 4000},
{6346, 400, 100, 4000},
{5709, 420, 100, 4000},
{4680, 440, 100, 4000},
{3456, 460, 100, 4000},
{1792, 480, 100, 4000},
{18429, 0, 105, 4000},
{17238, 20, 105, 4000},
{16495, 40, 105, 4000},
{15845, 60, 105, 4000},
{15584, 80, 105, 4000},
{15136, 100, 105, 4000},
{14420, 120, 105, 4000},
{13951, 140, 105, 4000},
{13261, 160, 105, 4000},
{12909, 180, 105, 4000},
{12408, 200, 105, 4000},
{11978, 220, 105, 4000},
{11398, 240, 105, 4000},
{10838, 260, 105, 4000},
{9959, 280, 105, 4000},
{9374, 300, 105, 4000},
{8715, 320, 105, 4000},
{8337, 340, 105, 4000},
{7706, 360, 105, 4000},
{7121, 380, 105, 4000},
{6425, 400, 105, 4000},
{5650, 420, 105, 4000},
{4310, 440, 105, 4000},
{3021, 460, 105, 4000},
{1423, 480, 105, 4000},
{17632, 0, 110, 4000},
{16920, 20, 110, 4000},
{16340, 40, 110, 4000},
{15910, 60, 110, 4000},
{15402, 80, 110, 4000},
{14849, 100, 110, 4000},
{14244, 120, 110, 4000},
{13651, 140, 110, 4000},
{13058, 160, 110, 4000},
{12486, 180, 110, 4000},
{11900, 200, 110, 4000},
{11366, 220, 110, 4000},
{10819, 240, 110, 4000},
{10260, 260, 110, 4000},
{9596, 280, 110, 4000},
{9017, 300, 110, 4000},
{8371, 320, 110, 4000},
{7967, 340, 110, 4000},
{7381, 360, 110, 4000},
{6867, 380, 110, 4000},
{6177, 400, 110, 4000},
{5520, 420, 110, 4000},
{4648, 440, 110, 4000},
{3398, 460, 110, 4000},
{1789, 480, 110, 4000},
{16904, 0, 115, 4000},
{16411, 20, 115, 4000},
{15943, 40, 115, 4000},
{15500, 60, 115, 4000},
{15057, 80, 115, 4000},
{14601, 100, 115, 4000},
{14010, 120, 115, 4000},
{13366, 140, 115, 4000},
{12650, 160, 115, 4000},
{12082, 180, 115, 4000},
{11554, 200, 115, 4000},
{11138, 220, 115, 4000},
{10695, 240, 115, 4000},
{10200, 260, 115, 4000},
{9588, 280, 115, 4000},
{8970, 300, 115, 4000},
{8449, 320, 115, 4000},
{7889, 340, 115, 4000},
{7336, 360, 115, 4000},
{6686, 380, 115, 4000},
{6133, 400, 115, 4000},
{5494, 420, 115, 4000},
{4497, 440, 115, 4000},
{3266, 460, 115, 4000},
{1665, 480, 115, 4000},
{16498, 0, 120, 4000},
{16055, 20, 120, 4000},
{15631, 40, 120, 4000},
{15194, 60, 120, 4000},
{14751, 80, 120, 4000},
{14295, 100, 120, 4000},
{13624, 120, 120, 4000},
{12993, 140, 120, 4000},
{12343, 160, 120, 4000},
{11848, 180, 120, 4000},
{11333, 200, 120, 4000},
{10760, 220, 120, 4000},
{10213, 240, 120, 4000},
{9692, 260, 120, 4000},
{9204, 280, 120, 4000},
{8723, 300, 120, 4000},
{8189, 320, 120, 4000},
{7648, 340, 120, 4000},
{7049, 360, 120, 4000},
{6516, 380, 120, 4000},
{5917, 400, 120, 4000},
{5214, 420, 120, 4000},
{4328, 440, 120, 4000},
{3116, 460, 120, 4000},
{1658, 480, 120, 4000},
{14992, 0, 130, 4000},
{14608, 20, 130, 4000},
{14224, 40, 130, 4000},
{13891, 60, 130, 4000},
{13488, 80, 130, 4000},
{13065, 100, 130, 4000},
{12538, 120, 130, 4000},
{11940, 140, 130, 4000},
{11386, 160, 130, 4000},
{10872, 180, 130, 4000},
{10390, 200, 130, 4000},
{9856, 220, 130, 4000},
{9420, 240, 130, 4000},
{8873, 260, 130, 4000},
{8358, 280, 130, 4000},
{7739, 300, 130, 4000},
{7264, 320, 130, 4000},
{6789, 340, 130, 4000},
{6366, 360, 130, 4000},
{5903, 380, 130, 4000},
{5448, 400, 130, 4000},
{4797, 420, 130, 4000},
{3938, 440, 130, 4000},
{2830, 460, 130, 4000},
{1496, 480, 130, 4000},
{14152, 0, 140, 4000},
{13710, 20, 140, 4000},
{13286, 40, 140, 4000},
{12941, 60, 140, 4000},
{12505, 80, 140, 4000},
{12075, 100, 140, 4000},
{11554, 120, 140, 4000},
{11098, 140, 140, 4000},
{10604, 160, 140, 4000},
{10200, 180, 140, 4000},
{9692, 200, 140, 4000},
{9210, 220, 140, 4000},
{8748, 240, 140, 4000},
{8305, 260, 140, 4000},
{7850, 280, 140, 4000},
{7264, 300, 140, 4000},
{6750, 320, 140, 4000},
{6307, 340, 140, 4000},
{5884, 360, 140, 4000},
{5447, 380, 140, 4000},
{4952, 400, 140, 4000},
{4458, 420, 140, 4000},
{3729, 440, 140, 4000},
{2668, 460, 140, 4000},
{1366, 480, 140, 4000},
{13138, 0, 150, 4000},
{12792, 20, 150, 4000},
{12551, 40, 150, 4000},
{12212, 60, 150, 4000},
{11828, 80, 150, 4000},
{11469, 100, 150, 4000},
{11073, 120, 150, 4000},
{10643, 140, 150, 4000},
{10025, 160, 150, 4000},
{9536, 180, 150, 4000},
{9159, 200, 150, 4000},
{8767, 220, 150, 4000},
{8338, 240, 150, 4000},
{7843, 260, 150, 4000},
{7453, 280, 150, 4000},
{7062, 300, 150, 4000},
{6593, 320, 150, 4000},
{6072, 340, 150, 4000},
{5682, 360, 150, 4000},
{5226, 380, 150, 4000},
{4849, 400, 150, 4000},
{4354, 420, 150, 4000},
{3742, 440, 150, 4000},
{2759, 460, 150, 4000},
{1447, 480, 150, 4000},
{11320, 0, 160, 4000},
{11060, 20, 160, 4000},
{10787, 40, 160, 4000},
{10513, 60, 160, 4000},
{10083, 80, 160, 4000},
{9712, 100, 160, 4000},
{9485, 120, 160, 4000},
{9276, 140, 160, 4000},
{9037, 160, 160, 4000},
{8743, 180, 160, 4000},
{8437, 200, 160, 4000},
{8110, 220, 160, 4000},
{7680, 240, 160, 4000},
{7257, 260, 160, 4000},
{6867, 280, 160, 4000},
{6483, 300, 160, 4000},
{6053, 320, 160, 4000},
{5643, 340, 160, 4000},
{5220, 360, 160, 4000},
{4803, 380, 160, 4000},
{4379, 400, 160, 4000},
{3892, 420, 160, 4000},
{3241, 440, 160, 4000},
{2303, 460, 160, 4000},
{1189, 480, 160, 4000},
{11125, 0, 180, 4000},
{10799, 20, 180, 4000},
{10526, 40, 180, 4000},
{10213, 60, 180, 4000},
{9920, 80, 180, 4000},
{9484, 100, 180, 4000},
{9016, 120, 180, 4000},
{8553, 140, 180, 4000},
{8156, 160, 180, 4000},
{7843, 180, 180, 4000},
{7576, 200, 180, 4000},
{7297, 220, 180, 4000},
{6985, 240, 180, 4000},
{6659, 260, 180, 4000},
{6346, 280, 180, 4000},
{5994, 300, 180, 4000},
{5604, 320, 180, 4000},
{5161, 340, 180, 4000},
{4751, 360, 180, 4000},
{4335, 380, 180, 4000},
{3963, 400, 180, 4000},
{3410, 420, 180, 4000},
{2720, 440, 180, 4000},
{1802, 460, 180, 4000},
{891, 480, 180, 4000},
{9699, 0, 200, 4000},
{9426, 20, 200, 4000},
{9159, 40, 200, 4000},
{8892, 60, 200, 4000},
{8612, 80, 200, 4000},
{8293, 100, 200, 4000},
{7955, 120, 200, 4000},
{7590, 140, 200, 4000},
{7284, 160, 200, 4000},
{6938, 180, 200, 4000},
{6600, 200, 200, 4000},
{6229, 220, 200, 4000},
{5891, 240, 200, 4000},
{5578, 260, 200, 4000},
{5252, 280, 200, 4000},
{4947, 300, 200, 4000},
{4576, 320, 200, 4000},
{4244, 340, 200, 4000},
{3879, 360, 200, 4000},
{3599, 380, 200, 4000},
{3286, 400, 200, 4000},
{2928, 420, 200, 4000},
{2368, 440, 200, 4000},
{1633, 460, 200, 4000},
{813, 480, 200, 4000},
{8859, 0, 220, 4000},
{8637, 20, 220, 4000},
{8397, 40, 220, 4000},
{8130, 60, 220, 4000},
{7844, 80, 220, 4000},
{7563, 100, 220, 4000},
{7277, 120, 220, 4000},
{6932, 140, 220, 4000},
{6601, 160, 220, 4000},
{6255, 180, 220, 4000},
{5969, 200, 220, 4000},
{5682, 220, 220, 4000},
{5383, 240, 220, 4000},
{5057, 260, 220, 4000},
{4732, 280, 220, 4000},
{4439, 300, 220, 4000},
{4166, 320, 220, 4000},
{3840, 340, 220, 4000},
{3534, 360, 220, 4000},
{3214, 380, 220, 4000},
{2889, 400, 220, 4000},
{2512, 420, 220, 4000},
{2030, 440, 220, 4000},
{1444, 460, 220, 4000},
{741, 480, 220, 4000},
{7980, 0, 240, 4000},
{7766, 20, 240, 4000},
{7538, 40, 240, 4000},
{7323, 60, 240, 4000},
{7115, 80, 240, 4000},
{6848, 100, 240, 4000},
{6561, 120, 240, 4000},
{6196, 140, 240, 4000},
{5897, 160, 240, 4000},
{5604, 180, 240, 4000},
{5331, 200, 240, 4000},
{5064, 220, 240, 4000},
{4784, 240, 240, 4000},
{4491, 260, 240, 4000},
{4237, 280, 240, 4000},
{3969, 300, 240, 4000},
{3729, 320, 240, 4000},
{3436, 340, 240, 4000},
{3143, 360, 240, 4000},
{2850, 380, 240, 4000},
{2563, 400, 240, 4000},
{2277, 420, 240, 4000},
{1854, 440, 240, 4000},
{1294, 460, 240, 4000},
{636, 480, 240, 4000},
{7414, 0, 260, 4000},
{7180, 20, 260, 4000},
{6939, 40, 260, 4000},
{6712, 60, 260, 4000},
{6471, 80, 260, 4000},
{6262, 100, 260, 4000},
{6008, 120, 260, 4000},
{5734, 140, 260, 4000},
{5402, 160, 260, 4000},
{5142, 180, 260, 4000},
{4920, 200, 260, 4000},
{4712, 220, 260, 4000},
{4458, 240, 260, 4000},
{4178, 260, 260, 4000},
{3943, 280, 260, 4000},
{3696, 300, 260, 4000},
{3449, 320, 260, 4000},
{3189, 340, 260, 4000},
{2949, 360, 260, 4000},
{2701, 380, 260, 4000},
{2434, 400, 260, 4000},
{2121, 420, 260, 4000},
{1698, 440, 260, 4000},
{1157, 460, 260, 4000},
{565, 480, 260, 4000},
{6652, 0, 280, 4000},
{6503, 20, 280, 4000},
{6301, 40, 280, 4000},
{6093, 60, 280, 4000},
{5884, 80, 280, 4000},
{5689, 100, 280, 4000},
{5441, 120, 280, 4000},
{5193, 140, 280, 4000},
{4934, 160, 280, 4000},
{4693, 180, 280, 4000},
{4452, 200, 280, 4000},
{4217, 220, 280, 4000},
{3996, 240, 280, 4000},
{3775, 260, 280, 4000},
{3567, 280, 280, 4000},
{3326, 300, 280, 4000},
{3098, 320, 280, 4000},
{2850, 340, 280, 4000},
{2616, 360, 280, 4000},
{2382, 380, 280, 4000},
{2134, 400, 280, 4000},
{1880, 420, 280, 4000},
{1502, 440, 280, 4000},
{1033, 460, 280, 4000},
{499, 480, 280, 4000},
{5226, 0, 320, 4000},
{5142, 20, 320, 4000},
{5038, 40, 320, 4000},
{4933, 60, 320, 4000},
{4797, 80, 320, 4000},
{4654, 100, 320, 4000},
{4498, 120, 320, 4000},
{4335, 140, 320, 4000},
{4172, 160, 320, 4000},
{4002, 180, 320, 4000},
{3821, 200, 320, 4000},
{3580, 220, 320, 4000},
{3365, 240, 320, 4000},
{3157, 260, 320, 4000},
{2987, 280, 320, 4000},
{2785, 300, 320, 4000},
{2596, 320, 320, 4000},
{2414, 340, 320, 4000},
{2212, 360, 320, 4000},
{1978, 380, 320, 4000},
{1757, 400, 320, 4000},
{1528, 420, 320, 4000},
{1242, 440, 320, 4000},
{812, 460, 320, 4000},
{376, 480, 320, 4000},
{4874, 0, 360, 4000},
{4718, 20, 360, 4000},
{4556, 40, 360, 4000},
{4380, 60, 360, 4000},
{4211, 80, 360, 4000},
{4062, 100, 360, 4000},
{3906, 120, 360, 4000},
{3723, 140, 360, 4000},
{3527, 160, 360, 4000},
{3358, 180, 360, 4000},
{3202, 200, 360, 4000},
{3032, 220, 360, 4000},
{2863, 240, 360, 4000},
{2720, 260, 360, 4000},
{2571, 280, 360, 4000},
{2388, 300, 360, 4000},
{2193, 320, 360, 4000},
{2004, 340, 360, 4000},
{1847, 360, 360, 4000},
{1678, 380, 360, 4000},
{1509, 400, 360, 4000},
{1294, 420, 360, 4000},
{1047, 440, 360, 4000},
{688, 460, 360, 4000},
{337, 480, 360, 4000},
{4153, 0, 400, 4000},
{4042, 20, 400, 4000},
{3911, 40, 400, 4000},
{3787, 60, 400, 4000},
{3657, 80, 400, 4000},
{3533, 100, 400, 4000},
{3364, 120, 400, 4000},
{3189, 140, 400, 4000},
{3006, 160, 400, 4000},
{2889, 180, 400, 4000},
{2752, 200, 400, 4000},
{2635, 220, 400, 4000},
{2465, 240, 400, 4000},
{2316, 260, 400, 4000},
{2166, 280, 400, 4000},
{2024, 300, 400, 4000},
{1874, 320, 400, 4000},
{1711, 340, 400, 4000},
{1561, 360, 400, 4000},
{1418, 380, 400, 4000},
{1262, 400, 400, 4000},
{1093, 420, 400, 4000},
{892, 440, 400, 4000},
{598, 460, 400, 4000},
{292, 480, 400, 4000},
{3664, 0, 440, 4000},
{3541, 20, 440, 4000},
{3437, 40, 440, 4000},
{3332, 60, 440, 4000},
{3221, 80, 440, 4000},
{3091, 100, 440, 4000},
{2948, 120, 440, 4000},
{2805, 140, 440, 4000},
{2674, 160, 440, 4000},
{2544, 180, 440, 4000},
{2420, 200, 440, 4000},
{2303, 220, 440, 4000},
{2192, 240, 440, 4000},
{2062, 260, 440, 4000},
{1919, 280, 440, 4000},
{1782, 300, 440, 4000},
{1652, 320, 440, 4000},
{1522, 340, 440, 4000},
{1398, 360, 440, 4000},
{1255, 380, 440, 4000},
{1105, 400, 440, 4000},
{923, 420, 440, 4000},
{766, 440, 440, 4000},
{512, 460, 440, 4000},
{265, 480, 440, 4000},
{3077, 0, 500, 4000},
{2974, 20, 500, 4000},
{2948, 40, 500, 4000},
{2877, 60, 500, 4000},
{2799, 80, 500, 4000},
{2629, 100, 500, 4000},
{2498, 120, 500, 4000},
{2446, 140, 500, 4000},
{2388, 160, 500, 4000},
{2349, 180, 500, 4000},
{2161, 200, 500, 4000},
{1991, 220, 500, 4000},
{1815, 240, 500, 4000},
{1704, 260, 500, 4000},
{1593, 280, 500, 4000},
{1470, 300, 500, 4000},
{1354, 320, 500, 4000},
{1236, 340, 500, 4000},
{1125, 360, 500, 4000},
{1001, 380, 500, 4000},
{871, 400, 500, 4000},
{741, 420, 500, 4000},
{611, 440, 500, 4000},
{379, 460, 500, 4000},
    {171, 480, 500, 4000}};

  for(int i = 0; i < N_SETTINGS; i++)
    for(int j = 0; j < 4; j++)
      settings[i][j] = hc_settings[i][j];

}

max_stars_repo_name: heinekmp/MLBPF

max_issues_repo_name: heinekmp/MLBPF

max_forks_repo_name: heinekmp/MLBPF


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef batch_8317ae94_87f4_4cc2_8dab_1eef2919a461_h
#define batch_8317ae94_87f4_4cc2_8dab_1eef2919a461_h

#include <ariel/config.h>
#include <gslib/rtree.h>
#include <ariel/loopblinn.h>

__ariel_begin__

struct bat_triangle;
struct bat_line;
struct bat_batch;

typedef vector<bat_triangle*> bat_triangles;
typedef vector<bat_line*> bat_lines;
typedef rtree_entity<bat_triangle*> bat_rtree_entity;
typedef rtree_node<bat_rtree_entity> bat_rtree_node;
typedef _tree_allocator<bat_rtree_node> bat_rtree_alloc;
typedef tree<bat_rtree_entity, bat_rtree_node, bat_rtree_alloc> bat_tree;
typedef rtree<bat_rtree_entity, quadratic_split_alg<16, 6, bat_tree>, bat_rtree_node, bat_rtree_alloc> bat_rtree;
typedef vector<bat_batch*> bat_batches;

enum bat_type
{
    bf_start,
    bf_cr = bf_start,       /* pt, cr */
    bf_klm_cr,              /* pt, klm, cr */
    bf_klm_tex,             /* pt, klm, tex */
    bf_end = bf_klm_tex,

    bs_start,
    bs_coef_cr = bs_start,  /* pt, coef, cr */
    bs_coef_tex,            /* pt, coef, tex */
    bs_end = bs_coef_tex,
};

struct bat_triangle
{
    lb_joint*           _joints[3];
    vec2                _reduced[3];        /* get a reduced triangle so that no overlapping caused by error would be detected. */
    float               _zorder;
    bool                _is_reduced;

public:
    bat_triangle();
    bat_triangle(lb_joint* i, lb_joint* j, lb_joint* k);
    ~bat_triangle() {}
    void make_rect(rectf& rc) const;
    void set_joint(int i, lb_joint* p) { _joints[i] = p; }
    lb_joint* get_joint(int i) const { return _joints[i]; }
    const vec2& get_point(int i) const { return _joints[i]->get_point(); }
    void* get_lb_binding(int i) const { return _joints[i]->get_binding(); }
    bat_type decide(uint brush_tag) const;
    bool has_klm_coords() const;
    const vec2& get_reduced_point(int i) const;
    bool is_overlapped(const bat_triangle& other) const;
    void ensure_make_reduced();
    void set_zorder(float z) { _zorder = z; }
    float get_zorder() const { return _zorder; }
    void get_center(vec2& c) const;
    void tracing() const;
    void trace_reduced_points() const;
};

struct bat_line
{
    vec2                _points[2];
    vec2                _contourpt[4];      /* calculated contour pts */
    lb_joint*           _srcjs[2];
    vec3                _coef;
    float               _width;
    float               _zorder;
    uint                _tag;               /* pen tag */
    bool                _half;              /* is half line? */
    bool                _recalc;            /* need recalculation? */

public:
    bat_line();
    void set_start_point(const vec2& p) { _points[0] = p; }
    void set_end_point(const vec2& p) { _points[1] = p; }
    void set_source_joint(int i, lb_joint* j) { _srcjs[i] = j; }
    lb_joint* get_source_joint(int i) const { return _srcjs[i]; }
    const vec2& get_start_point() const { return _points[0]; }
    const vec2& get_end_point() const { return _points[1]; }
    void set_pen_tag(uint t) { _tag = t; }
    void setup_coef();
    void set_coef(const vec3& c) { _coef = c; }
    const vec3& get_coef() const { return _coef; }
    void set_zorder(float z) { _zorder = z; }
    float get_zorder() const { return _zorder; }
    void set_line_width(float w) { _width = w; }
    float get_line_width() const { return _width; }
    void set_half_line(bool hl) { _half = hl; }
    bool is_half_line() const { return _half; }
    bat_type decide() const;
    void get_bound_rect(rectf& rc) const;
    int clip_triangle(bat_line output[2], const bat_triangle* triangle) const;
    void calc_contour_points();
    const vec2& get_contour_point(int i) const { return _contourpt[i]; }
    void set_contour_point(int i, const vec2& p) { _contourpt[i] = p; }
    bool set_need_recalc(bool b) { _recalc = b; }
    bool need_recalc() const { return _recalc; }
    void trim_contour(bat_line& line);
    void tracing() const;

public:
    static bat_line* create_line(lb_joint* i, lb_joint* j, float w, float z, uint t, bool half);
    static bat_line* create_half_line(lb_joint* i, lb_joint* j, const vec2& p1, const vec2& p2, float w, float z, uint t);
};

struct bat_batch
{
    bat_type            _type;

public:
    bat_batch(bat_type t): _type(t) {}
    virtual ~bat_batch() {}
    bat_type get_type() const { return _type; }
};

struct bat_fill_batch:
    public bat_batch
{
    bat_rtree           _rtree;

public:
    bat_fill_batch(bat_type t): bat_batch(t) {}
    bat_rtree& get_rtree() { return _rtree; }
    const bat_rtree& const_rtree() const { return _rtree; }
};

struct bat_stroke_batch:
    public bat_batch
{
    bat_lines           _lines;

public:
    bat_stroke_batch(bat_type t): bat_batch(t) {}
    void add_line(bat_line* p) { _lines.push_back(p); }
    bat_lines& get_lines() { return _lines; }
    const bat_lines& const_lines() const { return _lines; }
};

struct bat_stroke_host_batch:
    public bat_stroke_batch
{
public:
    bat_stroke_host_batch(bat_type t): bat_stroke_batch(t) {}
    virtual ~bat_stroke_host_batch();
};

class batch_processor
{
public:
    typedef bat_batches::iterator bat_iter;
    typedef bat_batches::const_iterator bat_const_iter;
    typedef bat_batches::reverse_iterator bat_reversed_iter;
    friend class rose;

public:
    batch_processor();
    ~batch_processor();
    void set_antialias(bool b) { _antialias = b; }
    bool is_aa_enabled() const { return _antialias; }
    void add_non_tex_polygon(lb_polygon* poly, float z, uint brush_tag);
    bat_batch* add_tex_polygons(lb_polygon_list& polys, float z);
    bat_line* add_line(lb_joint* i, lb_joint* j, float w, float z, uint pen_tag);
    bat_line* add_aa_border(lb_joint* i, lb_joint* j, float z, uint pen_tag);
    void finish_batching();
    bat_batches& get_batches() { return _batches; }
    void clear_batches();

protected:
    bat_triangles       _triangles;
    bat_lines           _lines;
    bat_batches         _batches;
    bool                _antialias;

protected:
    template<class _batch>
    _batch* create_batch(bat_type t);
    bat_triangle* create_triangle(lb_joint* i, lb_joint* j, lb_joint* k, float z);
    bat_line* create_line(lb_joint* i, lb_joint* j, float w, float z, uint t, bool half);
    bat_line* create_half_line(lb_joint* i, lb_joint* j, const vec2& p1, const vec2& p2, float w, float z, uint t);
    void add_triangle(lb_joint* i, lb_joint* j, lb_joint* k, bool b[3], float z, uint brush_tag);
    void collect_aa_borders(bat_triangle* triangle, bool b[3], uint pen_tag);
    void proceed_line_batch();
    void gather_tex_triangles(bat_triangles& triangles, lb_polygon* poly, float z);
    bat_batch* find_containable_tex_batch(const bat_triangles& triangles);
    bat_stroke_batch* find_associated_tex_stroke_batch(const bat_batch* bat);
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/**
 *
 * @file core_dgeqp3_norms.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include <math.h>
#include <cblas.h>
#include "common.h"

#define A(m,n) BLKADDR( A, double, m, n )

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dgeqp3_norms computes the 2-norm of each column of A[ ioff:m, joff:n ]
 *  that is marked with norms2[j] == -1 on entry. Entries that are not marked
 *  are assumed to already contain the correct 2-norm, so that the same routine
 *  can be used for computing the initial norms and for updating bad norms.
 *  The result is stored duplicated in norms1 and norms2.
 *
 *******************************************************************************
 *
 *  @param[in] A
 *          PLASMA descriptor of the matrix A.
 *          On entry, the M-by-N matrix described by the descriptor.
 *
 *  @param[in] ioff
 *          Row offset.
 *
 *  @param[in] joff
 *          Column offset.
 *
 *  @param[in,out] norms1
 *          Vector of size A.n.
 *          On exit, norms1[j] is 2-norm of column j, for j >= joff.
 *
 *  @param[in,out] norms2
 *          Vector of size A.n.
 *          On entry, if norms2[j] == -1, re-compute norm of column j.
 *          On exit, norms2[j] is 2-norm of column j, for j >= joff.
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgeqp3_norms = PCORE_dgeqp3_norms
#define CORE_dgeqp3_norms PCORE_dgeqp3_norms
#define CORE_dgessq PCORE_dgessq
int
CORE_dgessq(int M, int N,
            const double *A, int LDA,
            double *scale, double *sumsq);
#endif
void CORE_dgeqp3_norms( PLASMA_desc A, int ioff, int joff, double *norms1, double *norms2 )
{
    const double *Ai;
    int j, ii, ioff2, len, mb, nb, lda;
    double sumsq, scale;

    if ( A.nt != 1 ) {
        coreblas_error(1, "Illegal value of A.nt");
        return;
    }

    nb = min( A.nb, A.n );
    for( j = joff; j < nb; ++j ) {
        if ( norms2[j] == -1. ) {
            scale = 0.;
            sumsq = 1.;
            ioff2 = ioff;
            for( ii = 0; ii < A.mt; ++ii ) {
                mb = min( A.mb, A.m - ii*A.mb );
                Ai = A(ii,0);
                lda = BLKLDD( A, ii );
                len = mb - ioff2;
                CORE_dgessq( len, 1, Ai + j*lda + ioff2, lda, &scale, &sumsq );
                ioff2 = 0;
            }
            norms2[j] = scale * sqrt( sumsq );
            norms1[j] = norms2[j];
        }
    }
}


text:
#pragma once

#include <imageview/ImageRowView.h>
#include <imageview/IsPixelFormat.h>
#include <imageview/internal/ImageViewIterator.h>
#include <imageview/internal/ImageViewStorage.h>
#include <imageview/internal/PixelRef.h>

#include <gsl/assert>
#include <gsl/span>

#include <cstddef>
#include <type_traits>

namespace imageview {

// Non-owning view into a bitmap image with the specified pixel format.
//
// \param PixelFormat - specifies how colors are stored in the bitmap.
//        ContinuousImageView only supports images with fixed color depth (bits
//        per pixel). PixelFormat should satisfy IsPixelFormat trait, i.e. be
//        like
//          class MyPixelFormat {
//           public:
//            using color_type = MyColor;
//            static constexpr int kBytesPerPixel = N;
//            color_type read(gsl::span<const std::byte, kBytesPerPixel> pixel_data) const;
//            void write(const color_type& color, gsl::span<std::byte, kBytesPerPixel> pixel_data) const;
//          };
// \param Mutable - if true, ContinuousImageView provides write access to the
//        bitmap (naturally, this requires that ContinuousImageView is
//        constructed from a non-const pointer to the data). Otherwise, only
//        read-only access is provided.
template <class PixelFormat, bool Mutable = false>
class ContinuousImageView {
 public:
  static_assert(IsPixelFormat<PixelFormat>::value, "Not a PixelFormat.");

  using byte_type = std::conditional_t<Mutable, std::byte, const std::byte>;
  using value_type = typename PixelFormat::color_type;
  // If `Mutable == false`, then `reference` is an alias to `value_type`.
  // Otherwise, it is a proxy class, which mimics `value_type&`:
  // * is implicitly convertible to `value_type`.
  // * you can assign a value_type to it, changing the color of the referenced pixel.
  using reference = std::conditional_t<Mutable, detail::PixelRef<PixelFormat>, value_type>;
  // Constant LegacyInputIterator whose value_type is @value_type.
  // It is not, however, a LegacyForwardIterator, because
  //   std::iterator_traits<iterator>::reference
  // is not const value_type&.
  // Nevertheless, this iterator supports all arithmetic operations of LegacyRandomAccessIterator, e.g.,
  //   `it += n`, ` --it`, `it[n]`.
  using const_iterator = detail::ImageViewIterator<PixelFormat, false>;
  // LegacyInputIterator whose value_type is @value_type. If Mutable == true, then it also
  // models LegacyOutputIterator. It is not, however, a LegacyForwardIterator, because
  //   std::iterator_traits<iterator>::reference
  // is neither value_type& nor const value_type&.
  // Nevertheless, this iterator supports all arithmetic operations of LegacyRandomAccessIterator, e.g.,
  //   `it += n`, ` --it`, `it[n]`.
  using iterator = detail::ImageViewIterator<PixelFormat, Mutable>;

  // Construct an empty view.
  // This constructor is only available if PixelFormat is default-constructible.
  template <class Enable = std::enable_if_t<std::is_default_constructible_v<PixelFormat>>>
  constexpr ContinuousImageView() noexcept(noexcept(std::is_nothrow_default_constructible_v<PixelFormat>)) {}

  // Construct a view into an image.
  // This constructor is only available if PixelFormat is default-constructible
  // \param height - height of the image.
  // \param width - width of the image.
  // \param data - pointer to the uncompressed RGB24 bitmap data. The size of
  //        the array should be exactly
  //          height * width * PixelFormat::kBytesPerPixel.
  //        Pixels are assumed to be stored contiguously, with each pixel
  //        occupying exactly PixelFormat::kBytesPerPixel.
  template <class Enable = std::enable_if_t<std::is_default_constructible_v<PixelFormat>>>
  constexpr ContinuousImageView(unsigned int height, unsigned int width, gsl::span<byte_type> data) noexcept(
      std::is_nothrow_default_constructible_v<PixelFormat>);

  // Construct a view into an image.
  // \param height - height of the image.
  // \param width - width of the image.
  // \param data - pointer to the uncompressed RGB24 bitmap data. The size of
  //        the array should be exactly
  //          height * width * PixelFormat::kBytesPerPixel.
  //        Pixels are assumed to be stored contiguously, with each pixel
  //        occupying exactly PixelFormat::kBytesPerPixel.
  // \param pixel_format - PixelFormat instance to use.
  constexpr ContinuousImageView(unsigned int height, unsigned int width, gsl::span<byte_type> data,
                                const PixelFormat& pixel_format);

  // Construct a view into an image.
  // \param height - height of the image.
  // \param width - width of the image.
  // \param data - pointer to the uncompressed RGB24 bitmap data. The size of
  //        the array should be exactly
  //          height * width * PixelFormat::kBytesPerPixel.
  //        Pixels are assumed to be stored contiguously, with each pixel
  //        occupying exactly PixelFormat::kBytesPerPixel.
  // \param pixel_format - PixelFormat instance to use.
  constexpr ContinuousImageView(unsigned int height, unsigned int width, gsl::span<byte_type> data,
                                PixelFormat&& pixel_format);

  // Construct a read-only view from a mutable view.
  template <class Enable = std::enable_if_t<!Mutable>>
  constexpr ContinuousImageView(ContinuousImageView<PixelFormat, !Mutable> image);

  // Returns the height of the image.
  constexpr unsigned int height() const noexcept;

  // Returns the width of the image.
  constexpr unsigned int width() const noexcept;

  // Returns the total number of pixels.
  constexpr std::size_t area() const noexcept;

  // Returns true if the image has zero area, false otherwise.
  constexpr bool empty() const noexcept;

  // Returns the pixel format used by this image.
  constexpr const PixelFormat& pixelFormat() const noexcept;

  // Returns the pointer to the bitmap data.
  constexpr gsl::span<byte_type> data() const noexcept;

  // Returns the total size of the bitmap in bytes.
  constexpr std::size_t size_bytes() const noexcept;

  // Returns an iterator to the first pixel.
  constexpr iterator begin() const;

  // Returns a const iterator to the first pixel.
  constexpr const_iterator cbegin() const;

  // Returns an iterator past the last pixel.
  constexpr iterator end() const;

  // Returns a const iterator past the last pixel.
  constexpr const_iterator cend() const;

  // Access the specified pixel.
  // \param y - Y coordinate of the pixel. Should be within [0; height()).
  // \param x - X coordinate of the pixel. Should be within [0; width()).
  // \return the color of the specified pixel.
  constexpr reference operator()(unsigned int y, unsigned int x) const;

  // Returns a non-owning view into the specified row of the image.
  // \param y - 0-based index of the row. Should be within [0; height()).
  // \return a non-owning view into the row @y.
  constexpr ImageRowView<PixelFormat, Mutable> row(unsigned int y) const;

 private:
  constexpr gsl::span<byte_type, PixelFormat::kBytesPerPixel> getPixelData(unsigned int y, unsigned int x) const;

  detail::ImageViewStorage<PixelFormat, Mutable> storage_;
  unsigned int height_ = 0;
  unsigned int width_ = 0;
};

// Convert a ContinuousImageView into an ImageRowView.
// \param image - input image.
// \return an ImageRowView referring to the same data as image; the number of
//         elements in the returned view equals image.area().
template <class PixelFormat, bool Mutable>
constexpr ImageRowView<PixelFormat, Mutable> flatten(ContinuousImageView<PixelFormat, Mutable> image) {
  return ImageRowView<PixelFormat, Mutable>(image.data(), image.area(), image.pixelFormat());
}

template <class PixelFormat, bool Mutable>
template <class Enable>
constexpr ContinuousImageView<PixelFormat, Mutable>::ContinuousImageView(
    unsigned int height, unsigned int width,
    gsl::span<byte_type> data) noexcept(std::is_nothrow_default_constructible_v<PixelFormat>)
    : storage_(data.data()), height_(height), width_(width) {
  Expects(data.size() == height * width * PixelFormat::kBytesPerPixel);
}

template <class PixelFormat, bool Mutable>
constexpr ContinuousImageView<PixelFormat, Mutable>::ContinuousImageView(unsigned int height, unsigned int width,
                                                                         gsl::span<byte_type> data,
                                                                         const PixelFormat& pixel_format)
    : storage_(data.data(), pixel_format), height_(height), width_(width) {
  Expects(data.size() == height * width * PixelFormat::kBytesPerPixel);
}

template <class PixelFormat, bool Mutable>
constexpr ContinuousImageView<PixelFormat, Mutable>::ContinuousImageView(unsigned int height, unsigned int width,
                                                                         gsl::span<byte_type> data,
                                                                         PixelFormat&& pixel_format)
    : storage_(data.data(), std::move(pixel_format)), height_(height), width_(width) {
  Expects(data.size() == height * width * PixelFormat::kBytesPerPixel);
}

template <class PixelFormat, bool Mutable>
template <class Enable>
constexpr ContinuousImageView<PixelFormat, Mutable>::ContinuousImageView(
    ContinuousImageView<PixelFormat, !Mutable> image)
    : ContinuousImageView(image.height(), image.width(), image.data(), image.pixelFormat()) {}

template <class PixelFormat, bool Mutable>
constexpr unsigned int ContinuousImageView<PixelFormat, Mutable>::height() const noexcept {
  return height_;
}

template <class PixelFormat, bool Mutable>
constexpr unsigned int ContinuousImageView<PixelFormat, Mutable>::width() const noexcept {
  return width_;
}

template <class PixelFormat, bool Mutable>
constexpr std::size_t ContinuousImageView<PixelFormat, Mutable>::area() const noexcept {
  return static_cast<std::size_t>(height_) * width_;
}

template <class PixelFormat, bool Mutable>
constexpr bool ContinuousImageView<PixelFormat, Mutable>::empty() const noexcept {
  return height_ == 0 || width_ == 0;
}

template <class PixelFormat, bool Mutable>
constexpr const PixelFormat& ContinuousImageView<PixelFormat, Mutable>::pixelFormat() const noexcept {
  return storage_.pixelFormat();
}

template <class PixelFormat, bool Mutable>
constexpr auto ContinuousImageView<PixelFormat, Mutable>::data() const noexcept -> gsl::span<byte_type> {
  return gsl::span<byte_type>(storage_.data(), size_bytes());
}

template <class PixelFormat, bool Mutable>
constexpr std::size_t ContinuousImageView<PixelFormat, Mutable>::size_bytes() const noexcept {
  return area() * PixelFormat::kBytesPerPixel;
}

template <class PixelFormat, bool Mutable>
constexpr auto ContinuousImageView<PixelFormat, Mutable>::begin() const -> iterator {
  return iterator(storage_.data(), pixelFormat());
}

template <class PixelFormat, bool Mutable>
constexpr auto ContinuousImageView<PixelFormat, Mutable>::cbegin() const -> const_iterator {
  return const_iterator(storage_.data(), pixelFormat());
}

template <class PixelFormat, bool Mutable>
constexpr auto ContinuousImageView<PixelFormat, Mutable>::end() const -> iterator {
  return iterator(storage_.data() + size_bytes(), pixelFormat());
}

template <class PixelFormat, bool Mutable>
constexpr auto ContinuousImageView<PixelFormat, Mutable>::cend() const -> const_iterator {
  return const_iterator(storage_.data() + size_bytes(), pixelFormat());
}

template <class PixelFormat, bool Mutable>
constexpr auto ContinuousImageView<PixelFormat, Mutable>::getPixelData(unsigned int y, unsigned int x) const
    -> gsl::span<byte_type, PixelFormat::kBytesPerPixel> {
  Expects(y < height_);
  Expects(x < width_);
  const std::size_t offset = (y * width_ + x) * PixelFormat::kBytesPerPixel;
  return gsl::span<byte_type, PixelFormat::kBytesPerPixel>(storage_.data() + offset, PixelFormat::kBytesPerPixel);
}

template <class PixelFormat, bool Mutable>
constexpr auto ContinuousImageView<PixelFormat, Mutable>::operator()(unsigned int y, unsigned int x) const
    -> reference {
  const gsl::span<byte_type, PixelFormat::kBytesPerPixel> pixel_data = getPixelData(y, x);
  if constexpr (Mutable) {
    return detail::PixelRef<PixelFormat>(pixel_data, storage_.pixelFormat());
  } else {
    return storage_.pixelFormat().read(pixel_data);
  }
}

template <class PixelFormat, bool Mutable>
constexpr ImageRowView<PixelFormat, Mutable> ContinuousImageView<PixelFormat, Mutable>::row(unsigned int y) const {
  Expects(y < height_);
  const std::size_t bytes_per_row = static_cast<std::size_t>(width_) * PixelFormat::kBytesPerPixel;
  const gsl::span<byte_type> row_data(storage_.data() + y * bytes_per_row, bytes_per_row);
  return ImageRowView<PixelFormat, Mutable>(row_data, width_, pixelFormat());
}

}  // namespace imageview

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview