text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_interp.h>

#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spce_PET.h"
#include "inout_aper.h"
#include "trace_conf.h"
#include "aper_conf.h"
#include "spc_sex.h"
#include "disp_conf.h"
#include "spc_wl_calib.h"
#include "spce_binning.h"
#include "spc_optimum.h"
#include "spc_spc.h"
#include "fringe_conf.h"
#include "spc_resp.h"
#include "spc_FITScards.h"

#define AXE_IMAGE_PATH   "AXE_IMAGE_PATH"
#define AXE_OUTPUT_PATH  "AXE_OUTPUT_PATH"
#define AXE_DRIZZLE_PATH "AXE_DRIZZLE_PATH"
#define AXE_CONFIG_PATH  "AXE_CONFIG_PATH"

int
main (int argc, char *argv[])
{
  char *opt;

  char grism_file[MAXCHAR];
  char grism_file_path[MAXCHAR];

  char aper_file[MAXCHAR];
  char aper_file_path[MAXCHAR];

  char obj_PET_file[MAXCHAR];
  char obj_PET_file_path[MAXCHAR];

  char bck_PET_file[MAXCHAR];
  char bck_PET_file_path[MAXCHAR];

  char conf_file[MAXCHAR];
  char conf_file_path[MAXCHAR];

  char SPC_file[MAXCHAR];
  char SPC_file_path[MAXCHAR];

  char SPC_opt_file[MAXCHAR];
  char SPC_opt_file_path[MAXCHAR];

  char WHT_file[MAXCHAR];
  char WHT_file_path[MAXCHAR];

  char label[MAXCHAR];

  int i, index, dobck = 0, noflux = 1;

  object **oblist;

  FITScards      *cards;

  ap_pixel *obj_PET = NULL, *bck_PET = NULL;

  observation *obs;

  //tracestruct *trace;
  aperture_conf *conf;

  spectrum *obj_spec = NULL, *bck_spec = NULL, *sobj_spec = NULL;
  spectrum *resp;
  response_function *resp_func;
  calib_function *wl_calibration;

  fitsfile *OPET_ptr, *BPET_ptr;
  int f_status = 0;

  fitsfile *SPC_ptr, *SPC_opt_ptr, *WHT_ptr;

  gsl_matrix   *weights;
  drzstamp     *modvar;

  int obj_aperID, obj_beamID, objindex;
  int bck_aperID, bck_beamID;

  char table[MAXCHAR], table_path[MAXCHAR];
  //char comment[FLEN_COMMENT];
  int empty;
  int drizzle;
  int quant_cont=0;
  int opt_weights=0;
  int for_grism=0;
  int smooth_conv=0;

  d_point smooth_params;

  double exptime;
  double sky_cps;

  drzstamp_dim  dimension;
  gsl_matrix *coverage = NULL;

  if (((argc < 3))
      || (opt = get_online_option ("help", argc, argv)))
    {
      fprintf (stdout,
	       "aXe_PET2SPC Version %s:\n"
	       "             aXe task that produces 1-D, binned spectra using information\n"
	       "             contained in a OPET and a BPET. A 1-D spectrum is generated\n"
	       "             for each beam in each of both the OPET and the BPET files\n"
	       "             The binned background spectra of each beam (order) is then\n"
	       "             subtraced from the corresponding spectra form the OPET. Th\ne"
	       "             background subtraction (and reading a BPET altogether can be\n"
	       "             avoided by using the -noBPET option\n"
	       "             An SPC file, a multi-extension FITS file containing binned\n"
	       "             spectra is produced. Each extension (named after the beam ID,\n"
	       "             e.g. 11B for aperture (object) 11,beam (order) B) contains the\n"
	       "             following columns:\n"
	       "\n"
	       "              LAMBDA      ; the wavelength (in A)\n"
	       "              TCOUNT      ; the total number of counts (in DN)\n"
	       "              TERROR      ; the error in TERRORS (in DN)\n"
	       "              BCOUNT      ; the estimated number of counts from the\n"
	       "                            background (in DN)\n"
	       "              BERROR      ; the error in BCOUNTS (in DN)\n"
	       "              COUNT       ; the estimated number of counts from the\n"
	       "                            object (in DN)\n"
	       "              ERROR       ; the error in COUNTS (in DN)\n"
	       "              FLUX        ; the estimated flux (in erg/s/cm^2/A)\n"
	       "              FERROR      ; the error in FLUX (in erg/s/cm^s/A)\n"
	       "              WEIGHT      ; weight (in pixels)\n"
	       "\n"
	       "             Input FITS mages are looked for in $AXE_IMAGE_PATH\n"
	       "             aXe config file is looked for in $AXE_CONFIG_PATH\n"
	       "             All outputs are writen to $AXE_OUTPUT_PATH\n"
	       "\n"
	       "Usage:\n"
	       "     aXe_PET2SPC [g/prism image filename] [aXe config file name] [options]\n"
	       "\n"
	       "Options:\n"
	       "             -noBPET         - to disable the use of a BPET file\n"
	       "             -noflux         - to disable the flux calibration\n"
	       "             -drz            - use $AXE_DRIZZLE_PATH to locate the grism-, OAF-\n"
	       "                             - and PET-files instead of $AXE_IMAGE/OUTPUT_PATH\n"
	       "             -in_AF=[string] - overwrite the default input Aperture File name\n"
	       "             -OPET=[string]  - overwrite the default input Object PET file name\n"
	       "             -BPET=[string]  - overwrite the default input Background PET\n"
	       "                               file name\n"
	       "             -out_SPC=[string] - overwrite the default output SPC file name\n"
	       "\n"
	       "Example:\n"
	       "       ./aXe_PET2SPC slim_grism.fits SLIM.conf.A.0\n"
	       "\n",RELEASE);
      exit (1);
    }

  // make a general opening statement
  fprintf (stdout, "aXe_PET2SPC: Starting...\n");

  // get the name of the flt-file
  index = 0;
  strcpy (grism_file, argv[++index]);
  if ((opt = get_online_option ("drz", argc, argv)))
    {
      build_path (AXE_DRIZZLE_PATH, grism_file, grism_file_path);
      drizzle=1;
    }
  else
    {
      build_path (AXE_IMAGE_PATH, grism_file, grism_file_path);
      drizzle=0;
    }

  // get the name of the configuration file
  strcpy (conf_file, argv[++index]);
  build_path (AXE_CONFIG_PATH, conf_file, conf_file_path);

  // load the configuration file
  conf = get_aperture_descriptor (conf_file_path);

  // Determine where the various extensions are in the FITS file
  get_extension_numbers(grism_file_path, conf,conf->optkey1,conf->optval1);

  // Build aperture file name
  replace_file_extension (grism_file, aper_file, ".fits",
			  ".OAF", conf->science_numext);
  if (drizzle)
    build_path (AXE_DRIZZLE_PATH, aper_file, aper_file_path);
  else
    build_path (AXE_OUTPUT_PATH, aper_file, aper_file_path);

  // Build object PET file name
  replace_file_extension (grism_file, obj_PET_file, ".fits",
			  ".PET.fits", conf->science_numext);
  // make the total filename
  if (drizzle)
    build_path (AXE_DRIZZLE_PATH, obj_PET_file, obj_PET_file_path);
  else
    build_path (AXE_OUTPUT_PATH, obj_PET_file, obj_PET_file_path);

  // Build background PET file name
  replace_file_extension (grism_file, bck_PET_file, ".fits",
			  ".BCK.PET.fits", conf->science_numext);

  // make the total filename
  if (drizzle)
    build_path (AXE_DRIZZLE_PATH, bck_PET_file, bck_PET_file_path);
  else
    build_path (AXE_OUTPUT_PATH, bck_PET_file, bck_PET_file_path);

  // make a non-standard AF name if necessary
  if ((opt = get_online_option ("in_AF", argc, argv)))
    {
      strcpy(aper_file,opt);
      strcpy(aper_file_path,opt);
    }

  // make a non-standard PET name if necessary
  if ((opt = get_online_option ("OPET", argc, argv)))
    {
      strcpy(obj_PET_file,opt);
      strcpy(obj_PET_file_path,opt);
    }

  // make a non-standard BPET name if necessary
  if ((opt = get_online_option ("BPET", argc, argv)))
    {
      strcpy(bck_PET_file,opt);
      strcpy(bck_PET_file_path,opt);
    }

  // set the flagg for no background subtraction
  if ((opt = get_online_option ("noBPET",argc,argv)))
    {
      dobck = 0;
      strcpy(bck_PET_file,"None");
      strcpy(bck_PET_file_path, "None");
    }
  else
    {
      dobck = 1;
    }

  // set the flagg for flux
  if ((opt = get_online_option ("noflux",argc,argv)))
    noflux = 1;
  else
    noflux = 0;

  // check for the weights flagg,
  // set the file name if the flagg is set
  if ((opt = get_online_option ("opt_weights",argc,argv)))
    opt_weights = 1;
  else
    opt_weights = 0;

  // read the trigger for smoothing the sensitivity
  if ((opt = get_online_option ("smooth_conv",argc,argv)))
    smooth_conv = 1;
  else
    smooth_conv = 0;

  if ((opt = get_online_option ("out_SPC", argc, argv)))
    {
      strcpy (SPC_file, opt);
      strcpy (SPC_file_path, opt);

      if (opt_weights)
	{
	  replace_file_extension (SPC_file, SPC_opt_file, ".fits",
			      "_opt.SPC.fits", -1);
	  replace_file_extension (SPC_file_path, SPC_opt_file_path, ".fits",
			      "_opt.SPC.fits", -1);
	  replace_file_extension (SPC_opt_file, WHT_file, ".SPC.fits",
				  ".WHT.fits", -1);
	  replace_file_extension (SPC_opt_file_path, WHT_file_path, ".SPC.fits",
				  ".WHT.fits", -1);
	}
    }
  else
    {
      replace_file_extension (obj_PET_file, SPC_file, ".PET.fits",
			      ".SPC.fits", -1);
      if (drizzle)
	build_path (AXE_DRIZZLE_PATH, SPC_file, SPC_file_path);
      else
	build_path (AXE_OUTPUT_PATH, SPC_file, SPC_file_path);

      if (opt_weights)
	{
	  replace_file_extension (obj_PET_file, SPC_opt_file, ".PET.fits",
				  "_opt.SPC.fits", -1);
	  replace_file_extension (SPC_opt_file, WHT_file, ".SPC.fits",
				  ".WHT.fits", -1);
	  if (drizzle)
	    {
	      build_path (AXE_DRIZZLE_PATH, SPC_opt_file, SPC_opt_file_path);
	      build_path (AXE_DRIZZLE_PATH, WHT_file, WHT_file_path);
	    }
	  else
	    {
	      build_path (AXE_OUTPUT_PATH, SPC_opt_file, SPC_opt_file_path);
	      build_path (AXE_OUTPUT_PATH, WHT_file, WHT_file_path);
	    }
	}
    }

  // check the configuration file for the smoothin keywords
  if (!check_conf_for_smoothing(conf, smooth_conv))
    aXe_message(aXe_M_FATAL, __FILE__, __LINE__,
        "aXe_PET2SP: Either the configuration file %s does not contain\n"
        "the necessary keywords for the smoothing (POBJSIZE, SMFACTOR),\n"
        "or one of these keywords has an unreasonable value < 0.0!\n",
        conf_file_path);

  // give feedback onto the screen:
  // report on input and output
  // and also on specific parameter
  fprintf (stdout, "aXe_PET2SPC: Input configuration file name:   %s\n",
	   conf_file_path);
  fprintf (stdout, "aXe_PET2SPC: Input Object Aperture file name: %s\n",
	   aper_file_path);
  fprintf (stdout, "aXe_PET2SPC: Input Object PET file name:      %s\n",
	   obj_PET_file_path);
  if (dobck)
    fprintf (stdout, "aXe_PET2SPC: Input Background PET file name:  %s\n",
	     bck_PET_file_path);
  fprintf (stdout, "aXe_PET2SPC: Output SPC file name:            %s\n",
	   SPC_file_path);
  if (opt_weights)
    {
      fprintf (stdout, "aXe_PET2SPC: Computing optimal weights.\n");
      fprintf (stdout, "aXe_PET2SPC: Optimized SPC file name:         %s\n",
	       SPC_opt_file_path);
      fprintf (stdout, "aXe_PET2SPC: Output WHT file name:            %s\n",
	       WHT_file_path);
    }
  if (!noflux)
    {
      fprintf (stdout, "aXe_PET2SPC: Performing flux calibration.\n");
      if (smooth_conv)
	fprintf (stdout, "aXe_PET2SPC: Using smoothed sensitivity curves.\n");
    }
  fprintf (stdout, "\n\n");

  //
  // try to get the descriptor 'exptime' from the 'sci'-extension
  //
  exptime = (double)get_float_from_keyword(grism_file_path, conf->science_numext, conf->exptimekey);
  if (isnan(exptime))
    exptime = (double)get_float_from_keyword(grism_file_path, 1, conf->exptimekey);
  if (isnan(exptime))
    exptime = 1.0;

  //
  // try to get the descriptor 'SKY_CPS' from the 'sci'-extension
  //
  sky_cps = (double)get_float_from_keyword(grism_file_path, conf->science_numext, "SKY_CPS");
  if (isnan(sky_cps))
    sky_cps = 0.0;

  //  Open the OPET file for reading
  fits_open_file (&OPET_ptr, obj_PET_file_path, READONLY, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_PET2SPC: Could not open file: %s\n",
		   obj_PET_file_path);
    }

  // check whether the contamination is quantitative
  quant_cont = check_quantitative_contamination(OPET_ptr);

  if (opt_weights && !quant_cont)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "aXe_PET2SPC: The optimal extractions needs quantitative contamination! "
		 " Please re-run aXe with Gauss or Fluxcube contamination!");

  obs = load_dummy_observation ();
  /* Loading the object list */
  fprintf (stdout, "aXe_PET2SPC: Loading object aperture list...");
  oblist = file_to_object_list_seq (aper_file_path, obs);
  fprintf (stdout,"%d objects loaded.\n",object_list_size(oblist));

  if (dobck)
    {
      //  Open the file for reading
      fits_open_file (&BPET_ptr, bck_PET_file_path, READONLY, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_PET2SP: Could not open file: %s\n",
		       bck_PET_file_path);
	}
    }

  /* Copy the header info from the grism image */
  SPC_ptr = create_SPC_opened (SPC_file_path,1);
  cards = get_FITS_cards_opened(OPET_ptr);
  put_FITS_cards_opened(SPC_ptr, cards);
  free_FITScards(cards);

  if (opt_weights)
    {
      cards = get_FITS_cards_opened(OPET_ptr);

      // open the WHT file and add the header keywords
      WHT_ptr = create_FITSimage_opened (WHT_file_path, 1);
      put_FITS_cards_opened(WHT_ptr, cards);

      // open the opt_SPC and add the header keywords
      SPC_opt_ptr = create_SPC_opened (SPC_opt_file_path,1);
      put_FITS_cards_opened(SPC_opt_ptr, cards);

      // delete the header keywords
      free_FITScards(cards);
    }

  // do something only if there
  // exist valid objects
  i = 0;
  if (oblist!=NULL)
    {

      // turn until the end of the PET is reached
      while (1)
	{
	  empty=0;

	  // Get the PET for this object
	  obj_PET = get_ALL_from_next_in_PET(OPET_ptr, &obj_aperID, &obj_beamID);

	  // load the background PET if requested
	  if (dobck)
	    {
	      bck_PET = get_ALL_from_next_in_PET(BPET_ptr, &bck_aperID, &bck_beamID);
	      if ((bck_aperID!=obj_aperID)||(bck_beamID!=obj_beamID))
		aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			     "Background PET and Object PET extensions are not"
			     " in the same order and cannot be combined.\n");
	    }

	  // end of PET reached: the break condition
	  if ((obj_aperID==-1) && (obj_beamID==-1))
	    break;

	  // signal an empty PET
	  if (obj_PET==NULL)
	    empty=1;

	  // give feedback to the screen
	  fprintf (stdout, "aXe_PET2SPC: BEAM %d; %d%c\n", i, obj_aperID, BEAM(obj_beamID));fflush(stdout);

	  // identify the object which matches the PET
	  objindex =  find_object_in_object_list(oblist,obj_aperID);

	  // look whether we are for grisms or prisms
	  for_grism = check_for_grism (conf_file_path, obj_beamID);
	  wl_calibration  = get_calfunc_for_beam(oblist[objindex]->beams[obj_beamID], for_grism, conf_file_path, conf);

	  // compute the object spectrum
	  obj_spec = bin_naive (obj_PET, oblist[objindex]->beams[obj_beamID].width,
				oblist[objindex]->beams[obj_beamID].orient,	quant_cont);

	  // check for the existence of a background PET
	  if (dobck)
	    {
	      // compute the background spectrum
	      bck_spec = bin_naive (bck_PET,
				    oblist[objindex]->beams[bck_beamID].width,
				    oblist[objindex]->beams[bck_beamID].orient, quant_cont);
	    }
	  else
	    {
	      // create a dummy background spectrum
	      bck_spec = empty_counts_spectrum_copy(obj_spec);
	    }

	  // subtract the background spectrum from the
	  // object (or forground) spectrum
	  sobj_spec = subtract_spectra (obj_spec, bck_spec);

	  if(!noflux)
	    {
	      get_troughput_table_name(conf_file_path,
				       oblist[objindex]->beams[obj_beamID].ID,
				       table);
	      if (strcmp(table,"None"))
		{
		  build_path (AXE_CONFIG_PATH, table, table_path);
		  resp=get_response_function_from_FITS(table_path,2);
		  resp_func = create_response_function(table_path);
		  if (resp->spec_len <2)
		    {
		      aXe_message (aXe_M_WARN1, __FILE__, __LINE__,
				   "Throughput table %s contains only %d"
				   " values. No sensitivity curve was applied.\n",
				   table_path,resp->spec_len);
		    }
		  else
		    {
		      fprintf(stdout,"aXe_PET2SPC: Applying sensitivity contained in %s\n",table_path);
		      smooth_params = get_smooth_pars_for_beam(conf, smooth_conv, oblist[objindex]->beams[obj_beamID]);
		      if (smooth_params.x > 0.0)
		        {
		          // apply a smoothed flux conversion
			        apply_smoothed_response(wl_calibration, for_grism, quant_cont, resp_func, smooth_params, sobj_spec);
			       }
		      else
			   {
			       // apply a normal flux conversion
			       apply_response_function(sobj_spec, resp, quant_cont);
			   }
		    }
		  // free the memory of the
		  // response functions
		  free_spectrum(resp);
		  free_response_function(resp_func);
		}
	    }
	  if (empty!=1)
	    {
	      add_spectra_to_SPC_opened (SPC_ptr, obj_spec, bck_spec,
					 sobj_spec, oblist[objindex]->ID, oblist[objindex]->beams[obj_beamID].ID);

	      /* Copy header from OPET extension into this SPC extension */
	      cards = get_FITS_cards_opened(OPET_ptr);
	      put_FITS_cards_opened(SPC_ptr,cards);
	      free_FITScards(cards);

	    }

	  free_spectrum (bck_spec);
	  free_spectrum (sobj_spec);
	  free_spectrum (obj_spec);

	  if (opt_weights)
	    {
	      // get the dimension in trace length
	      // and crossdispersion
	      dimension = get_all_dims(obj_PET, bck_PET,
				       oblist[objindex]->beams[obj_beamID], dobck);

	      // check for empty PET
	      if (!dimension.resolution)
		{
		  // create dummies in case of empty PET's
		  weights = get_default_weight();
		  modvar  = get_default_modvar();
		}
	      else
		{
		  // prepare the PET's by computing the inverse variance.
		  // Also the trace distances are shifted by 0.5
		  // to get a sampling comparable to the unweighted
		  // extraction
		  prepare_inv_variance(obj_PET, bck_PET, dobck, conf, exptime, sky_cps, 0.0);
		  // compute the inverse variance and the profile
		  // image in the trace distance - crossdispersion plane
		  modvar = compute_modvar(obj_PET, oblist[objindex]->beams[obj_beamID], dimension);

		  // compute the optimal weights
		  weights = comp_allweight(modvar);

		}

	      if (dimension.resolution && empty != 1)
		{
		  sprintf (label, "WHT_%d%c", obj_aperID, BEAM (obj_beamID));
		  gsl_to_FITSimage_opened (weights, WHT_ptr ,0,label);

		  // make and store the default header
		  cards = beam_to_FITScards(oblist[objindex],obj_beamID);
		  put_FITS_cards_opened(WHT_ptr,cards);
		  free_FITScards(cards);
		}

	      // create the optimal weighted
	      // foreground spectrum
	      obj_spec = bin_optimal (obj_PET,oblist[objindex]->beams[obj_beamID],
				      quant_cont, weights, dimension, coverage);

	      // check for the presence of a background PET
	      if (dobck)
		{
		  // create the optimal weighted
		  // background spectrum
		  bck_spec = bin_optimal (bck_PET,oblist[objindex]->beams[obj_beamID],
					  quant_cont, weights, dimension, coverage);
		}
	      else
		{
		  // make a dummy background spectrum
		  bck_spec = empty_counts_spectrum_copy(obj_spec);
		}

	      // release memory
	      gsl_matrix_free(weights);
	      free_drzstamp(modvar);

	  // subtract the background spectrum from the
	  // object (or forground) spectrum
	  sobj_spec = subtract_spectra (obj_spec, bck_spec);

	  if(!noflux)
	    {
	      get_troughput_table_name(conf_file_path,
				       oblist[objindex]->beams[obj_beamID].ID,
				       table);
	    if (strcmp(table,"None"))
	      {
		build_path (AXE_CONFIG_PATH, table, table_path);
		resp=get_response_function_from_FITS(table_path,2);
		resp_func = create_response_function(table_path);
		if (resp->spec_len <2)
		  {
		    aXe_message (aXe_M_WARN1, __FILE__, __LINE__,
				 "Throughput table %s contains only %d"
				 " values. No sensitivity curve was applied.\n",
				 table_path,resp->spec_len);
		  }
		else
		  {
		    fprintf(stdout,"aXe_PET2SPC: Applying sensitivity contained in %s\n",table_path);
		    smooth_params = get_smooth_pars_for_beam(conf, smooth_conv, oblist[objindex]->beams[obj_beamID]);
        if (smooth_params.x > 0.0)
          {
            apply_smoothed_response(wl_calibration, for_grism, quant_cont, resp_func, smooth_params, sobj_spec);
          }
		    else
		      {
			apply_response_function(sobj_spec, resp, quant_cont);
		      }
		  }
		// free the memory
		free_spectrum(resp);
		free_response_function(resp_func);
	      }
	    }
	  if (empty!=1)
	    {
	      add_spectra_to_SPC_opened (SPC_opt_ptr, obj_spec, bck_spec,
					 sobj_spec, oblist[objindex]->ID, oblist[objindex]->beams[obj_beamID].ID);

	      /* Copy header from OPET extension into this SPC extension */
	      cards = get_FITS_cards_opened(OPET_ptr);
	      put_FITS_cards_opened(SPC_opt_ptr,cards);
	      free_FITScards(cards);

	    }

	  free_spectrum (bck_spec);
	  free_spectrum (sobj_spec);
	  free_spectrum (obj_spec);
	    }

	  // free memory
	  free_calib (wl_calibration);
	  if (bck_PET!=NULL)
	    free (bck_PET);
	  if (obj_PET!=NULL)
	    free (obj_PET);
	  i++;
        }
    }

  fits_close_file (SPC_ptr, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_PET2SPC: " "Error closing SPC: %s\n",
		   SPC_file_path);
    }

  if (opt_weights)
    {
      fits_close_file (WHT_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_PET2SPC: " "Error closing WHT: %s\n",
		       WHT_file_path);
	}
      fits_close_file (SPC_opt_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_PET2SPC: " "Error closing PSC: %s\n",
		       SPC_opt_file_path);
	}
    }

  if (oblist!=NULL)
    free_oblist (oblist);

  fprintf (stdout, "aXe_PET2SPC: Done...\n");
  exit (0);
}


text:
/**
 *
 * @file qwrapper_claset2.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:59 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_claset2(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int M, int N,
                       PLASMA_Complex32_t alpha, PLASMA_Complex32_t *A, int LDA)
{
    DAG_CORE_LASET;
    QUARK_Insert_Task(quark, CORE_claset2_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(PLASMA_Complex32_t),         &alpha, VALUE,
        sizeof(PLASMA_Complex32_t)*M*N,     A,      OUTPUT,
        sizeof(int),                        &LDA,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_claset2_quark = PCORE_claset2_quark
#define CORE_claset2_quark PCORE_claset2_quark
#endif
void CORE_claset2_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    PLASMA_Complex32_t alpha;
    PLASMA_Complex32_t *A;
    int LDA;

    quark_unpack_args_6(quark, uplo, M, N, alpha, A, LDA);
    CORE_claset2(uplo, M, N, alpha, A, LDA);
}

ext: c


text:
/* matrix/gsl_matrix_ushort.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_USHORT_H__
#define __GSL_MATRIX_USHORT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_ushort.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned short * data;
  gsl_block_ushort * block;
  int owner;
} gsl_matrix_ushort;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_view;

typedef _gsl_matrix_ushort_view gsl_matrix_ushort_view;

typedef struct
{
  gsl_matrix_ushort matrix;
} _gsl_matrix_ushort_const_view;

typedef const _gsl_matrix_ushort_const_view gsl_matrix_ushort_const_view;

/* Allocation */

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_block (gsl_block_ushort * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_ushort *
gsl_matrix_ushort_alloc_from_matrix (gsl_matrix_ushort * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_ushort *
gsl_vector_ushort_alloc_row_from_matrix (gsl_matrix_ushort * m,
                                        const size_t i);

GSL_FUN gsl_vector_ushort *
gsl_vector_ushort_alloc_col_from_matrix (gsl_matrix_ushort * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_ushort_free (gsl_matrix_ushort * m);

/* Views */

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_submatrix (gsl_matrix_ushort * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_row (gsl_matrix_ushort * m, const size_t i);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_column (gsl_matrix_ushort * m, const size_t j);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_diagonal (gsl_matrix_ushort * m);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_subdiagonal (gsl_matrix_ushort * m, const size_t k);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_superdiagonal (gsl_matrix_ushort * m, const size_t k);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_subrow (gsl_matrix_ushort * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_ushort_view
gsl_matrix_ushort_subcolumn (gsl_matrix_ushort * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_array (unsigned short * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_array_with_tda (unsigned short * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector (gsl_vector_ushort * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_ushort_view
gsl_matrix_ushort_view_vector_with_tda (gsl_vector_ushort * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_submatrix (const gsl_matrix_ushort * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_row (const gsl_matrix_ushort * m,
                            const size_t i);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_column (const gsl_matrix_ushort * m,
                               const size_t j);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_diagonal (const gsl_matrix_ushort * m);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subdiagonal (const gsl_matrix_ushort * m,
                                    const size_t k);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_superdiagonal (const gsl_matrix_ushort * m,
                                      const size_t k);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subrow (const gsl_matrix_ushort * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_ushort_const_view
gsl_matrix_ushort_const_subcolumn (const gsl_matrix_ushort * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array (const unsigned short * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_array_with_tda (const unsigned short * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector (const gsl_vector_ushort * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_ushort_const_view
gsl_matrix_ushort_const_view_vector_with_tda (const gsl_vector_ushort * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_ushort_set_zero (gsl_matrix_ushort * m);
GSL_FUN void gsl_matrix_ushort_set_identity (gsl_matrix_ushort * m);
GSL_FUN void gsl_matrix_ushort_set_all (gsl_matrix_ushort * m, unsigned short x);

GSL_FUN int gsl_matrix_ushort_fread (FILE * stream, gsl_matrix_ushort * m) ;
GSL_FUN int gsl_matrix_ushort_fwrite (FILE * stream, const gsl_matrix_ushort * m) ;
GSL_FUN int gsl_matrix_ushort_fscanf (FILE * stream, gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_fprintf (FILE * stream, const gsl_matrix_ushort * m, const char * format);

GSL_FUN int gsl_matrix_ushort_memcpy(gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);
GSL_FUN int gsl_matrix_ushort_swap(gsl_matrix_ushort * m1, gsl_matrix_ushort * m2);

GSL_FUN int gsl_matrix_ushort_swap_rows(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ushort_swap_columns(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ushort_swap_rowcol(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ushort_transpose (gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_transpose_memcpy (gsl_matrix_ushort * dest, const gsl_matrix_ushort * src);

GSL_FUN unsigned short gsl_matrix_ushort_max (const gsl_matrix_ushort * m);
GSL_FUN unsigned short gsl_matrix_ushort_min (const gsl_matrix_ushort * m);
GSL_FUN void gsl_matrix_ushort_minmax (const gsl_matrix_ushort * m, unsigned short * min_out, unsigned short * max_out);

GSL_FUN void gsl_matrix_ushort_max_index (const gsl_matrix_ushort * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_ushort_min_index (const gsl_matrix_ushort * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_ushort_minmax_index (const gsl_matrix_ushort * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_ushort_isnull (const gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_ispos (const gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_isneg (const gsl_matrix_ushort * m);
GSL_FUN int gsl_matrix_ushort_isnonneg (const gsl_matrix_ushort * m);

GSL_FUN int gsl_matrix_ushort_add (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_sub (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_mul_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_div_elements (gsl_matrix_ushort * a, const gsl_matrix_ushort * b);
GSL_FUN int gsl_matrix_ushort_scale (gsl_matrix_ushort * a, const double x);
GSL_FUN int gsl_matrix_ushort_add_constant (gsl_matrix_ushort * a, const double x);
GSL_FUN int gsl_matrix_ushort_add_diagonal (gsl_matrix_ushort * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_ushort_get_row(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t i);
GSL_FUN int gsl_matrix_ushort_get_col(gsl_vector_ushort * v, const gsl_matrix_ushort * m, const size_t j);
GSL_FUN int gsl_matrix_ushort_set_row(gsl_matrix_ushort * m, const size_t i, const gsl_vector_ushort * v);
GSL_FUN int gsl_matrix_ushort_set_col(gsl_matrix_ushort * m, const size_t j, const gsl_vector_ushort * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL unsigned short   gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x);
GSL_FUN INLINE_DECL unsigned short * gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const unsigned short * gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
unsigned short
gsl_matrix_ushort_get(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_ushort_set(gsl_matrix_ushort * m, const size_t i, const size_t j, const unsigned short x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
unsigned short *
gsl_matrix_ushort_ptr(gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (unsigned short *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const unsigned short *
gsl_matrix_ushort_const_ptr(const gsl_matrix_ushort * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const unsigned short *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_USHORT_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/* ode-initval2/cscal.c
 *
 * Copyright (C) 2002, 2007 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv2.h>
#include "control_utils.c"

typedef struct
{
  double eps_abs;
  double eps_rel;
  double a_y;
  double a_dydt;
  double *scale_abs;
}
sc_control_state_t;

static void *
sc_control_alloc (void)
{
  sc_control_state_t *s =
    (sc_control_state_t *) malloc (sizeof (sc_control_state_t));

  if (s == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for sc_control_state",
                      GSL_ENOMEM);
    }

  return s;
}

static int
sc_control_init (void *vstate,
                 double eps_abs, double eps_rel, double a_y, double a_dydt)
{
  sc_control_state_t *s = (sc_control_state_t *) vstate;

  if (eps_abs < 0)
    {
      GSL_ERROR ("eps_abs is negative", GSL_EINVAL);
    }
  else if (eps_rel < 0)
    {
      GSL_ERROR ("eps_rel is negative", GSL_EINVAL);
    }
  else if (a_y < 0)
    {
      GSL_ERROR ("a_y is negative", GSL_EINVAL);
    }
  else if (a_dydt < 0)
    {
      GSL_ERROR ("a_dydt is negative", GSL_EINVAL);
    }

  s->eps_rel = eps_rel;
  s->eps_abs = eps_abs;
  s->a_y = a_y;
  s->a_dydt = a_dydt;

  return GSL_SUCCESS;
}

static int
sc_control_hadjust (void *vstate, size_t dim, unsigned int ord,
                    const double y[], const double yerr[], const double yp[],
                    double *h)
{
  sc_control_state_t *state = (sc_control_state_t *) vstate;

  const double eps_abs = state->eps_abs;
  const double eps_rel = state->eps_rel;
  const double a_y = state->a_y;
  const double a_dydt = state->a_dydt;
  const double *scale_abs = state->scale_abs;

  const double S = 0.9;
  const double h_old = *h;

  double rmax = DBL_MIN;
  size_t i;

  for (i = 0; i < dim; i++)
    {
      const double D0 =
        eps_rel * (a_y * fabs (y[i]) + a_dydt * fabs (h_old * yp[i]))
        + eps_abs * scale_abs[i];
      const double r = fabs (yerr[i]) / fabs (D0);
      rmax = GSL_MAX_DBL (r, rmax);
    }

  if (rmax > 1.1)
    {
      /* decrease step, no more than factor of 5, but a fraction S more
         than scaling suggests (for better accuracy) */
      double r = S / pow (rmax, 1.0 / ord);

      if (r < 0.2)
        r = 0.2;

      *h = r * h_old;

      return GSL_ODEIV_HADJ_DEC;
    }
  else if (rmax < 0.5)
    {
      /* increase step, no more than factor of maxscale */
      double r = S / pow (rmax, 1.0 / (ord + 1.0));
      const double maxscale = 4.9;

      if (r > maxscale)
        r = maxscale;

      if (r < 1.0)              /* don't allow any decrease caused by S<1 */
        r = 1.0;

      *h = r * h_old;

      return GSL_ODEIV_HADJ_INC;
    }
  else
    {
      /* no change */
      return GSL_ODEIV_HADJ_NIL;
    }
}

static int
sc_control_errlevel (void *vstate, const double y, const double dydt,
                     const double h, const size_t ind, double *errlev)
{
  sc_control_state_t *state = (sc_control_state_t *) vstate;

  const double eps_abs = state->eps_abs;
  const double eps_rel = state->eps_rel;
  const double a_y = state->a_y;
  const double a_dydt = state->a_dydt;
  const double *scale_abs = state->scale_abs;

  *errlev = eps_rel * (a_y * fabs (y) + a_dydt * fabs (h * dydt))
    + eps_abs * scale_abs[ind];

  if (*errlev <= 0.0)
    {
      GSL_ERROR ("errlev <= zero", GSL_ESANITY);
    }

  return GSL_SUCCESS;
}

static void
sc_control_free (void *vstate)
{
  sc_control_state_t *state = (sc_control_state_t *) vstate;
  free (state->scale_abs);
  free (state);
}

static const gsl_odeiv2_control_type sc_control_type = { "scaled",      /* name */
  &sc_control_alloc,
  &sc_control_init,
  &sc_control_hadjust,
  &sc_control_errlevel,
  &control_set_driver_null,
  &sc_control_free
};

const gsl_odeiv2_control_type *gsl_odeiv2_control_scaled = &sc_control_type;

gsl_odeiv2_control *
gsl_odeiv2_control_scaled_new (double eps_abs, double eps_rel,
                               double a_y, double a_dydt,
                               const double scale_abs[], size_t dim)
{
  gsl_odeiv2_control *c =
    gsl_odeiv2_control_alloc (gsl_odeiv2_control_scaled);

  int status = gsl_odeiv2_control_init (c, eps_abs, eps_rel, a_y, a_dydt);

  if (status != GSL_SUCCESS)
    {
      gsl_odeiv2_control_free (c);
      GSL_ERROR_NULL ("error trying to initialize control", status);
    }

  {
    sc_control_state_t *s = (sc_control_state_t *) c->state;

    s->scale_abs = (double *) malloc (dim * sizeof (double));

    if (s->scale_abs == 0)
      {
        free (s);
        GSL_ERROR_NULL ("failed to allocate space for scale_abs", GSL_ENOMEM);
      }

    memcpy (s->scale_abs, scale_abs, dim * sizeof (double));
  }

  return c;
}

max_stars_repo_name: igormcoelho/optstats

max_issues_repo_name: igormcoelho/optstats

max_forks_repo_name: igormcoelho/optstats


text:
/*
 * Gates.h
 *
 *  Created on: 28 Jan 2019
 *      Author: jobinrjohnson
 */

#ifndef IMPLEMENTATION_GATES_H_
#define IMPLEMENTATION_GATES_H_

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>

#define ccout std::cout
#define cendl std::endl

#include <qubit.h>

#define GATE_HARDAMARD 1
#define GATE_PAULI_X 2
#define GATE_PAULI_Y 6
#define GATE_PAULI_Z 7
#define GATE_CNOT 3
#define GATE_SWAP 4
#define GATE_T 5

#include <iostream>

class Gates
{

  protected:
	gsl_matrix_complex *get_gate_matrix(int);

  public:
	Gates();
	virtual ~Gates();

	static bool get_measure(Qubit &q)
	{

		// TODO Better random and measure implementation

		gsl_vector_complex *state = q.get_state();
		double amp_0 = gsl_complex_abs2(gsl_vector_complex_get(state, 0));

		const gsl_rng_type *T;
		gsl_rng *r;
		gsl_rng_env_setup();
		gsl_rng_default_seed = time(NULL);
		T = gsl_rng_default;
		r = gsl_rng_alloc(T);

		double u = gsl_rng_uniform(r);

		gsl_rng_free(r);

		if (u >= amp_0)
		{
			return true;
		}

		return false;
	}
};

#endif /* IMPLEMENTATION_GATES_H_ */

max_stars_repo_name: jobinrjohnson/QSim

max_issues_repo_name: jobinrjohnson/QSim

max_forks_repo_name: jobinrjohnson/QSim


text:
#include <stdlib.h>

#include <gsl/gsl_statistics_double.h>
#include <hdf5.h>
#include <hdf5_hl.h>

#include "qdm.h"

qdm_ijk *
qdm_ijk_alloc(
    size_t size1,
    size_t size2,
    size_t size3
)
{
  qdm_ijk *t = malloc(sizeof(qdm_ijk));

  t->size1 = size1;
  t->size2 = size2;
  t->size3 = size3;

  t->data = malloc(sizeof(double) * size1 * size2 * size3);
  t->owner = true;

  return t;
}

qdm_ijk *
qdm_ijk_calloc(
    size_t size1,
    size_t size2,
    size_t size3
)
{
  qdm_ijk *t = malloc(sizeof(qdm_ijk));

  t->size1 = size1;
  t->size2 = size2;
  t->size3 = size3;

  t->data = calloc(sizeof(double), size1 * size2 * size3);
  t->owner = true;

  return t;
}

void
qdm_ijk_free(
    qdm_ijk *t
)
{
  if (t == NULL) {
    return;
  }

  if (t->owner) {
    free(t->data);
  }

  free(t);
}

qdm_ijk_view
qdm_ijk_view_array(
    double *base,
    size_t size1,
    size_t size2,
    size_t size3
)
{
  qdm_ijk_view v = {
    .ijk = {
      .size1 = size1,
      .size2 = size2,
      .size3 = size3,

      .data = base,
      .owner = false,
    },
  };

  return v;
}

gsl_vector_view
qdm_ijk_get_k(
    qdm_ijk *t,
    size_t i,
    size_t j
)
{
  return gsl_vector_view_array_with_stride(
      &t->data[i * t->size2 + j],
      t->size1 * t->size2,
      t->size3
  );
}

gsl_matrix_view
qdm_ijk_get_ij(
    qdm_ijk *t,
    size_t k
)
{
  return gsl_matrix_view_array(
      &t->data[k * t->size1 * t->size2],
      t->size1,
      t->size2
  );
}

double
qdm_ijk_get(
    qdm_ijk *t,
    size_t i,
    size_t j,
    size_t k
)
{
  gsl_matrix_view ij = qdm_ijk_get_ij(t, k);

  return gsl_matrix_get(&ij.matrix, i, j);
}

gsl_matrix *
qdm_ijk_cov(
    qdm_ijk *t
)
{
  size_t d = t->size1 * t->size2;
  size_t dr = 0;

  gsl_matrix *r = gsl_matrix_calloc(d, d);

  for (size_t i0 = 0; i0 < t->size1; i0++) {
    for (size_t j0 = 0; j0 < t->size2; j0++) {
      gsl_vector_view k0 = qdm_ijk_get_k(t, i0, j0);

      /* TODO: When we have time (and if it matters), this is symmetric and we
       * don't need to do a bunch of these if we get a little clever about
       * copying the mirror'd results.
       */
      for (size_t ip = 0; ip < t->size1; ip++) {
        for (size_t jp = 0; jp < t->size2; jp++) {
          gsl_vector_view kp = qdm_ijk_get_k(t, ip, jp);

          double cov = gsl_stats_covariance(
              k0.vector.data,
              k0.vector.stride,
              kp.vector.data,
              kp.vector.stride,
              kp.vector.size
          );

          gsl_matrix_set(r, dr, ip * t->size2 + jp, cov);
        }
      }

      dr++;
    }
  }

  return r;
}

int
qdm_ijk_write(
    hid_t id,
    const char *name,
    const qdm_ijk *t
)
{
  int status = 0;

  hid_t datatype  = -1;
  hid_t dataspace = -1;
  hid_t dataset   = -1;

  if (t == NULL) {
    goto cleanup;
  }

  datatype = H5Tcopy(H5T_NATIVE_DOUBLE);

  status = H5Tset_order(datatype, H5T_ORDER_LE);
  if (status != 0) {
    goto cleanup;
  }

  hsize_t dims[3] = {
    t->size3,
    t->size1,
    t->size2,
  };
  dataspace = H5Screate_simple(3, dims, NULL);
  if (dataspace < 0) {
    status = dataspace;
    goto cleanup;
  }

  dataset = H5Dcreate(id, name, datatype, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  if (dataset < 0) {
    status = dataset;
    goto cleanup;
  }

  status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, t->data);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  if (dataset >= 0) {
    H5Dclose(dataset);
  }

  if (dataspace >= 0) {
    H5Sclose(dataspace);
  }

  if (datatype >= 0) {
    H5Tclose(datatype);
  }

  H5Oflush(id);

  return status;
}

int
qdm_ijk_read(
    hid_t id,
    const char *name,
    qdm_ijk **t
)
{
  int status = 0;

  int rank = 0;

  status = H5LTget_dataset_ndims(id, name, &rank);
  if (status < 0) {
    return status;
  }

  if (rank != 3) {
    return -1;
  }

  hsize_t dims[3];

  status = H5LTget_dataset_info(id, name, dims, NULL, NULL);
  if (status < 0) {
    return status;
  }

  size_t size3 = dims[0];
  size_t size1 = dims[1];
  size_t size2 = dims[2];

  qdm_ijk *tmp = qdm_ijk_alloc(size1, size2, size3);

  status = H5LTread_dataset_double(id, name, tmp->data);
  if (status < 0) {
    qdm_ijk_free(tmp);

    return status;
  }

  *t = tmp;

  return status;
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
//**************************************************************************\
//* This file is property of and copyright by the ALICE Project            *\
//* ALICE Experiment at CERN, All rights reserved.                         *\
//*                                                                        *\
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *\
//*                  for The ALICE HLT Project.                            *\
//*                                                                        *\
//* Permission to use, copy, modify and distribute this software and its   *\
//* documentation strictly for non-commercial purposes is hereby granted   *\
//* without fee, provided that the above copyright notice appears in all   *\
//* copies and that both the copyright notice and this permission notice   *\
//* appear in the supporting documentation. The authors make no claims     *\
//* about the suitability of this software for any purpose. It is          *\
//* provided "as is" without express or implied warranty.                  *\
//**************************************************************************

/// \file GPUQA.h
/// \author David Rohr

#ifndef GPUQA_H
#define GPUQA_H

#include "GPUQAConfig.h"
struct AliHLTTPCClusterMCWeight;
class TH1F;
class TH2F;
class TCanvas;
class TPad;
class TLegend;
class TPad;
class TH1;
class TFile;
class TH1D;
typedef short int Color_t;

#if !defined(GPUCA_BUILD_QA) || defined(GPUCA_GPUCODE)

namespace GPUCA_NAMESPACE
{
namespace gpu
{
class GPUChainTracking;

class GPUQA
{
 public:
  GPUQA(GPUChainTracking* rec) {}
  ~GPUQA() = default;

  typedef structConfigQA configQA;

  int InitQA() { return 1; }
  void RunQA(bool matchOnly = false) {}
  int DrawQAHistograms() { return 1; }
  void SetMCTrackRange(int min, int max) {}
  bool SuppressTrack(int iTrack) const { return false; }
  bool SuppressHit(int iHit) const { return false; }
  bool HitAttachStatus(int iHit) const { return false; }
  int GetMCTrackLabel(unsigned int trackId) const { return -1; }
  bool clusterRemovable(int cid, bool prot) const { return false; }
  static bool QAAvailable() { return false; }
};
} // namespace gpu
} // namespace GPUCA_NAMESPACE

#else

#include "GPUTPCDef.h"
#include <cmath>

#ifdef GPUCA_TPC_GEOMETRY_O2
#include <gsl/span>
#endif

namespace o2
{
class MCCompLabel;
}

class AliHLTTPCClusterMCLabel;

namespace GPUCA_NAMESPACE
{
namespace gpu
{
class GPUChainTracking;
class GPUTPCMCInfo;

class GPUQA
{
 public:
  GPUQA(GPUChainTracking* rec);
  ~GPUQA();

  typedef GPUQAConfig configQA;

  int InitQA();
  void RunQA(bool matchOnly = false);
  int DrawQAHistograms();
  void SetMCTrackRange(int min, int max);
  bool SuppressTrack(int iTrack) const;
  bool SuppressHit(int iHit) const;
  bool HitAttachStatus(int iHit) const;
  int GetMCTrackLabel(unsigned int trackId) const;
  bool clusterRemovable(int cid, bool prot) const;
  static bool QAAvailable() { return true; }

 private:
  struct additionalMCParameters {
    float pt, phi, theta, eta, nWeightCls;
  };

  struct additionalClusterParameters {
    int attached, fakeAttached, adjacent, fakeAdjacent;
    float pt;
  };

  void SetAxisSize(TH1F* e);
  void SetLegend(TLegend* l);
  double* CreateLogAxis(int nbins, float xmin, float xmax);
  void ChangePadTitleSize(TPad* p, float size);
  void DrawHisto(TH1* histo, char* filename, char* options);
  void doPerfFigure(float x, float y, float size);
  void GetName(char* fname, int k);
  template <class T>
  T* GetHist(T*& ee, std::vector<TFile*>& tin, int k, int nNewInput);

#ifdef GPUCA_TPC_GEOMETRY_O2
  using mcLabels_t = gsl::span<const o2::MCCompLabel>;
  using mcLabel_t = o2::MCCompLabel;
  using mcLabelI_t = mcLabel_t;
  using mcInfo_t = GPUTPCMCInfo;
  mcLabels_t GetMCLabel(unsigned int i);
  mcLabel_t GetMCLabel(unsigned int i, unsigned int j);
#else
  using mcLabels_t = AliHLTTPCClusterMCLabel;
  using mcLabel_t = AliHLTTPCClusterMCWeight;
  struct mcLabelI_t {
    int getTrackID() const { return AbsLabelID(track); }
    int getEventID() const { return 0; }
    bool isFake() const { return track < 0; }
    bool isValid() const { return track != MC_LABEL_INVALID; }
    void invalidate() { track = MC_LABEL_INVALID; }
    void setFakeFlag(bool v = true) { track = v ? FakeLabelID(track) : AbsLabelID(track); }
    void setNoise() { track = MC_LABEL_INVALID; }
    bool operator==(const mcLabel_t& l);
    bool operator!=(const mcLabel_t& l) { return !(*this == l); }
    mcLabelI_t() = default;
    mcLabelI_t(const mcLabel_t& l);
    int track = MC_LABEL_INVALID;
  };
  using mcInfo_t = GPUTPCMCInfo;
  const mcLabels_t& GetMCLabel(unsigned int i);
  const mcLabel_t& GetMCLabel(unsigned int i, unsigned int j);
  const mcInfo_t& GetMCTrack(const mcLabelI_t& label);
  static int FakeLabelID(const int id);
  static int AbsLabelID(const int id);
#endif
  template <class T>
  static auto& GetMCTrackObj(T& obj, const mcLabelI_t& l);

  unsigned int GetNMCCollissions();
  unsigned int GetNMCTracks(int iCol);
  unsigned int GetNMCLabels();
  const mcInfo_t& GetMCTrack(unsigned int iTrk, unsigned int iCol);
  const mcInfo_t& GetMCTrack(const mcLabel_t& label);
  int GetMCLabelNID(const mcLabels_t& label);
  int GetMCLabelNID(unsigned int i);
  int GetMCLabelID(unsigned int i, unsigned int j);
  int GetMCLabelCol(unsigned int i, unsigned int j);
  static int GetMCLabelID(const mcLabels_t& label, unsigned int j);
  static int GetMCLabelID(const mcLabel_t& label);
  float GetMCLabelWeight(unsigned int i, unsigned int j);
  float GetMCLabelWeight(const mcLabels_t& label, unsigned int j);
  float GetMCLabelWeight(const mcLabel_t& label);
  bool mcPresent();

  static bool MCComp(const mcLabel_t& a, const mcLabel_t& b);

  GPUChainTracking* mTracking;
  const configQA& mConfig;

  //-------------------------: Some compile time settings....
  static const constexpr bool PLOT_ROOT = 0;
  static const constexpr bool FIX_SCALES = 0;
  static const constexpr bool PERF_FIGURE = 0;
  static const constexpr float FIXED_SCALES_MIN[5] = {-0.05, -0.05, -0.2, -0.2, -0.5};
  static const constexpr float FIXED_SCALES_MAX[5] = {0.4, 0.7, 5, 3, 6.5};
  static const constexpr float LOG_PT_MIN = -1.;

  const char* str_perf_figure_1 = "ALICE Performance 2018/03/20";
  // const char* str_perf_figure_2 = "2015, MC pp, #sqrt{s} = 5.02 TeV";
  const char* str_perf_figure_2 = "2015, MC Pb-Pb, #sqrt{s_{NN}} = 5.02 TeV";
  //-------------------------

  std::vector<mcLabelI_t> mTrackMCLabels;
#ifdef GPUCA_TPC_GEOMETRY_O2
  std::vector<std::vector<int>> mTrackMCLabelsReverse;
  std::vector<std::vector<int>> mRecTracks;
  std::vector<std::vector<int>> mFakeTracks;
  std::vector<std::vector<additionalMCParameters>> mMCParam;
  std::vector<std::vector<mcInfo_t>> mMCInfos;
  std::vector<int> mNColTracks;
#else
  std::vector<int> mTrackMCLabelsReverse[1];
  std::vector<int> mRecTracks[1];
  std::vector<int> mFakeTracks[1];
  std::vector<additionalMCParameters> mMCParam[1];
#endif
  std::vector<additionalClusterParameters> mClusterParam;
  int mNTotalFakes = 0;

  TH1F* mEff[4][2][2][5][2]; // eff,clone,fake,all - findable - secondaries - y,z,phi,eta,pt - work,result
  TCanvas* mCEff[6];
  TPad* mPEff[6][4];
  TLegend* mLEff[6];

  TH1F* mRes[5][5][2]; // y,z,phi,lambda,pt,ptlog res - param - res,mean
  TH2F* mRes2[5][5];
  TCanvas* mCRes[7];
  TPad* mPRes[7][5];
  TLegend* mLRes[6];

  TH1F* mPull[5][5][2]; // y,z,phi,lambda,pt,ptlog res - param - res,mean
  TH2F* mPull2[5][5];
  TCanvas* mCPull[7];
  TPad* mPPull[7][5];
  TLegend* mLPull[6];

  static constexpr int N_CLS_HIST = 8;
  static constexpr int N_CLS_TYPE = 3;
  enum CL_types { CL_attached = 0,
                  CL_fake = 1,
                  CL_att_adj = 2,
                  CL_fakeAdj = 3,
                  CL_tracks = 4,
                  CL_physics = 5,
                  CL_prot = 6,
                  CL_all = 7 };
  TH1D* mClusters[N_CLS_TYPE * N_CLS_HIST - 1]; // attached, fakeAttached, attach+adjacent, fakeAdjacent, physics, protected, tracks, all / count, rel, integral
  TCanvas* mCClust[N_CLS_TYPE];
  TPad* mPClust[N_CLS_TYPE];
  TLegend* mLClust[N_CLS_TYPE];

  long long int mNRecClustersRejected = 0, mNRecClustersTube = 0, mNRecClustersTube200 = 0, mNRecClustersLoopers = 0, mNRecClustersLowPt = 0, mNRecClusters200MeV = 0, mNRecClustersPhysics = 0, mNRecClustersProt = 0, mNRecClustersUnattached = 0, mNRecClustersTotal = 0, mNRecClustersHighIncl = 0,
                mNRecClustersAbove400 = 0, mNRecClustersFakeRemove400 = 0, mNRecClustersFullFakeRemove400 = 0, mNRecClustersBelow40 = 0, mNRecClustersFakeProtect40 = 0;
  double mNRecClustersUnaccessible = 0;

  TH1F* mTracks;
  TCanvas* mCTracks;
  TPad* mPTracks;
  TLegend* mLTracks;

  TH1F* mNCl;
  TCanvas* mCNCl;
  TPad* mPNCl;
  TLegend* mLNCl;

  int mNEvents = 0;
  bool mQAInitialized = false;
  std::vector<std::vector<int>> mcEffBuffer;
  std::vector<std::vector<int>> mcLabelBuffer;
  std::vector<std::vector<bool>> mGoodTracks;
  std::vector<std::vector<bool>> mGoodHits;

  static constexpr float Y_MAX = 40;
  static constexpr float Z_MAX = 100;
  static constexpr float PT_MIN = GPUCA_MIN_TRACK_PT_DEFAULT;
  static constexpr float PT_MIN2 = 0.1;
  static constexpr float PT_MIN_PRIM = 0.1;
  static constexpr float PT_MIN_CLUST = GPUCA_MIN_TRACK_PT_DEFAULT;
  static constexpr float PT_MAX = 20;
  static constexpr float ETA_MAX = 1.5;
  static constexpr float ETA_MAX2 = 0.9;

  static constexpr float MIN_WEIGHT_CLS = 40;
  static constexpr float FINDABLE_WEIGHT_CLS = 70;

  static constexpr int MC_LABEL_INVALID = -1e9;

  static constexpr bool CLUST_HIST_INT_SUM = false;

  static constexpr const int COLORCOUNT = 12;
  Color_t* mColorNums;

  static const constexpr char* EFF_TYPES[4] = {"Rec", "Clone", "Fake", "All"};
  static const constexpr char* FINDABLE_NAMES[2] = {"", "Findable"};
  static const constexpr char* PRIM_NAMES[2] = {"Prim", "Sec"};
  static const constexpr char* PARAMETER_NAMES[5] = {"Y", "Z", "#Phi", "#lambda", "Relative #it{p}_{T}"};
  static const constexpr char* PARAMETER_NAMES_NATIVE[5] = {"Y", "Z", "sin(#Phi)", "tan(#lambda)", "q/#it{p}_{T} (curvature)"};
  static const constexpr char* VSPARAMETER_NAMES[6] = {"Y", "Z", "Phi", "Eta", "Pt", "Pt_log"};
  static const constexpr char* EFF_NAMES[3] = {"Efficiency", "Clone Rate", "Fake Rate"};
  static const constexpr char* EFFICIENCY_TITLES[4] = {"Efficiency (Primary Tracks, Findable)", "Efficiency (Secondary Tracks, Findable)", "Efficiency (Primary Tracks)", "Efficiency (Secondary Tracks)"};
  static const constexpr double SCALE[5] = {10., 10., 1000., 1000., 100.};
  static const constexpr double SCALE_NATIVE[5] = {10., 10., 1000., 1000., 1.};
  static const constexpr char* XAXIS_TITLES[5] = {"#it{y}_{mc} (cm)", "#it{z}_{mc} (cm)", "#Phi_{mc} (rad)", "#eta_{mc}", "#it{p}_{Tmc} (GeV/#it{c})"};
  static const constexpr char* AXIS_TITLES[5] = {"#it{y}-#it{y}_{mc} (mm) (Resolution)", "#it{z}-#it{z}_{mc} (mm) (Resolution)", "#phi-#phi_{mc} (mrad) (Resolution)", "#lambda-#lambda_{mc} (mrad) (Resolution)", "(#it{p}_{T} - #it{p}_{Tmc}) / #it{p}_{Tmc} (%) (Resolution)"};
  static const constexpr char* AXIS_TITLES_NATIVE[5] = {"#it{y}-#it{y}_{mc} (mm) (Resolution)", "#it{z}-#it{z}_{mc} (mm) (Resolution)", "sin(#phi)-sin(#phi_{mc}) (Resolution)", "tan(#lambda)-tan(#lambda_{mc}) (Resolution)", "q*(q/#it{p}_{T} - q/#it{p}_{Tmc}) (Resolution)"};
  static const constexpr char* AXIS_TITLES_PULL[5] = {"#it{y}-#it{y}_{mc}/#sigma_{y} (Pull)", "#it{z}-#it{z}_{mc}/#sigma_{z} (Pull)", "sin(#phi)-sin(#phi_{mc})/#sigma_{sin(#phi)} (Pull)", "tan(#lambda)-tan(#lambda_{mc})/#sigma_{tan(#lambda)} (Pull)",
                                                      "q*(q/#it{p}_{T} - q/#it{p}_{Tmc})/#sigma_{q/#it{p}_{T}} (Pull)"};
  static const constexpr char* CLUSTER_NAMES[N_CLS_HIST] = {"Correctly attached clusters", "Fake attached clusters", "Attached + adjacent clusters", "Fake adjacent clusters", "Clusters of reconstructed tracks", "Used in Physics", "Protected", "All clusters"};
  static const constexpr char* CLUSTER_TITLES[N_CLS_TYPE] = {"Clusters Pt Distribution / Attachment", "Clusters Pt Distribution / Attachment (relative to all clusters)", "Clusters Pt Distribution / Attachment (integrated)"};
  static const constexpr char* CLUSTER_NAMES_SHORT[N_CLS_HIST] = {"Attached", "Fake", "AttachAdjacent", "FakeAdjacent", "FoundTracks", "Physics", "Protected", "All"};
  static const constexpr char* CLUSTER_TYPES[N_CLS_TYPE] = {"", "Ratio", "Integral"};
  static const constexpr int COLORS_HEX[COLORCOUNT] = {0xB03030, 0x00A000, 0x0000C0, 0x9400D3, 0x19BBBF, 0xF25900, 0x7F7F7F, 0xFFD700, 0x07F707, 0x07F7F7, 0xF08080, 0x000000};

  static const constexpr int CONFIG_DASHED_MARKERS = 0;

  static const constexpr float AXES_MIN[5] = {-Y_MAX, -Z_MAX, 0.f, -ETA_MAX, PT_MIN};
  static const constexpr float AXES_MAX[5] = {Y_MAX, Z_MAX, 2.f * M_PI, ETA_MAX, PT_MAX};
  static const constexpr int AXIS_BINS[5] = {51, 51, 144, 31, 50};
  static const constexpr int RES_AXIS_BINS[] = {1017, 113}; // Consecutive bin sizes, histograms are binned down until the maximum entry is 50, each bin size should evenly divide its predecessor.
  static const constexpr float RES_AXES[5] = {1., 1., 0.03, 0.03, 1.0};
  static const constexpr float RES_AXES_NATIVE[5] = {1., 1., 0.1, 0.1, 5.0};
  static const constexpr float PULL_AXIS = 10.f;

  int mMCTrackMin = -1, mMCTrackMax = -1;
};

inline bool GPUQA::SuppressTrack(int iTrack) const { return (mConfig.matchMCLabels.size() && !mGoodTracks[mNEvents][iTrack]); }
inline bool GPUQA::SuppressHit(int iHit) const { return (mConfig.matchMCLabels.size() && !mGoodHits[mNEvents - 1][iHit]); }
inline bool GPUQA::HitAttachStatus(int iHit) const { return (mClusterParam.size() && mClusterParam[iHit].fakeAttached ? (mClusterParam[iHit].attached ? 1 : 2) : 0); }

} // namespace gpu
} // namespace GPUCA_NAMESPACE

#endif
#endif

max_stars_repo_name: mpoghos/AliRoot

max_issues_repo_name: mpoghos/AliRoot

max_forks_repo_name: mpoghos/AliRoot


text:
#pragma once

#include <cstdint>
#include <array>
#include <gsl/gsl>

#include "halley/data_structures/maybe.h"
#include "halley/maths/vector2.h"
#include "halley/maths/colour.h"
#include "halley/file_formats/image.h"

namespace Halley {
	enum class AsepriteBlendMode;

	enum class AsepriteDepth
	{
		RGBA32,
		Greyscale16,
		Indexed8
	};

	struct AsepriteCel
	{
	public:
		Vector2i pos;
		Vector2i size;
		uint16_t layer = 0;
		uint16_t linkedFrame = 0;
		uint8_t opacity = 255;
		bool linked = false;

		Bytes rawData;
		std::unique_ptr<Image> imgData;

		void loadImage(AsepriteDepth depth, const Vector<uint32_t>& palette);
		void drawAt(Image& image, uint8_t opacity, AsepriteBlendMode blendMode) const;
	};

	struct AsepriteFrame
	{
	public:
		AsepriteFrame(uint16_t duration);

		uint16_t duration;
		Vector<AsepriteCel> cels;
	};

	enum class AsepriteLayerType
	{
		Normal,
		Group
	};

	enum class AsepriteBlendMode
	{
		Normal,
		Multiply,
		Screen,
		Overlay,
		Darken,
		Lighten,
		ColorDodge,
		ColorBurn,
		HardLight,
		SoftLight,
		Difference,
		Exclusion,
		Hue,
		Saturation,
		Color,
		Luminosity,
		Addition,
		Subtract,
		Divide
	};

	struct AsepriteLayer
	{
		AsepriteLayerType type = AsepriteLayerType::Normal;
		AsepriteBlendMode blendMode = AsepriteBlendMode::Normal;
		int childLevel = 0;
		bool visible = true;
		bool editable = true;
		bool lockMovement = false;
		bool background = false;
		bool preferLinkedCels = false;
		bool layerGroupDisplaysCollapsed = false;
		bool referenceLayer = false;
		uint8_t opacity = 255;
		String layerName;

		OptionalLite<uint32_t> parentIdx;
		bool visibleInHierarchy = true;
	};

	enum class AsepriteAnimationDirection
	{
		Forward,
		Reverse,
		PingPong
	};

	struct AsepriteTag
	{
		int fromFrame = -1;
		int toFrame = -1;
		AsepriteAnimationDirection animDirection = AsepriteAnimationDirection::Forward;

		String name;
	};

    class AsepriteFile {
    public:
		AsepriteFile();

		void load(gsl::span<const gsl::byte> data);

		const Vector<AsepriteTag>& getTags() const;
		std::map<String, std::unique_ptr<Image>> makeGroupFrameImages(int frameNumber, bool groupSeparated);
	    const AsepriteFrame& getFrame(int n) const;
	    size_t getNumberOfFrames() const;

    private:
	    void addFrame(uint16_t duration);
	    void addChunk(uint16_t chunkType, gsl::span<const gsl::byte> data);

	    void addLayerChunk(gsl::span<const gsl::byte> span);
		void addCelChunk(gsl::span<const gsl::byte> span);
		void addCelExtraChunk(gsl::span<const gsl::byte> span);
		void addPaletteChunk(gsl::span<const gsl::byte> span);
		void addTagsChunk(gsl::span<const gsl::byte> span);
    	void postProcessLayers();

		template <typename T>
		void readData(T& dst, gsl::span<const gsl::byte>& data) const
		{
			if (data.size() < T::size) {
				throw Exception("Insufficient data to decode Aseprite entry", HalleyExceptions::Tools);
			}
			memcpy(&dst, data.data(), T::size);
			data = data.subspan(T::size);
		}

		String readString(gsl::span<const gsl::byte>& data) const
		{
			if (data.size() < 2) {
				throw Exception("Insufficient data to decode Aseprite entry", HalleyExceptions::Tools);
			}
			uint16_t size;
			memcpy(&size, data.data(), 2);
			data = data.subspan(2);
			if (data.size() < size) {
				throw Exception("Insufficient data to decode Aseprite entry", HalleyExceptions::Tools);
			}
			String result = String(reinterpret_cast<const char*>(data.data()), size);
			data = data.subspan(size);
			return result;
		}

	    AsepriteCel* getCelAt(int frameNumber, int layerNumber);
	    size_t getBPP() const;

		Vector2i size;
		uint32_t flags;
		AsepriteDepth colourDepth;
		uint8_t transparentEntry;
		uint16_t numOfColours;

		Vector<AsepriteFrame> frames;
		Vector<AsepriteLayer> layers;
		Vector<AsepriteTag> tags;
		Vector<uint32_t> paletteBg;
		Vector<uint32_t> paletteTransparent;
    };
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley