void nmeeg_cleanup_compute_fitting_matrix(int np, int ne, double *P, double *Q);
  /* Stores into {Q} the inverse of {P*P'}.  Also checks whether the inverse is 
    valid. */

void nmeeg_cleanup_print_channel_stats(FILE *wr, int nt, int nc, int ne, double **val, char *chname[]);
  /* Gathers and prints the per-channel statistics of the data set {val[0..nt-1][0..nc-1]}.
    Assumes that the first {ne} channels are electrode voltages (real or synthetic). */

string_vec_t nmeeg_cleanup_parse_chnames(argparser_t *pp);
  /* Parses a list of channel names from the command line, stopping at the next 
     keyword (any argument that starts with "-"). */


#include <neuromat_eeg.h>
#include <neuromat_eeg_io.h>
#include <neuromat_eeg_header.h>
#include <neuromat_eeg_channel_stats.h>


void nclup_extract_blink_patterns
  ( int nt, 
    int ne, 
    double **vin, 
    int np, 
    double *P, 
    double *Q, 
    double **vot
  );
  /* Converts {ne} electrode potential channels from the dataset {vin}
    to {ne+np} channels in the new dataset {vot}, using the pattern matrix {P}.

    Assumes that {vin[it]} and {vot[it]} are corresponding input and
    output data frames for each {it} in {0..nt-1}; which have at least
    {ne} and {ne+np} electrode values, respectively (possibly followed
    by other channels).
    
    Also assumes that {P} is a linearized array with {np} rows and
    {ne} columns, such that {P[ip*ne + ie]} is the value of
    electrode {ie} in pattern {ip}.  Assumes that {Q} is a linearized
    matrix with {np} rows and {np} columns, the inverse of {P*P'}
    where {P'} is the transpose of {P}. 
    
    For each input frame {vin[it]}, the program considers its first {ne}
    samples to be a row vector {vi[0..ne-1]}, and computes {np}
    coefficients {vo[0..np-1]} by the product {Q*P*vi'}. These
    coefficients are stored into the output frame {vot[it]} as channels
    {ne..ne+np-1}.
    
    Then the program subtracts from the vector {vi} the fitted components, namely {P*vo'},
    and stores that residual into {vot[it]} as channels {0..ne-1}. */

void nclup_recombine_patterns(int nt, int ne, int np, double **vsp, int mp, int ip[], double *P, double **vrc);
  /* Separates a specified pattern from an EEG dataset.
  
    Assumes that {vsp} is an EEG dataset with {nt} frames {vsp[0..nt-1]},
    each frame containing {ne} residual electrode readings followed by {np} pattern
    coefficients (and possibly other channels), as extracted by
    {nclup_extract_patterns}. Assumes that {P} is a linearized array with
    {np} rows and {ne} columns, where {P[ip*ne + ie]} is the value of
    electrode {ie} in pattern {ip}.  Assumes that {ip[0..mp-1]} are 
    integers in the range {0..np-1}, Finally assumes that {vrc} is an 
    array of {nt} data frames, each having {nc} channels.
    
    The procedure stores into channels {0..ne-1} of each frame
    {vrc[it]} the EEG electrode potentials corresponding to the
    linear combination of rows {ip[0..mp-1]} of {P}, with their
    respective coefficients as found by {nclup_extract_patterns};
    namely, sets
    
      {vrc[it][ie] = SUM{ vsp[it][ne+ip[k]]*P[ip[k]*ne + ie] : k \in 0..mp-1}}
      
    for {ie}in {0..ne-1}.  The remaining {nc-ne} channels of each frame
    {vrc[it]} are left unchanged. */
          
void nclup_append_norm_channel(int nt, int nc, double **val, int mc, int ic[]);
  /* Assumes that {val} is an array of {nt} frames, each with at least {nc+1} channels
    of which the first {nc} are defined.  Also assumes that {ic[0..mc-1]} is
    a list of indices in {0..nc-1}.  Sets channel {nc} of every frame {v=val[it]}
    to the Euclidean norm of the values {v[ic[0..mc-1]}. */

void nclup_delete_channels(int nt, int *ncP, int *neP, double **val, char *chname[], bool_t keep[]);
  /* Assumes that {val} is an array with {nt} frames, each with {nc} channels,
    of thich the first {ne} are electrodes (real or synthetic);
    where {nc} is the input value of {*ncP}, and {ne} is the 
    input value of {*neP}.  Also assumes that {chname[0..nc-1]}
    are the (unshared) names of those channels, and that {keep} is an array of {nc}
    booleans.  Deletes from {chname} and from every frame of {val} 
    those channels with indices {ic} such taht {keep[ic]} is false.
    Compacts the remaining channels, keeping their order.  Then 
    updates {*ncP} and {*neP}with the number of surviving channels
    and surviving electrodes, respectively. */

    neuromat_eeg_header_t *h = NULL; /* Parameters from the input file header. */
    double **vin = NULL;
    nclup_read_run_data(stdin, &h, &vin);
    int nt = h->nt;
    int nc_in = h->nc;
    int ne_in = h->ne;
    
    fprintf(stderr, "input channel statistics\n");
    nclup_print_channel_stats(stderr, nt, nc_in, ne_in, vin, h->chname);
    
    int np = o->pattern.ne; /* Number of pattern frames: */
    int nr = o->norm.ne;    /* Number of norm channels. */
    int nd = o->delete.ne;  /* Number of channels to be deleted. */
    
    int it, ie;

    /* Allocate the output dataset will all channels (before channel deletion): */
    int nc_tp = nc_in + np + nr; /* Number of temporary channels. */
    double **vot = notnull(malloc(nt*sizeof(double*)), "no mem"); /* Index {[0..nt-1][0..nc_tp-1]}. */
    for (it = 0; it < nt; it++) { vot[it] = notnull(malloc(nc_tp*sizeof(double)), "no mem"); }
        
    /* Allocate the list of output channel names (with space also for channels to be deleted): */
    char *chname_ot[nc_tp];
    
    /* Initialize the list of output channels with the input electrode names: */
    int nc_ot = ne_in; /* Number of channels in output file. */
    int ne_ot = ne_in; /* Number of electrodes (real or synthetic) in output file. */
    for (ie = 0; ie < ne_in; ie++) { chname_ot[ie] = txtcat(h->chname[ie], ""); }
    
    if (np > 0)
      { /* Read the patterns into {P[0..np*ne-1]}, save their names in {chname_ot[ne_in..ne_in+np-1]}: */
        demand(np <= ne_in, "too many patterns");
        double *P = rmxn_alloc(np,ne_in); /* Each row is a pattern frame. */
        char **patnames = &(chname_ot[nc_ot]); /* Section of {chname_ot} with the pattern names. */
        int ip;
        for (ip = 0; ip < np; ip++)
          { pat_spec_t *Sip = &(o->pattern.e[ip]);
            double *Pip = &(P[ip*ne_in]);
            fprintf(stderr, "reading pattern %s from file %s\n", Sip->patname, Sip->fname);
            nclup_read_pattern(Sip->fname, ne_in, h->chname, o->normalize, Pip);
            /* Check for repeated names: */
            int jp = neuromat_eeg_find_channel_by_name(Sip->patname, 0, nc_ot-1, chname_ot, FALSE);
            demand(jp == -1, "repeated pattern or electrode name");
            assert(nc_ot == ne_in + ip); /* Paranoia. */
            patnames[ip] = txtcat(Sip->patname, "");
            nc_ot++;
            ne_ot++;
          }
          
        fprintf(stderr, "== debug 0 ==\n");
        int ic; for (ic = 0; ic < nc_ot; ic++) { fprintf(stderr, "channel %3d = %-4s\n", ic, chname_ot[ic]); }
          
        /* Compute the array {Q[0..np*np-1]}, the inverse of {P*P'}: */
        double *Q = rmxn_alloc(np,np); /* The array {(P*P')^-1}. */
        nclup_compute_fitting_matrix(np,ne_in,P,Q);

        fprintf(stderr, "== debug 0.4 ==\n");
        for (ic = 0; ic < nc_ot; ic++) { fprintf(stderr, "channel %3d = %-4s\n", ic, chname_ot[ic]); }
        
        /* Extract the patterns: */
        nclup_extract_patterns(nt, ne_in, vin, np, P, Q, vot);

        fprintf(stderr, "== debug 0.5 ==\n");
        for (ic = 0; ic < nc_ot; ic++) { fprintf(stderr, "channel %3d = %-4s\n", ic, chname_ot[ic]); }

        /* Write the fitted pattern combinations: */
        int nw = o->writeComp.ne;
        if (nw > 0)
          { int iw;
            for (iw = 0; iw < nw; iw++)
              { wcomp_spec_t *Siw = &(o->writeComp.e[iw]);  /* Spec of one output pattern combination file. */
                fprintf(stderr, "writing component %s to file %s\n", Siw->compname, Siw->fname);
                int mp = Siw->patnames.ne; /* Number of patterns to combine into this file. */
                int ip[mp]; /* Indices of patterns to combine. */
                int k;
                for (k = 0; k < mp; k++) 
                  { ip[k] = neuromat_eeg_find_channel_by_name(Siw->patnames.e[k], 0, np-1, patnames, TRUE); }
                nclup_recombine_patterns(nt, ne_in, np, vot, mp, ip, P, vin);
                h->component = Siw->compname;
                FILE *wr = open_write(Siw->fname, TRUE);
                nclup_write_run_data(wr, nt, nc_in, vin, h);
                fclose(wr);
              }
          }
      }

    fprintf(stderr, "== debug 1 ==\n");
    int ic; for (ic = 0; ic < nc_ot; ic++) { fprintf(stderr, "channel %3d = %-4s\n", ic, chname_ot[ic]); }

    if (nr > 0)
      { /* Add the norm channels to {vot}: */
        int ir;
        for (ir = 0; ir < nr; ir++)
          { norm_spec_t *Sir = &(o->norm.e[ir]);;
            int mc = Sir->chname.ne; /* Number of channels to be combined into this norm */
            int ic[mc]; /* Indices of channels to combine. */
            int k;
            for (k = 0; k < mc; k++) 
              { ic[k] = neuromat_eeg_find_channel_by_name(Sir->chname.e[k], 0, nc_ot-1, chname_ot, TRUE); }
            nclup_append_norm_channel(nt, nc_ot, vot, mc, ic);
            int jc = neuromat_eeg_find_channel_by_name(Sir->normname, 0, nc_ot-1, chname_ot, FALSE);
            demand(jc < 0, "repeated norm, electrode, or pattern name");
            chname_ot[nc_ot] = txtcat(Sir->normname, "");
            nc_ot++;
            ne_ot++;
          }
      }

    assert(ne_ot == nc_ot);

    fprintf(stderr, "== debug 2 ==\n");
    for (ic = 0; ic < nc_ot; ic++) { fprintf(stderr, "channel %3d = %-4s\n", ic, chname_ot[ic]); }

    if (nc_in > ne_in)
      { /* Append marker/trigger (non-electrode) channels: */
        int it, ic, jc;
        for (it = 0; it < nt; it++) 
          { double *vi = vin[it];
            double *vo = vot[it];
            for (ic = ne_in, jc = nc_ot; ic < nc_in; ic++, jc++) { vo[jc] = vi[ic]; }
          }
        for (ic = ne_in; ic < nc_in; ic++)
          { char *mkname = h->chname[ic];
            int kc = neuromat_eeg_find_channel_by_name(mkname, 0, nc_ot-1, chname_ot, FALSE);
            demand(kc < 0, "collision with marker channel name");
            chname_ot[nc_ot] = txtcat(mkname, "");
            nc_ot++;
          }
      }

    assert(nc_ot == nc_tp);

    fprintf(stderr, "== debug 3 ==\n");
    for (ic = 0; ic < nc_ot; ic++) { fprintf(stderr, "channel %3d = %-4s\n", ic, chname_ot[ic]); }

    if (nd > 0)
      { /* Delete channels: */
        bool_t keep[nc_ot]; /* TRUE for channels to keep. */
        int ic, id;
        for (ic = 0; ic < nc_in; ic++) { keep[ic] = TRUE; }
        for (id = 0; id < nd; id++)
          { char *delname = o->delete.e[id];
            int ick = neuromat_eeg_find_channel_by_name(delname, 0, nc_ot-1, chname_ot, TRUE);
            demand(keep[ick], "duplicate channel in delete list");
            keep[ick] = FALSE;
          }
        nclup_delete_channels(nt, &nc_ot, &ne_ot, vot, chname_ot, keep);
      }

    fprintf(stderr, "== debug 4 ==\n");
    for (ic = 0; ic < nc_ot; ic++) { fprintf(stderr, "channel %3d = %-4s\n", ic, chname_ot[ic]); }

    fprintf(stderr, "output channel statistics\n");
    nclup_print_channel_stats(stderr, nt, nc_ot, ne_ot, vot, chname_ot);
    
    /* Write the output EEG dataset: */
    h->nc = nc_ot; /* Principal component strengths and triggers. */
    h->chname = chname_ot;
    h->ne = ne_ot;
    h->component = NULL;
    nclup_write_run_data(stdout, nt, nc_ot, vot, h);

    /* Dont's bother to free storage: */

void nclup_delete_channels(int nt, int *ncP, int *neP, double **val, char *chname[], bool_t keep[])
  {
    int nc_in = (*ncP);
    int ne_in = (*neP);
    int ic, nc_ot, ne_ot;
    
    /* Delete the channels in all frames: */
    int it;
    for (it = 0; it < nt; it++)
      { double *v = val[it];
        for (ic = 0, nc_ot = 0; ic < nc_in; ic++)
          { if (keep[ic]) { v[nc_ot] = v[ic]; nc_ot++; } }
      }
    
    /* Delete the channel names: */
    for (ic = 0, nc_ot = 0, ne_ot = 0; ic < nc_in; ic++)
      { if (keep[ic]) 
          { chname[nc_ot] = chname[ic]; 
            if (ic < ne_in) { ne_ot++; }
            nc_ot++;
          } 
        else
          { free(chname[ic]); }
      }

    /* Update channel and electrode counts: */
    (*ncP) = nc_ot;
    (*neP) = ne_ot;
  }

void nclup_read_pattern(char *fname, int ne, char *chname[], bool_t normalize, double pat[])
  {
    FILE *rd = open_read(fname, TRUE); 
    neuromat_eeg_header_t *h = NULL;
    double **val = NULL;
    nclup_read_run_data(rd, &h, &val);
    demand(h->ne == ne, "electrodes counts in pattern and file do not match");
    demand(h->nt == 1, "pattern should have a single data frame");
    demand (h->chname != NULL, "missing channel names in pattern header");
    double sum2 = 0;
    int ie;
    for (ie = 0; ie < ne; ie++) 
      { demand(strcmp(chname[ie],h->chname[ie])== 0, "channel names in pattern and file do not match"); 
        double vi = val[0][ie];
        pat[ie] = vi;
        sum2 += vi*vi;
      }
    double norm = sqrt(sum2);
    fprintf(stderr, "pattern norm = %12.6f\n", norm);
    demand(norm > 1.0e-5, "pattern norm too small\n");
    if (normalize)
      { fprintf(stderr, "rescaling to unit rms norm\n");
        rn_scale(ne, 1/norm, pat, pat);
      }
    free(val);
    free(h);
    fclose(rd);
  }
    
void nmeeg_cleanup_append_norm_channel(int nt, int nc, double **val, int mc, int ic[]) 
  {
    int it, k;
    for (k = 0; k < mc; k++) { demand((ic[k] >= 0) && (ic[k] < nc), "invalid channel index"); }
    for (it = 0; it < nt; it++)
      { double *v = val[it];
        assert(v != NULL);
        int k;
        double sum2 = 0;
        for (k = 0; k < mc; k++)
          { double vk = v[ic[k]];
            sum2 += vk*vk;
          }
        v[nc] = sqrt(sum2);
      }
  }
           
void nmeeg_cleanup_compute_fitting_matrix(int np, int ne, double *P, double *Q)
  {
    double *R = rmxn_alloc(np,np); /* The array {P*P'}. */
    rmxn_mul_tr(np, np, ne, P, P, R);
    rmxn_inv(np,R,Q);
    /* Check the inverse: */
    double *S = rmxn_alloc(np,np); /* Should be the identity. */
    rmxn_mul(np,np,np,R,Q,S);
    int i, j;
    for (i = 0; i < np; i++)
      { for (j = 0; j < np; j++)
          { double Sij = S[i*np + j];
            double Iij = (i == j ? 1.0 : 0.0);
            if (fabs(Sij - Iij) > 1.0e-6)
              { fprintf(stderr, "  S[%d,%d] = %24.16e\n", i, j, Sij);
                demand(FALSE, "inversion failed, aborted");
              }
          }
      }
    free(R);
    free(S);
  }
    


void nmeeg_cleanup_extract_patterns(int nt, int ne, double **vin, int np, double *P, double *Q, double **vot)
  {
    demand(ne >= 0, "invalid {ne}");
    demand(np >= 0, "invalid {np}");
    demand(np <= ne, "more components than electrodes");
    int it;
    for (it = 0; it < nt; it++)
      { double *vi = vin[it]; /* Original electrode values {vi[0..ne-1]}. */
        double *cf = &(vot[it][ne]); /* Coefficients of each pattern {cf[0..np-1]}. */
        double *vo = vot[it]; /* Residual is {vo[0..ne-1]}. */
        /* Map channels {valt[0..ne-1]} by first {np} rows of {P} to {vpct[0..np-1]}: */
        double b[np];
        rmxn_map_col(np, ne, P, vi, b);
        rmxn_map_col(np, np, Q, b, cf);
        /* Compute residual: */
        double vr[ne]; /* Reconstructed frame from fitted patterns. */
        rmxn_tr_mul(np, ne, 1, P, cf, vr);
        rn_sub(ne, vi, vr, vo);
      }
  }

void nmeeg_cleanup_recombine_patterns(int nt, int ne, int np, double **vsp, int mp, int ip[], double *P, double **vrc)
  {
    demand((0 <= np) && (np <= ne), "invalid {np,ne}");
    int it;
    for (it = 0; it < nt; it++)
      { double *cf = &(vsp[it][ne]); /* Fitted pattern coefficients. */
        double *vr = vrc[it];        /* Reconstructed electrodes. */
        /* Combine the requested components: */
        rn_zero(ne, vr);
        int k;
        for (k = 0; k < mp; k++)
          { int ipk = ip[k];
            demand((0 <= ipk) && (ipk < np), "invalid pattern index");
            double *Pi = &(P[ipk*ne]); /* Selected pattern. */
            double cfi = cf[ipk]; /* Its coefficient. */
            rn_mix_in (ne, cfi, Pi, vr);
          }
      }
  }

string_vec_t nmeeg_cleanup_parse_chnames(argparser_t *pp)
  { 
    string_vec_t chname = string_vec_new(10);
    int nc = 0;
    while (argparser_next_is_non_keyword(pp))
      { char *name = argparser_get_next_non_keyword(pp);
        string_vec_expand(&chname, nc);
        chname.e[nc] = name;
        nc++;
      }
    string_vec_trim(&chname, nc);
    return chname;
  }


void nmeeg_cleanup_print_channel_stats(FILE *wr, int nt, int nc, int ne, double **val, char *chname[])
  {
    /* Compute the channel means and variances over the selected frames: */
    neuromat_eeg_channel_stats_t *st = neuromat_eeg_channel_stats_new(nc);
    neuromat_eeg_channel_stats_t *stg = neuromat_eeg_channel_stats_new(1);
    double eps = 0.01; /* Assumed uncertainty of measurement (µV). */
    neuromat_eeg_channel_stats_gather_all(nt, nc, val, NULL, eps, st, ne, stg);
    fprintf(wr, "  --- channel statistics ---\n");
    neuromat_eeg_channel_stats_print_all(wr, 2, nc, chname, FALSE, st, ne, stg);
    fprintf(wr, "\n");
    free(st);
    free(stg);
  }

vec_typeimpl(pat_spec_vec_t,pat_spec_vec,pat_spec_t);
vec_typeimpl(norm_spec_vec_t,norm_spec_vec,norm_spec_t);
vec_typeimpl(wcomp_spec_vec_t,wcomp_spec_vec,wcomp_spec_t);