text: /*
   * Copyright (c) 2011-2013, 2018-2019 Genome Research Ltd.
   * Author(s): James Bonfield
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions are met:
   *
   *    1. Redistributions of source code must retain the above copyright notice,
   *       this list of conditions and the following disclaimer.
   *
   *    2. Redistributions in binary form must reproduce the above
   *       copyright notice, this list of conditions and the following
   *       disclaimer in the documentation and/or other materials provided
   *       with the distribution.
   *
   *    3. Neither the names Genome Research Ltd and Wellcome Trust Sanger
   *       Institute nor the names of its contributors may be used to endorse
   *       or promote products derived from this software without specific
   *       prior written permission.
   *
   * THIS SOFTWARE IS PROVIDED BY GENOME RESEARCH LTD AND CONTRIBUTORS "AS
   * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
   * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
   * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL GENOME RESEARCH
   * LTD OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
   * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
   * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
   * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
   * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
   */
  
  // We use generic maps to turn 0-M into 0-N where N <= M
  // before adding these into the context.  These are used
  // for positions, running-diffs and quality values.
  //
  // This can be used as a simple divisor, eg pos/24 to get
  // 2 bits of positional data for each quarter along a 100bp
  // read, or it can be tailored for specific such as noting
  // the first 5 cycles are poor, then we have stability and
  // a gradual drop off in the last 20 or so.  Perhaps we then
  // map pos 0-4=0, 5-79=1, 80-89=2, 90-99=3.
  //
  // We don't need to specify how many bits of data we are
  // using (2 in the above example), as that is just implicit
  // in the values in the map.  Specify not to use a map simply
  // disables that context type (our map is essentially 0-M -> 0).
  
  // Example of command line usage:
  //
  // f=~/scratch/data/q4
  // cc -Wall -DTEST_MAIN -O3 -g fqzcomp_qual2.c -lm
  // ./a.out $f > /tmp/_ && ./a.out -d < /tmp/_ > /tmp/__ && cmp /tmp/__ $f
  
  #include "config.h"
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <assert.h>
  #include <string.h>
  #include <math.h>
  #include <limits.h>
  #include <ctype.h>
  #include <math.h>
  #include <inttypes.h>
  #include <sys/types.h>
  
  //#define NO_THREADS
  #ifndef NO_THREADS
  #include <pthread.h>
  #endif
  
  #include "fqzcomp_qual.h"
  #include "varint.h"
  
  #define CTX_BITS 16
  #define CTX_SIZE (1<<CTX_BITS)
  
  #ifndef MIN
  #  define MIN(a,b) ((a)<(b)?(a):(b))
  #  define MAX(a,b) ((a)>(b)?(a):(b))
  #endif
  
  #define QMAX 256
  #define QBITS 12
  #define QSIZE (1<<QBITS)
  
  #define NSYM 2
  #include "c_simple_model.h"
  
  #undef NSYM
  #define NSYM QMAX
  //#include "c_escape_model.h"
  #include "c_simple_model.h"
  //#include "c_cdf_model.h"
  //#include "c_cdf16_model.h"
  
  // An array of 0,0,0, 1,1,1,1, 3, 5,5
  // is turned into a run-length of 3x0, 4x1, 0x2, 1x4, 0x4, 2x5,
  // which then becomes 3 4 0 1 0 2.
  //
  // NB: size > 255 therefore means we need to repeatedly read to find
  // the actual run length.
  // Alternatively we could bit-encode instead of byte encode, eg BETA.
  static int store_array(unsigned char *out, unsigned int *array, int size) {
      unsigned char tmp[2048];
  
      int i, j, k;
      for (i = j = k = 0; i < size; j++) {
  	int run_len = i;
  	while (i < size && array[i] == j)
  	    i++;
  	run_len = i-run_len;
  
  	int r;
  	do {
  	    r = MIN(255, run_len);
  	    tmp[k++] = r;
  	    run_len -= r;
  	} while (r == 255);
      }
      while (i < size)
  	tmp[k++] = 0, j++;
  
      // RLE on out.
      //    1 2 3 3 3 3 3 4 4    5
      // => 1 2 3 3 +3... 4 4 +0 5
      int last = -1;
      for (i = j = 0; j < k; i++) {
  	out[i] = tmp[j++];
  	if (out[i] == last) {
  	    int n = j;
  	    while (j < k && tmp[j] == last)
  		j++;
  	    out[++i] = j-n;
  	} else {
  	    last = out[i];
  	}
      }
      k = i;
  
  //    fprintf(stderr, "Store_array %d => %d {", size, k);
  //    for (i = 0; i < k; i++)
  //	fprintf(stderr, "%d,", out[i]);
  //    fprintf(stderr, "}\n");
      return k;
  }
  
  static int read_array(unsigned char *in, size_t in_size, unsigned int *array, int size) {
      unsigned char R[1024];
      int i, j, z, last = -1, nb = 0;
  
      size = MIN(1024, size);
  
      // Remove level one of run-len encoding
      for (i = j = z = 0; z < size && i < in_size; i++) {
  	int run = in[i];
  	R[j++] = run;
  	z += run;
  	if (run == last) {
  	    if (i+1 >= in_size)
  		return -1;
  	    int copy = in[++i];
  	    z += run * copy;
  	    while (copy-- && z < size && j < 1024)
  		R[j++] = run;
  	}
  	if (j >= 1024)
  	    return -1;
  	last = run;
      }
      nb = i;
  
      // Now expand inner level of run-length encoding
      int R_max = j;
      for (i = j = z = 0; j < size; i++) {
  	int run_len = 0;
  	int run_part;
  	if (z >= R_max)
  	    return -1;
  	do {
  	    run_part = R[z++];
  	    run_len += run_part;
  	} while (run_part == 255 && z < R_max);
  	if (run_part == 255)
  	    return -1;
  
  	while (run_len && j < size)
  	    run_len--, array[j++] = i;
      }
  
      return nb;
  }
  
  // FIXME: how to auto-tune these rather than trial and error?
  // r2 = READ2
  // qa = qual avg (0, 2, 4)
  static int strat_opts[][12] = {
  //   qb  qs pb ps db ds ql sl pl  dl  r2 qa
      {10, 5, 4,-1, 2, 1, 0, 14, 10, 14, 0,-1}, // basic options (level < 7)
      {8,  5, 7, 0, 0, 0, 0, 14, 8,  14, 1,-1}, // e.g. HiSeq 2000
      {12, 6, 2, 0, 2, 3, 0, 9,  12, 14, 0, 0}, // e.g. MiSeq
      {12, 6, 0, 0, 0, 0, 0, 12, 0,  0,  0, 0}, // e.g. IonTorrent; adaptive O1
      {0,  0, 0, 0, 0, 0, 0, 0,  0,  0,  0, 0}, // custom
  };
  static int nstrats = sizeof(strat_opts) / sizeof(*strat_opts);
  
  #ifdef __SSE__
  #   include <xmmintrin.h>
  #else
  #   define _mm_prefetch(a,b)
  #endif
  
  typedef struct {
      unsigned int qctx;  // quality sub-context
      unsigned int p;     // pos (bytes remaining)
      unsigned int add_d; // whether to update delta (skip first cycle)
      unsigned int delta; // delta running total
      unsigned int prevq; // previous quality
      unsigned int s;     // selector
      unsigned int qtot, qlen;
      unsigned int first_len;
  } fqz_state;
  
  static void dump_table(unsigned int *tab, int size, char *name) {
      int i, last = -99, run = 0;
      fprintf(stderr, "\t%s\t{", name);
      for (i = 0; i < size; i++) {
  	if (tab[i] == last) {
  	    run++;
  	} else if (run == 1 && tab[i] == last+1) {
  	    int first = last;
  	    do {
  		last = tab[i];
  		i++;
  	    } while (i < size && tab[i] == last+1);
  	    i--;
  
  	    // Want 0,1,2,3,3,3 as 0..2 3x3, not 0..3 3x2
  	    if (tab[i] == tab[i+1])
  		i--;
  	    if (tab[i] != first)
  		fprintf(stderr, "..%d", tab[i]);
  	    run = 1;
  	    last = -99;
  	} else {
  	    if (run > 1)
  		fprintf(stderr, " x %d%s%d", run, i?", ":"", tab[i]);
  	    else
  		fprintf(stderr, "%s%d", i?", ":"", tab[i]);
  	    run = 1;
  	    last = tab[i];
  	}
      }
      if (run > 1)
  	fprintf(stderr, " x %d", run);
      fprintf(stderr, "}\n");
  }
  
  static void dump_map(unsigned int *map, int size, char *name) {
      int i, c = 0;
      fprintf(stderr, "\t%s\t{", name);
      for (i = 0; i < size; i++)
  	if (map[i] != INT_MAX)
  	    fprintf(stderr, "%s%d=%d", c++?", ":"", i, map[i]);
      fprintf(stderr, "}\n");
  }
  
  #pragma GCC diagnostic ignored "-Wunused-function"
  static void dump_params(fqz_gparams *gp) {
      fprintf(stderr, "Global params = {\n");
      fprintf(stderr, "\tvers\t%d\n", gp->vers);
      fprintf(stderr, "\tgflags\t0x%02x\n", gp->gflags);
      fprintf(stderr, "\tnparam\t%d\n", gp->nparam);
      fprintf(stderr, "\tmax_sel\t%d\n", gp->max_sel);
      fprintf(stderr, "\tmax_sym\t%d\n", gp->max_sym);
      if (gp->gflags & GFLAG_HAVE_STAB)
  	dump_table(gp->stab, 256, "stab");
      fprintf(stderr, "}\n");
  
      int i;
      for (i = 0; i < gp->nparam; i++) {
  	fqz_param *pm = &gp->p[i];
  	fprintf(stderr, "\nParam[%d] = {\n", i);
  	fprintf(stderr, "\tcontext\t0x%04x\n", pm->context);
  	fprintf(stderr, "\tpflags\t0x%02x\n",  pm->pflags);
  	fprintf(stderr, "\tmax_sym\t%d\n",  pm->max_sym);
  	fprintf(stderr, "\tqbits\t%d\n",   pm->qbits);
  	fprintf(stderr, "\tqshift\t%d\n",  pm->qshift);
  	fprintf(stderr, "\tqloc\t%d\n",    pm->qloc);
  	fprintf(stderr, "\tsloc\t%d\n",    pm->sloc);
  	fprintf(stderr, "\tploc\t%d\n",    pm->ploc);
  	fprintf(stderr, "\tdloc\t%d\n",    pm->dloc);
  
  	if (pm->pflags & PFLAG_HAVE_QMAP)
  	    dump_map(pm->qmap, 256, "qmap");
  
  	if (pm->pflags & PFLAG_HAVE_QTAB)
  	    dump_table(pm->qtab, 256, "qtab");
  	if (pm->pflags & PFLAG_HAVE_PTAB)
  	    dump_table(pm->ptab, 1024, "ptab");
  	if (pm->pflags & PFLAG_HAVE_DTAB)
  	    dump_table(pm->dtab, 256, "dtab");
  	fprintf(stderr, "}\n");
      }
  }
  
  typedef struct {
      SIMPLE_MODEL(QMAX,_) *qual;
      SIMPLE_MODEL(256,_)   len[4];
      SIMPLE_MODEL(2,_)     revcomp;
      SIMPLE_MODEL(256,_)   sel;
      SIMPLE_MODEL(2,_)     dup;
  } fqz_model;
  
  #ifndef NO_THREADS
  /*
   * Thread local storage, used to avoid repeated malloc/free calls.
   */
  pthread_once_t fqz_once = PTHREAD_ONCE_INIT;
  pthread_key_t fqz_key;
  
  static void fqz_tls_init(void) {
      pthread_key_create(&fqz_key, free);
  }
  #endif
  
  static int fqz_create_models(fqz_model *m, fqz_gparams *gp) {
      int i;
  
  #ifndef NO_THREADS
      pthread_once(&fqz_once, fqz_tls_init);
  
      m->qual = pthread_getspecific(fqz_key);
  
      if (!m->qual) {
          if (!(m->qual = malloc(sizeof(*m->qual) * CTX_SIZE)))
  	    return -1;
  	pthread_setspecific(fqz_key, m->qual);
      }
  #else
      if (!(m->qual = malloc(sizeof(*m->qual) * CTX_SIZE)))
  	return -1;
  #endif
  
      for (i = 0; i < CTX_SIZE; i++)
  	SIMPLE_MODEL(QMAX,_init)(&m->qual[i], gp->max_sym+1);
  
      for (i = 0; i < 4; i++)
  	SIMPLE_MODEL(256,_init)(&m->len[i],256);
  
      SIMPLE_MODEL(2,_init)(&m->revcomp,2);
      SIMPLE_MODEL(2,_init)(&m->dup,2);
      if (gp->max_sel > 0)
  	SIMPLE_MODEL(256,_init)(&m->sel, gp->max_sel+1);
  
      return 0;
  }
  
  static void fqz_destroy_models(fqz_model *m) {
  #ifdef NO_THREADS
      free(m->qual);
  #endif
  }
  
  static inline unsigned int fqz_update_ctx(fqz_param *pm, fqz_state *state, int q) {
      unsigned int last = 0; // pm->context
      state->qctx = (state->qctx << pm->qshift) + pm->qtab[q];
      last += (state->qctx & pm->qmask) << pm->qloc;
  
      // The final shifts have been factored into the tables already.
      last += pm->ptab[MIN(1023, state->p)];      // << pm->ploc
      last += pm->dtab[MIN(255,  state->delta)];  // << pm->dloc
      last += state->s << pm->sloc;
  
      // On the fly average is slow work.
      // However it can be slightly better than using a selector bit
      // as it's something we can compute on the fly and thus doesn't
      // consume output bits for storing the selector itself.
      //
      // Q4 (novaseq.bam)
      // qtot+=q*q -DQ1=8.84 -DQ2=8.51 -DQ3=7.70; 7203598 (-0.7%)
      // qtot+=q   -DQ1=2.96 -DQ2=2.85 -DQ3=2.69; 7207315
      // vs old delta;                            7255614 (default params)
      // vs 2 bit selector (no delta)             7203006 (-x 0x8261000e80)
      // vs 2 bit selector (no delta)             7199153 (-x 0x7270000e70) -0.8%
      // vs 2 bit selector (no delta)             7219668 (-x 0xa243000ea0)
      //{
      //	double qa = state->qtot / (state->qlen+.01);
      //	//fprintf(stderr, "%f\n", qa);
      //	int x = 0;
      //	if (qa>=Q1) x=3;
      //	else if (qa>=Q2) x=2;
      //	else if (qa>=Q3) x=1;
      //	else x=0;
      //	last += x << pm->dloc; // tmp reuse of delta pos
      //  state->qtot += q*q;
      //  state->qlen++;
      //}
  
      // Only update delta after 1st base.
      //state->delta += state->add_d * (state->prevq != q);
      //state->add_d = 1;
      state->delta += (state->prevq != q);
      state->prevq = q;
  
      state->p--;
  
      return last & (CTX_SIZE-1);
  }
  
  
  // Build quality stats for qhist and set nsym, do_dedup and do_sel params.
  // One_param is -1 to gather stats on all data, or >= 0 to gather data
  // on one specific selector parameter.  Used only in TEST_MAIN via
  // fqz_manual_parameters at the moment.
  void fqz_qual_stats(fqz_slice *s,
  		    unsigned char *in, size_t in_size,
  		    fqz_param *pm,
  		    uint32_t qhist[256],
  		    int one_param) {
  #define NP 128
      uint32_t qhistb[NP][256] = {{0}};  // both
      uint32_t qhist1[NP][256] = {{0}};  // READ1 only
      uint32_t qhist2[NP][256] = {{0}};  // READ2 only
      uint64_t t1[NP] = {0};             // Count for READ1
      uint64_t t2[NP] = {0};             // COUNT for READ2
      uint32_t avg[2560] = {0};          // Avg qual *and later* avg-to-selector map.
  
      int dir = 0;
      int last_len = 0;
      int do_dedup = 0;
      size_t rec;
      size_t i, j;
      int num_rec = 0;
  
      // See what info we've been given.
      // Do we have READ1 / READ2?
      // Do we have selector hidden in the top bits of flag?
      int max_sel = 0;
      int has_r2 = 0;
      for (rec = 0; rec < s->num_records; rec++) {
  	if (one_param >= 0 && (s->flags[rec] >> 16) != one_param)
  	    continue;
  	num_rec++;
  	if (max_sel < (s->flags[rec] >> 16))
  	    max_sel = (s->flags[rec] >> 16);
  	if (s->flags[rec] & FQZ_FREAD2)
  	    has_r2 = 1;
      }
  
      // Dedup detection and histogram stats gathering
      int *avg_qual = calloc((s->num_records+1), sizeof(int));
      if (!avg_qual)
  	return;
  
      rec = i = j = 0;
      while (i < in_size) {
  	if (one_param >= 0 && (s->flags[rec] >> 16) != one_param) {
  	    avg_qual[rec] = 0;
  	    i += s->len[rec++];
  	    continue;
  	}
  	if (rec < s->num_records) {
  	    j = s->len[rec];
  	    dir = s->flags[rec] & FQZ_FREAD2 ? 1 : 0;
  	    if (i > 0 && j == last_len
  		&& !memcmp(in+i-last_len, in+i, j))
  		do_dedup++; // cache which records are dup?
  	} else {
  	    j = in_size - i;
  	    dir = 0;
  	}
  	last_len = j;
  
  	uint32_t (*qh)[256] = dir ? qhist2 : qhist1;
  	uint64_t *th        = dir ? t2     : t1;
  
  	uint32_t tot = 0;
  	for (; i < in_size && j > 0; i++, j--) {
  	    tot += in[i];
  	    qhist[in[i]]++;
  	    qhistb[j & (NP-1)][in[i]]++;
  	    qh[j & (NP-1)][in[i]]++;
  	    th[j & (NP-1)]++;
  	}
  	tot = last_len ? (tot*10.0)/last_len+.5 : 0;
  
  	avg_qual[rec] = tot;
  	avg[MIN(2559, tot)]++;
  
  	rec++;
      }
      pm->do_dedup = ((rec+1)/(do_dedup+1) < 500);
  
      last_len = 0;
  
      // Unique symbol count
      for (i = pm->max_sym = pm->nsym = 0; i < 256; i++) {
  	if (qhist[i])
  	    pm->max_sym = i, pm->nsym++;
      }
  
  
      // Auto tune: does average quality helps us?
      if (pm->do_qa != 0) {
  	// Histogram of average qual in avg[]
  	// NB: we convert avg[] from count to selector index
  
  	// Few symbols means high compression which means
  	// selector bits become more significant fraction.
  	// Reduce selector bits by skewing the distribution
  	// to not be even binning.
  	double qf0 = pm->nsym > 8 ? 0.2 : 0.05;
  	double qf1 = pm->nsym > 8 ? 0.5 : 0.22;
  	double qf2 = pm->nsym > 8 ? 0.8 : 0.60;
  
  	int total = 0;
  	i = 0;
  	while (i < 2560) {
  	    total += avg[i];
  	    if (total > qf0 * num_rec) {
  		//fprintf(stderr, "Q1=%d\n", (int)i);
  		break;
  	    }
  	    avg[i++] = 0;
  	}
  	while (i < 2560) {
  	    total += avg[i];
  	    if (total > qf1 * num_rec) {
  		//fprintf(stderr, "Q2=%d\n", (int)i);
  		break;
  	    }
  	    avg[i++] = 1;
  	}
  	while (i < 2560) {
  	    total += avg[i];
  	    if (total > qf2 * num_rec) {
  		//fprintf(stderr, "Q3=%d\n", (int)i);
  		break;
  	    }
  	    avg[i++] = 2;
  	}
  	while (i < 2560)
  	    avg[i++] = 3;
  
  	// Compute simple entropy of merged signal vs split signal.
          i = 0;
  	rec = 0;
  
  	int qbin4[4][NP][256] = {{{0}}};
  	int qbin2[2][NP][256] = {{{0}}};
  	int qbin1   [NP][256] = {{0}};
  	int qcnt4[4][NP] = {{0}};
  	int qcnt2[4][NP] = {{0}};
  	int qcnt1   [NP] = {0};
          while (i < in_size) {
  	    if (one_param >= 0 && (s->flags[rec] >> 16) != one_param) {
  		i += s->len[rec++];
  		continue;
  	    }
  	    if (rec < s->num_records)
  		j = s->len[rec];
  	    else
  		j = in_size - i;
  	    last_len = j;
  
  	    uint32_t tot = avg_qual[rec];
  	    int qb4 = avg[MIN(2559, tot)];
  	    int qb2 = qb4/2;
  
  	    for (; i < in_size && j > 0; i++, j--) {
  		int x = j & (NP-1);
  		qbin4[qb4][x][in[i]]++;  qcnt4[qb4][x]++;
  		qbin2[qb2][x][in[i]]++;  qcnt2[qb2][x]++;
  		qbin1     [x][in[i]]++;  qcnt1     [x]++;
  	    }
  	    rec++;
  	}
  
  	double e1 = 0, e2 = 0, e4 = 0;
  	for (j = 0; j < NP; j++) {
  	    for (i = 0; i < 256; i++) {
  		if (qbin1   [j][i]) e1 += qbin1   [j][i] * log(qbin1   [j][i] / (double)qcnt1   [j]);
  		if (qbin2[0][j][i]) e2 += qbin2[0][j][i] * log(qbin2[0][j][i] / (double)qcnt2[0][j]);
  		if (qbin2[1][j][i]) e2 += qbin2[1][j][i] * log(qbin2[1][j][i] / (double)qcnt2[1][j]);
  		if (qbin4[0][j][i]) e4 += qbin4[0][j][i] * log(qbin4[0][j][i] / (double)qcnt4[0][j]);
  		if (qbin4[1][j][i]) e4 += qbin4[1][j][i] * log(qbin4[1][j][i] / (double)qcnt4[1][j]);
  		if (qbin4[2][j][i]) e4 += qbin4[2][j][i] * log(qbin4[2][j][i] / (double)qcnt4[2][j]);
  		if (qbin4[3][j][i]) e4 += qbin4[3][j][i] * log(qbin4[3][j][i] / (double)qcnt4[3][j]);
  	    }
  	}
  	e1 /= -log(2)/8;
  	e2 /= -log(2)/8;
  	e4 /= -log(2)/8;
  	//fprintf(stderr, "E1=%f E2=%f E4=%f\n", e1, e2+s->num_records/8, e4+s->num_records/4);
  
  	// Note by using the selector we're robbing bits from elsewhere in
  	// the context, which may reduce compression better.
  	// We don't know how much by, so this is basically a guess!
  	// For now we just say need 5% saving here.
  	double qm = pm->do_qa > 0 ? 1 : 0.98;
  	if ((pm->do_qa == -1 || pm->do_qa >= 4) &&
  	    e4 + s->num_records/4 < e2*qm + s->num_records/8 &&
  	    e4 + s->num_records/4 < e1*qm) {
  	    //fprintf(stderr, "do q4\n");
  	    for (i = 0; i < s->num_records; i++) {
  		//fprintf(stderr, "%d -> %d -> %d, %d\n", (int)i, avg_qual[i], avg[MIN(2559, avg_qual[i])], s->flags[i]>>16);
  		s->flags[i] |= avg[MIN(2559, avg_qual[i])] <<16;
  	    }
  	    pm->do_sel = 1;
  	    max_sel = 3;
  	} else if ((pm->do_qa == -1 || pm->do_qa >= 2) && e2 + s->num_records/8 < e1*qm) {
  	    //fprintf(stderr, "do q2\n");
  	    for (i = 0; i < s->num_records; i++)
  		s->flags[i] |= (avg[MIN(2559, avg_qual[i])]>>1) <<16;
  	    pm->do_sel = 1;
  	    max_sel = 1;
  	}
  
  	if (pm->do_qa == -1) {
  	    // assume qual, pos, delta in that order.
  	    if (pm->pbits > 0 && pm->dbits > 0) {
  		// 1 from pos/delta
  		pm->sloc = pm->dloc-1;
  		pm->pbits--;
  		pm->dbits--;
  		pm->dloc++;
  	    } else if (pm->dbits >= 2) {
  		// 2 from delta
  		pm->sloc = pm->dloc;
  		pm->dbits -= 2;
  		pm->dloc += 2;
  	    } else if (pm->qbits >= 2) {
  		pm->qbits -= 2;
  		pm->ploc -= 2;
  		pm->sloc = 16-2 - pm->do_r2;
  		if (pm->qbits == 6 && pm->qshift == 5)
  		    pm->qbits--;
  	    }
  	    pm->do_qa = 4;
  	}
      }
  
      // Auto tune: does splitting up READ1 and READ2 help us?
      if (has_r2 || pm->do_r2) { // FIXME: && but debug for now
  	double e1 = 0, e2 = 0; // entropy sum
  
  	for (j = 0; j < NP; j++) {
  	    if (!t1[j] || !t2[j]) continue;
  	    for (i = 0; i < 256; i++) {
  		if (!qhistb[j][i]) continue;
  		e1 -= (qhistb[j][i])*log(qhistb[j][i] / (double)(t1[j]+t2[j]));
  		if (qhist1[j][i])
  		    e2 -= qhist1[j][i] * log(qhist1[j][i] / (double)t1[j]);
  		if (qhist2[j][i])
  		    e2 -= qhist2[j][i] * log(qhist2[j][i] / (double)t2[j]);
  	    }
  	}
  	e1 /= log(2)*8; // bytes
  	e2 /= log(2)*8;
  
  	//fprintf(stderr, "read1/2 entropy merge %f split %f\n", e1, e2);
  
  	// Note by using the selector we're robbing bits from elsewhere in
  	// the context, which may reduce compression better.
  	// We don't know how much by, so this is basically a guess!
  	// For now we just say need 5% saving here.
  	double qm = pm->do_r2 > 0 ? 1 : 0.95;
  	if (e2 + (8+s->num_records/8) < e1*qm) {
  	    for (rec = 0; rec < s->num_records; rec++) {
  		if (one_param >= 0 && (s->flags[rec] >> 16) != one_param)
  		    continue;
  		int sel = s->flags[rec] >> 16;
  		s->flags[rec] =  (s->flags[rec] & 0xffff)
  		    | ((s->flags[rec] & FQZ_FREAD2)
  		       ? ((sel*2)+1) << 16
  		       : ((sel*2)+0) << 16);
  		if (max_sel < (s->flags[rec]>>16))
  		    max_sel = (s->flags[rec]>>16);
  	    }
  	}
      }
  
      // We provided explicit selector data or auto-tuned it
      if (max_sel > 0) {
  	pm->do_sel = 1;
  	pm->max_sel = max_sel;
      }
  
      free(avg_qual);
  }
  
  static inline
  int fqz_store_parameters1(fqz_param *pm, unsigned char *comp) {
      int comp_idx = 0, i, j;
  
      // Starting context
      comp[comp_idx++] = pm->context;
      comp[comp_idx++] = pm->context >> 8;
  
      comp[comp_idx++] = pm->pflags;
      comp[comp_idx++] = pm->max_sym;
  
      comp[comp_idx++] = (pm->qbits<<4)|pm->qshift;
      comp[comp_idx++] = (pm->qloc<<4)|pm->sloc;
      comp[comp_idx++] = (pm->ploc<<4)|pm->dloc;
  
      if (pm->store_qmap) {
  	for (i = j = 0; i < 256; i++)
  	    if (pm->qmap[i] != INT_MAX)
  		comp[comp_idx++] = i;
      }
  
      if (pm->qbits && pm->use_qtab)
  	// custom qtab
  	comp_idx += store_array(comp+comp_idx, pm->qtab, 256);
  
      if (pm->pbits && pm->use_ptab)
  	// custom ptab
  	comp_idx += store_array(comp+comp_idx, pm->ptab, 1024);
  
      if (pm->dbits && pm->use_dtab)
  	// custom dtab
  	comp_idx += store_array(comp+comp_idx, pm->dtab, 256);
  
      return comp_idx;
  }
  
  static 
  int fqz_store_parameters(fqz_gparams *gp, unsigned char *comp) {
      int comp_idx = 0;
      comp[comp_idx++] = gp->vers; // Format number
  
      comp[comp_idx++] = gp->gflags;
  
      if (gp->gflags & GFLAG_MULTI_PARAM)
  	comp[comp_idx++] = gp->nparam;
  
      if (gp->gflags & GFLAG_HAVE_STAB) {
  	comp[comp_idx++] = gp->max_sel;
  	comp_idx += store_array(comp+comp_idx, gp->stab, 256);
      }
  
      int i;
      for (i = 0; i < gp->nparam; i++)
  	comp_idx += fqz_store_parameters1(&gp->p[i], comp+comp_idx);
  
      //fprintf(stderr, "Encoded %d bytes of param\n", comp_idx);
      return comp_idx;
  }
  
  // Choose a set of parameters based on quality statistics and
  // some predefined options (selected via "strat").
  static inline
  int fqz_pick_parameters(fqz_gparams *gp,
  			int vers,
  			int strat,
  			fqz_slice *s,
  			unsigned char *in,
  			size_t in_size) {
      //approx sqrt(delta), must be sequential
      int dsqr[] = {
  	0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3,
  	4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5,
  	5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  	6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7
      };
      uint32_t qhist[256] = {0};
  
      if (strat >= nstrats) strat = nstrats-1;
  
      // Start with 1 set of parameters.
      // FIXME: add support for multiple params later.
      memset(gp, 0, sizeof(*gp));
      gp->vers = FQZ_VERS;
  
      if (!(gp->p = calloc(1, sizeof(fqz_param))))
  	return -1;
      gp->nparam = 1;
      gp->max_sel = 0;
  
      if (vers == 3) // V3.0 doesn't store qual in original orientation
  	gp->gflags |= GFLAG_DO_REV;
  
      fqz_param *pm = gp->p;
  
      // Programmed strategies, which we then amend based on our
      // statistical analysis of the quality stream.
      pm->qbits  = strat_opts[strat][0];
      pm->qshift = strat_opts[strat][1];
      pm->pbits  = strat_opts[strat][2];
      pm->pshift = strat_opts[strat][3];
      pm->dbits  = strat_opts[strat][4];
      pm->dshift = strat_opts[strat][5];
      pm->qloc   = strat_opts[strat][6];
      pm->sloc   = strat_opts[strat][7];
      pm->ploc   = strat_opts[strat][8];
      pm->dloc   = strat_opts[strat][9];
  
      // Params for controlling behaviour here.
      pm->do_r2 = strat_opts[strat][10];
      pm->do_qa = strat_opts[strat][11];
  
      // Validity check input lengths and buffer size
      size_t tlen = 0, i;
      for (i = 0; i < s->num_records; i++) {
  	if (tlen + s->len[i] > in_size)
  	    // Oversized buffer
  	    s->len[i] = in_size - tlen;
  	tlen += s->len[i];
      }
      if (s->num_records > 0 && tlen < in_size)
  	// Undersized buffer
  	s->len[s->num_records-1] += in_size - tlen;
  
      // Quality metrics, for all recs
      fqz_qual_stats(s, in, in_size, pm, qhist, -1);
  
      pm->store_qmap = (pm->nsym <= 8 && pm->nsym*2 < pm->max_sym);
  
      // Check for fixed length.
      uint32_t first_len = s->len[0];
      for (i = 1; i < s->num_records; i++) {
  	if (s->len[i] != first_len)
  	    break;
      }
      pm->fixed_len = (i == s->num_records);
      pm->use_qtab = 0; // unused by current encoder
  
      if (strat >= nstrats-1)
  	goto manually_set; // used in TEST_MAIN for debugging
  
      if (pm->pshift < 0)
  	pm->pshift = MAX(0, log((double)s->len[0]/(1<<pm->pbits))/log(2)+.5);
  
      if (pm->nsym <= 4) {
  	// NovaSeq
  	pm->qshift = 2; // qmax 64, although we can store up to 256 if needed
  	if (in_size < 5000000) {
  	    pm->pbits =2;
  	    pm->pshift=5;
  	}
      } else if (pm->nsym <= 8) {
  	// HiSeqX
  	pm->qbits =MIN(pm->qbits,9);
  	pm->qshift=3;
  	if (in_size < 5000000)
  	    pm->qbits =6;
      }
  
      if (in_size < 300000) {
  	pm->qbits=pm->qshift;
  	pm->dbits=2;
      }
  
   manually_set:
  //    fprintf(stderr, "-x 0x%x%x%x%x%x%x%x%x%x%x%x%x\n",
  //	    pm->qbits, pm->qshift,
  //	    pm->pbits, pm->pshift,
  //	    pm->dbits, pm->dshift,
  //	    pm->qloc, pm->sloc, pm->ploc, pm->dloc,
  //	    pm->do_r2, pm->do_qa);
  
      for (i = 0; i < sizeof(dsqr)/sizeof(*dsqr); i++)
  	if (dsqr[i] > (1<<pm->dbits)-1)
  	    dsqr[i] = (1<<pm->dbits)-1;
  
      if (pm->store_qmap) {
  	int j;
  	for (i = j = 0; i < 256; i++)
  	    if (qhist[i])
  		pm->qmap[i] = j++;
  	    else
  		pm->qmap[i] = INT_MAX;
  	pm->max_sym = pm->nsym;
      } else {
  	pm->nsym = 255;
  	for (i = 0; i < 256; i++)
  	    pm->qmap[i] = i;
      }
      if (gp->max_sym < pm->max_sym)
  	gp->max_sym = pm->max_sym;
  
      // Produce ptab from pshift.
      if (pm->qbits) {
  	for (i = 0; i < 256; i++) {
  	    pm->qtab[i] = i; // 1:1
  
  	    // Alternative mappings:
  	    //qtab[i] = i > 30 ? MIN(max_sym,i)-15 : i/2;  // eg for 9827 BAM
  	}
  
      }
      pm->qmask = (1<<pm->qbits)-1;
  
      if (pm->pbits) {
  	for (i = 0; i < 1024; i++)
  	    pm->ptab[i] = MIN((1<<pm->pbits)-1, i>>pm->pshift);
  
  	// Alternatively via analysis of quality distributions we
  	// may select a bunch of positions that are special and
  	// have a non-uniform ptab[].
  	// Manual experimentation on a NovaSeq run saved 2.8% here.
      }
  
      if (pm->dbits) {
  	for (i = 0; i < 256; i++)
  	    pm->dtab[i] = dsqr[MIN(sizeof(dsqr)/sizeof(*dsqr)-1, i>>pm->dshift)];
      }
  
      pm->use_ptab = (pm->pbits > 0);
      pm->use_dtab = (pm->dbits > 0);
  
      pm->pflags =
  	(pm->use_qtab   ?PFLAG_HAVE_QTAB :0)|
  	(pm->use_dtab   ?PFLAG_HAVE_DTAB :0)|
  	(pm->use_ptab   ?PFLAG_HAVE_PTAB :0)|
  	(pm->do_sel     ?PFLAG_DO_SEL    :0)|
  	(pm->fixed_len  ?PFLAG_DO_LEN    :0)|
  	(pm->do_dedup   ?PFLAG_DO_DEDUP  :0)|
  	(pm->store_qmap ?PFLAG_HAVE_QMAP :0);
  
      gp->max_sel = 0;
      if (pm->do_sel) {
  	// 2 selectors values, but 1 parameter block.
  	// We'll use the sloc instead to encode the selector bits into
  	// the context.
  	gp->max_sel = 1; // indicator to check recs
  	gp->gflags |= GFLAG_HAVE_STAB;
  	// NB: stab is already all zero
      }
  
      if (gp->max_sel) {
  	int max = 0;
  	for (i = 0; i < s->num_records; i++) {
  	    if (max < (s->flags[i] >> 16))
  		max = (s->flags[i] >> 16);
  	}
  	gp->max_sel = max;
      }
  
      return 0;
  }
  
  static void fqz_free_parameters(fqz_gparams *gp) {
      if (gp && gp->p) free(gp->p);
  }
  
  static
  unsigned char *compress_block_fqz2f(int vers,
  				    int strat,
  				    fqz_slice *s,
  				    unsigned char *in,
  				    size_t in_size,
  				    size_t *out_size,
  				    fqz_gparams *gp) {
      fqz_gparams local_gp;
      int free_params = 0;
  
      unsigned int last = 0;
      size_t i, j;
      ssize_t rec = 0;
  
      int last_len = 0;
      int comp_idx = 0;
      RangeCoder rc;
  
      unsigned char *comp = (unsigned char *)malloc(in_size*1.1+100000);
      unsigned char *compe = comp + (size_t)(in_size*1.1+100000);
      if (!comp)
  	return NULL;
  
      // Pick and store params
      if (!gp) {
  	gp = &local_gp;
  	if (fqz_pick_parameters(gp, vers, strat, s, in, in_size) < 0)
  	    return NULL;
  	free_params = 1;
      }
  
      //dump_params(gp);
      comp_idx = var_put_u32(comp, compe, in_size);
      comp_idx += fqz_store_parameters(gp, comp+comp_idx);
  
      fqz_param *pm;
  
      // Optimise tables to remove shifts in loop (NB: cannot do this in next vers)
      for (j = 0; j < gp->nparam; j++) {
  	pm = &gp->p[j];
  
  	for (i = 0; i < 1024; i++)
  	    pm->ptab[i] <<= pm->ploc;
  
  	for (i = 0; i < 256; i++)
  	    pm->dtab[i] <<= pm->dloc;
      }
  
      // Create models and initialise range coder
      fqz_model model;
      if (fqz_create_models(&model, gp) < 0)
  	return NULL;
  
      RC_SetOutput(&rc, (char *)comp+comp_idx);
      RC_StartEncode(&rc);
  
      // For CRAM3.1, reverse upfront if needed
      pm = &gp->p[0];
      if (gp->gflags & GFLAG_DO_REV) {
  	i = rec = j = 0;
  	while (i < in_size) {
  	    int len = rec < s->num_records-1
  		? s->len[rec] : in_size - i;
  
  	    if (s->flags[rec] & FQZ_FREVERSE) {
  		// Reverse complement sequence - note: modifies buffer
  		int I,J;
  		unsigned char *cp = in+i;
  		for (I = 0, J = len-1; I < J; I++, J--) {
  		    unsigned char c;
  		    c = cp[I];
  		    cp[I] = cp[J];
  		    cp[J] = c;
  		}
  	    }
  
  	    i += len;
  	    rec++;
  	}
  	rec = 0;
      }
  
      fqz_state state = {0};
      pm = &gp->p[0];
      state.p = 0;
      state.first_len = 1;
      int x;
  
      for (i = 0; i < in_size; i++) {
  	if (state.p == 0) {
  	    if (pm->do_sel || (gp->gflags & GFLAG_MULTI_PARAM)) {
  		state.s = rec < s->num_records
  		    ? s->flags[rec] >> 16 // reuse spare bits
  		    : 0;
  		SIMPLE_MODEL(256,_encodeSymbol)(&model.sel, &rc, state.s);
  		//fprintf(stderr, "State %d\n", state.s);
  	    } else {
  		state.s = 0;
  	    }
  	    x = (gp->gflags & GFLAG_HAVE_STAB) ? gp->stab[state.s] : state.s;
  	    pm = &gp->p[x];
  
  	    //fprintf(stderr, "sel %d param %d\n", state.s, x);
  
  	    int len = s->len[rec];
  	    if (!pm->fixed_len || state.first_len) {
  		SIMPLE_MODEL(256,_encodeSymbol)(&model.len[0], &rc, (len>> 0) & 0xff);
  		SIMPLE_MODEL(256,_encodeSymbol)(&model.len[1], &rc, (len>> 8) & 0xff);
  		SIMPLE_MODEL(256,_encodeSymbol)(&model.len[2], &rc, (len>>16) & 0xff);
  		SIMPLE_MODEL(256,_encodeSymbol)(&model.len[3], &rc, (len>>24) & 0xff);
  		//fprintf(stderr, "Len %d\n", len);
  		state.first_len = 0;
  	    }
  
  	    if (gp->gflags & GFLAG_DO_REV) {
  		// no need to reverse complement for V4.0 as the core format
  		// already has this feature.
  		if (s->flags[rec] & FQZ_FREVERSE)
  		    SIMPLE_MODEL(2,_encodeSymbol)(&model.revcomp, &rc, 1);
  		else
  		    SIMPLE_MODEL(2,_encodeSymbol)(&model.revcomp, &rc, 0);
  		//fprintf(stderr, "Rev %d\n", (s->flags[rec] & FQZ_FREVERSE) ? 1 : 0);
  	    }
  
  	    rec++;
  
  	    state.qtot = 0;
  	    state.qlen = 0;
  
  	    state.p = len;
  	    state.add_d = 0;
  	    state.delta = 0;
  	    state.qctx = 0;
  	    state.prevq = 0;
  
  	    last = pm->context;
  
  	    if (pm->do_dedup) {
  		// Possible dup of previous read?
  		if (i && len == last_len && !memcmp(in+i-last_len, in+i, len)) {
  		    SIMPLE_MODEL(2,_encodeSymbol)(&model.dup, &rc, 1);
  		    i += len-1;
  		    state.p = 0;
  		    //fprintf(stderr, "Dup 1\n");
  		    continue;
  		} else {
  		    SIMPLE_MODEL(2,_encodeSymbol)(&model.dup, &rc, 0);
  		    //fprintf(stderr, "Dup 0\n");
  		}
  
  		last_len = len;
  	    }
  	}
  
  	unsigned char q = in[i];
  	unsigned char qm = pm->qmap[q];
  
  	SIMPLE_MODEL(QMAX,_encodeSymbol)(&model.qual[last], &rc, qm);
  	//fprintf(stderr, "Sym %d with ctx %04x delta %d prevq %d q %d\n", qm, last, state.delta, state.prevq, qm);
  	//fprintf(stderr, "pos=%d, delta=%d\n", state.p, state.delta);
  	last = fqz_update_ctx(pm, &state, qm);
      }
  
      RC_FinishEncode(&rc);
  
      // For CRAM3.1, undo our earlier reversal step
      if (gp->gflags & GFLAG_DO_REV) {
  	i = rec = j = 0;
  	while (i < in_size) {
  	    int len = rec < s->num_records-1
  		? s->len[rec]
  		: in_size - i;
  
  	    if (s->flags[rec] & FQZ_FREVERSE) {
  		// Reverse complement sequence - note: modifies buffer
  		int I,J;
  		unsigned char *cp = in+i;
  		for (I = 0, J = len-1; I < J; I++, J--) {
  		    unsigned char c;
  		    c = cp[I];
  		    cp[I] = cp[J];
  		    cp[J] = c;
  		}
  	    }
  
  	    i += len;
  	    rec++;
  	}
      }
  
      // Clear selector abuse of flags
      for (rec = 0; rec < s->num_records; rec++)
  	s->flags[rec] &= 0xffff;
  
      *out_size = comp_idx + RC_OutSize(&rc);
      //fprintf(stderr, "%d -> %d\n", (int)in_size, (int)*out_size);
  
      fqz_destroy_models(&model);
      if (free_params)
  	fqz_free_parameters(gp);
  
      return comp;
  }
  
  // Read fqz paramaters.
  //
  // FIXME: pass in and check in_size.
  //
  // Returns number of bytes read on success,
  //         -1 on failure.
  static inline
  int fqz_read_parameters1(fqz_param *pm, unsigned char *in, size_t in_size) {
      int in_idx = 0;
      size_t i;
  
      if (in_size < 7)
  	return -1;
  
      // Starting context
      pm->context = in[in_idx] + (in[in_idx+1]<<8);
      in_idx += 2;
  
      // Bit flags
      pm->pflags     = in[in_idx++];
      pm->use_qtab   = pm->pflags & PFLAG_HAVE_QTAB;
      pm->use_dtab   = pm->pflags & PFLAG_HAVE_DTAB;
      pm->use_ptab   = pm->pflags & PFLAG_HAVE_PTAB;
      pm->do_sel     = pm->pflags & PFLAG_DO_SEL;
      pm->fixed_len  = pm->pflags & PFLAG_DO_LEN;
      pm->do_dedup   = pm->pflags & PFLAG_DO_DEDUP;
      pm->store_qmap = pm->pflags & PFLAG_HAVE_QMAP;
      pm->max_sym    = in[in_idx++];
  
      // Sub-context sizes and locations
      pm->qbits      = in[in_idx]>>4;
      pm->qmask      = (1<<pm->qbits)-1;
      pm->qshift     = in[in_idx++]&15;
      pm->qloc       = in[in_idx]>>4;
      pm->sloc       = in[in_idx++]&15;
      pm->ploc       = in[in_idx]>>4;
      pm->dloc       = in[in_idx++]&15;
  
      // Maps and tables
      if (pm->store_qmap) {
  	for (i = 0; i < 256; i++) pm->qmap[i] = INT_MAX; // so dump_map works
  	if (in_idx + pm->max_sym > in_size)
  	    return -1;
  	for (i = 0; i < pm->max_sym; i++)
  	    pm->qmap[i] = in[in_idx++];
      } else {
  	for (i = 0; i < 256; i++)
  	    pm->qmap[i] = i;
      }
  
      if (pm->qbits) {
  	if (pm->use_qtab)
  	    in_idx += read_array(in+in_idx, in_size-in_idx, pm->qtab, 256);
  	else
  	    for (i = 0; i < 256; i++)
  		pm->qtab[i] = i;
      }
  
      if (pm->use_ptab)
  	in_idx += read_array(in+in_idx, in_size-in_idx, pm->ptab, 1024);
      else
  	for (i = 0; i < 1024; i++)
  	    pm->ptab[i] = 0;
  
      if (pm->use_dtab)
          in_idx += read_array(in+in_idx, in_size-in_idx, pm->dtab, 256);
      else
  	for (i = 0; i < 256; i++)
  	    pm->dtab[i] = 0;
  
      return in_idx;
  }
  
  static
  int fqz_read_parameters(fqz_gparams *gp, unsigned char *in, size_t in_size) {
      int in_idx = 0;
      int i;
  
      if (in_size < 10)
  	return -1;
  
      // Format version
      gp->vers = in[in_idx++];
      if (gp->vers != FQZ_VERS)
  	return -1;
  
      // Global glags
      gp->gflags = in[in_idx++];
  
      // Number of param blocks and param selector details
      gp->nparam = (gp->gflags & GFLAG_MULTI_PARAM) ? in[in_idx++] : 1;
      if (gp->nparam <= 0)
  	return -1;
      gp->max_sel = gp->nparam > 1 ? gp->nparam : 0;
  
      if (gp->gflags & GFLAG_HAVE_STAB) {
  	gp->max_sel = in[in_idx++];
  	in_idx += read_array(in+in_idx, in_size-in_idx, gp->stab, 256);
      } else {
  	for (i = 0; i < gp->nparam; i++)
  	    gp->stab[i] = i;
  	for (; i < 256; i++)
  	    gp->stab[i] = gp->nparam-1;
      }
  
      // Load the individual parameter locks
      if (!(gp->p = malloc(gp->nparam * sizeof(*gp->p))))
  	return -1;
  
      gp->max_sym = 0;
      for (i = 0; i < gp->nparam; i++) {
  	int e = fqz_read_parameters1(&gp->p[i], in + in_idx, in_size-in_idx);
  	if (e < 0)
  	    goto err;
  	in_idx += e;
  
  	if (gp->max_sym < gp->p[i].max_sym)
  	    gp->max_sym = gp->p[i].max_sym;
      }
  
      //fprintf(stderr, "Decoded %d bytes of param\n", in_idx);
      return in_idx;
  
   err:
      fqz_free_parameters(gp);
      gp->nparam = 0;
      return -1;
  }
  
  static
  unsigned char *uncompress_block_fqz2f(fqz_slice *s,
  				      unsigned char *in,
  				      size_t in_size,
  				      size_t *out_size,
  				      int *lengths,
  				      int nlengths) {
      fqz_gparams gp;
      fqz_param *pm;
      char *rev_a = NULL;
      int *len_a = NULL;
      memset(&gp, 0, sizeof(gp));
  
      uint32_t len;
      ssize_t i, rec = 0, in_idx;
      in_idx = var_get_u32(in, in+in_size, &len);
      *out_size = len;
  
  #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
      if (len > 100000)
  	return NULL;
  #endif
  
      unsigned char *uncomp = NULL;
      RangeCoder rc;
      unsigned int last = 0;
  
      // Decode parameter blocks
      if ((i = fqz_read_parameters(&gp, in+in_idx, in_size-in_idx)) < 0)
  	return NULL;
      //dump_params(&gp);
      in_idx += i;
  
      // Optimisations to remove shifts from main loop
      for (i = 0; i < gp.nparam; i++) {
  	int j;
  	pm = &gp.p[i];
  	for (j = 0; j < 1024; j++)
  	    pm->ptab[j] <<= pm->ploc;
  	for (j = 0; j < 256; j++)
  	    pm->dtab[j] <<= pm->dloc;
      }
  
      // Initialise models and entropy coder
      fqz_model model;
      if (fqz_create_models(&model, &gp) < 0)
  	return NULL;
  
      RC_SetInput(&rc, (char *)in+in_idx, (char *)in+in_size);
      RC_StartDecode(&rc);
  
  
      // Allocate buffers
      uncomp = (unsigned char *)malloc(*out_size);
      if (!uncomp)
  	goto err;
  
      int nrec = 1000;
      rev_a = malloc(nrec);
      len_a = malloc(nrec * sizeof(int));
      if (!rev_a || !len_a)
  	goto err;
  
      // Main decode loop
      fqz_state state;
      state.delta = 0;
      state.prevq = 0;
      state.qctx = 0;
      state.p = 0;
      state.s = 0;
      state.first_len = 1;
  
      int rev = 0;
      int last_len = 0;
      int x = 0;
      pm = &gp.p[x];
      for (rec = i = 0; i < len; i++) {
  	if (rec >= nrec) {
  	    nrec *= 2;
  	    rev_a = realloc(rev_a, nrec);
  	    len_a = realloc(len_a, nrec*sizeof(int));
  	    if (!rev_a || !len_a)
  		goto err;
  	}
  
  	if (state.p == 0) {
  	    // New record
  	    if (pm->do_sel) {
  		state.s = SIMPLE_MODEL(256,_decodeSymbol)(&model.sel, &rc);
  		//fprintf(stderr, "State %d\n", state.s);
  	    } else {
  		state.s = 0;
  	    }
  	    x = (gp.gflags & GFLAG_HAVE_STAB) ? gp.stab[MIN(255, state.s)] : state.s;
  	    if (x >= gp.nparam)
  		goto err;
  	    pm = &gp.p[x];
  
  	    unsigned int len = last_len;
  	    if (!pm->fixed_len || state.first_len) {
  		len  = SIMPLE_MODEL(256,_decodeSymbol)(&model.len[0], &rc);
  		len |= SIMPLE_MODEL(256,_decodeSymbol)(&model.len[1], &rc)<<8;
  		len |= SIMPLE_MODEL(256,_decodeSymbol)(&model.len[2], &rc)<<16;
  		len |= ((unsigned)SIMPLE_MODEL(256,_decodeSymbol)(&model.len[3], &rc))<<24;
  		//fprintf(stderr, "Len %d\n", len);
  		state.first_len = 0;
  		last_len = len;
  	    }
  	    if (len > *out_size-i || len <= 0)
  		goto err;
  
  	    if (lengths && rec < nlengths)
  		lengths[rec] = len;
  
  	    if (gp.gflags & GFLAG_DO_REV) {
  		rev = SIMPLE_MODEL(2,_decodeSymbol)(&model.revcomp, &rc);
  		//fprintf(stderr, "rev %d\n", rev);
  		rev_a[rec] = rev;
  		len_a[rec] = len;
  	    }
  
  	    if (pm->do_dedup) {
  		if (SIMPLE_MODEL(2,_decodeSymbol)(&model.dup, &rc)) {
  		    //fprintf(stderr, "Dup 1\n");
  		    // Dup of last line
  		    if (len > i)
  			goto err;
  		    memcpy(uncomp+i, uncomp+i-len, len);
  		    i += len-1;
  		    state.p = 0;
  		    rec++;
  		    continue;
  		} else {
  		    //fprintf(stderr, "Dup 0\n");
  		}
  	    }
  
  	    rec++;
  
  	    state.p = len;
  	    state.add_d = 0;
  	    state.delta = 0;
  	    state.prevq = 0;
  	    state.qctx = 0;
  
  	    last = pm->context;
  	}
  
  	// Decode and output quality
  	unsigned char Q = SIMPLE_MODEL(QMAX,_decodeSymbol)(&model.qual[last], &rc);
  	unsigned char q = pm->qmap[Q];
  	//fprintf(stderr, "Sym %d with ctx %04x delta %d prevq %d q %d\n", Q, last, state.delta, state.prevq, Q);
          uncomp[i] = q;
  
  	// Compute new quality context
  	last = fqz_update_ctx(pm, &state, Q);
      }
  
      if (rec >= nrec) {
  	nrec *= 2;
  	rev_a = realloc(rev_a, nrec);
  	len_a = realloc(len_a, nrec*sizeof(int));
  	if (!rev_a || !len_a)
  	    goto err;
      }
      rev_a[rec] = rev;
      len_a[rec] = len;
  
      if (gp.gflags & GFLAG_DO_REV) {
  	for (i = rec = 0; i < len && rec < nrec; i += len_a[rec++]) {
  	    if (!rev_a[rec])
  		continue;
  
  	    int I, J;
  	    unsigned char *cp = uncomp+i;
  	    for (I = 0, J = len_a[rec]-1; I < J; I++, J--) {
  		unsigned char c;
  		c = cp[I];
  		cp[I] = cp[J];
  		cp[J] = c;
  	    }
  	}
      }
  
      RC_FinishDecode(&rc);
      fqz_destroy_models(&model);
      //free(model.qual);
      free(rev_a);
      free(len_a);
      fqz_free_parameters(&gp);
  
  #ifdef TEST_MAIN
      s->num_records = rec;
  #endif
  
      return uncomp;
  
   err:
      fqz_destroy_models(&model);
      free(rev_a);
      free(len_a);
      fqz_free_parameters(&gp);
      free(uncomp);
  
      return NULL;
  }
  
  char *fqz_compress(int vers, fqz_slice *s, char *in, size_t uncomp_size,
  		   size_t *comp_size, int strat, fqz_gparams *gp) {
      return (char *)compress_block_fqz2f(vers, strat, s, (unsigned char *)in,
  					uncomp_size, comp_size, gp);
  }
  
  char *fqz_decompress(char *in, size_t comp_size, size_t *uncomp_size,
  		     int *lengths, int nlengths) {
      return (char *)uncompress_block_fqz2f(NULL, (unsigned char *)in,
  					  comp_size, uncomp_size, lengths, nlengths);
  }
  