Does the work when creating peptide and protein quant lookups. This
    loops through storing options and parses columns, passing on to the
    storing functions
7395
def create_pep_protein_quant_lookup(fns, pqdb, poolnames, featcolnr, patterns,
                                    storefuns, isobqcolpattern=None,
                                    psmnrpattern=None):
    """"""
    tablefn_map = create_tablefn_map(fns, pqdb, poolnames)
    feat_map = pqdb.get_feature_map()
    for pattern, storefun in zip(patterns, storefuns):
        if pattern is None:
            continue
        colmap = get_colmap(fns, pattern, single_col=True)
        if colmap:
            store_single_col_data(fns, tablefn_map, feat_map,
                                  storefun, featcolnr, colmap)
    if isobqcolpattern is not None:
        isocolmap = get_colmap(fns, isobqcolpattern, antipattern=psmnrpattern)
    else:
        return
    if psmnrpattern is not None:
        psmcolmap = get_colmap(fns, psmnrpattern)
    else:
        psmcolmap = False
    create_isobaric_quant_lookup(fns, tablefn_map,
                                 feat_map, pqdb,
                                 featcolnr,
                                 isocolmap, psmcolmap)
652901
def store_psm_protein_relations(fn, header, pgdb, proteins):
    
    # TODO do we need an OrderedDict or is regular dict enough?
    # Sorting for psm_id useful?
    allpsms = OrderedDict()
    last_id, psmids_to_store = None, set()
    store_soon = False
    for psm in tsvreader.generate_tsv_psms(fn, header):
        psm_id, prots = tsvreader.get_pepproteins(psm)
        prots = [x for x in prots if x in proteins]
        try:
            # In case the PSMs are presented unrolled
            allpsms[psm_id].extend(prots)
        except KeyError:
            allpsms[psm_id] = prots
        if len(psmids_to_store) % DB_STORE_CHUNK == 0:
            store_soon = True
        if store_soon and last_id != psm_id:
            pgdb.store_peptides_proteins(allpsms, psmids_to_store)
            store_soon = False
            psmids_to_store = set()
        psmids_to_store.add(psm_id)
        last_id = psm_id
    if len(psmids_to_store) > 0:
        pgdb.store_peptides_proteins(allpsms, psmids_to_store)
    pgdb.index_protein_peptides()
    return allpsms
1177057
function validateColumns(columns) {
  if (!columns) {
    return;
  }

  // check IDs
  var ids = {};
  columns.forEach(function (column, i) {
    var id = column.id;

    if (!ids[id]) {
      ids[id] = [i];
    } else {
      ids[id].push(i);
    }
  });
  Object.keys(ids).forEach(function (id) {
    if (ids[id].length > 1) {
      console.warn('Column ID \'' + id + '\' used in multiple columns ' + ids[id].join(', '), ids[id].map(function (index) {
        return columns[index];
      }));
    }
  });
}
1151154
function __() {
    var us = arguments.length > 0 && arguments[0] !== undefined ? arguments[0] : [];
    var store = arguments.length > 1 && arguments[1] !== undefined ? arguments[1] : false;

    store = typeof us === "boolean" ? us : store;
    us = typeof us === "boolean" ? [] : us;

    var seq = []; //seq is vanilla JS array + features
    seq.id = getID();
    seq.store = store;
    seq.us = us;
    seq.ds = [];
    seq.updatedTo = {};
    seq.us.map(function (u) {
      //  seq is a member of u.ds
      u.ds[u.ds.length] = seq;
      u.updatedTo[seq.id] = 0; // non-interference dependency
    });
    seq.IndexOnTimestamp = {};
    seq.TimestampOnIndex = {};
    seq.propagating = 0;
    seq.done = 0;
    seq.eqs = [];
    seq.addEq = function (eq) {
      return seq.eqs[seq.eqs.length] = eq;
    };
    seq.evalEqs = function (val) {
      var val0 = val;
      seq.eqs.map(function (eq) {
        return val = eq(val);
      });
      return seq.us.length === 0 ? val0 : val;
    };
    //api-----------------
    Object.keys(__.api).map(function (api) {
      seq[api] = __.api[api](__, seq, store);
    });
    //-----------------
    Object.defineProperties(seq, //detect t update on each seqs
    {
      t: {
        //foo.t
        get: function get() {
          return seq.valOnT;
        },
        set: function set(tval) {
          var sanityCheck = seq.propagating === 0 && seq.us.length !== 0 ? function () {
            throw new Error("Do not set a value of the sequence that depends on other sequences!");
          }() : function () {
            seq.propagating = 0;
            var core = seq.done === 0 ? function () {
              seq.valOnT = seq.evalEqs(tval); //self eqs eval
              seq.T = Date.now();
              var core1 = store ? function () {
                seq.IndexOnTimestamp[seq.T] = seq.length;
                seq.TimestampOnIndex[seq.length] = seq.T;
                seq[seq.length] = seq.valOnT; //after funcs
              }() : true;
              Object.keys(seq.updatedTo).map(function (key) {
                seq.updatedTo[key] = 1;
              });
              //clear updated us in non-interference way
              seq.us.map(function (u) {
                return u.updatedTo[seq.id] = 0;
              });
              seq.ds.map(function (d) {
                var usAllUpdated = d.us.map(function (u) {
                  return u.updatedTo[d.id];
                }).reduce(function (a, b) {
                  return a * b;
                });
                var core2 = usAllUpdated === 1 ? function () {
                  d.propagating = 1;
                  d.t = d.us.map(function (u) {
                    return u.t;
                  });
                }() : true;
              });
            }() : true;
          }();
        }
      }
    });
    return seq;
  }
1344376
async function createColumns(workSheetInfo, fieldNames, startIndex) {
    if (!util.isArray(fieldNames)) {
        return Promise.reject(
            new Error('Field names should be specified in array format')
        );
    }

    let promises = fieldNames.map(fieldName => {
        let payload = api.converter.createFieldRequest(fieldName, ++startIndex);
        let requestPayload = util._extend(
            {body: payload, cellId: 'R1C' + startIndex},
            workSheetInfo
        );

        return executeRequest('create_field', requestPayload);
    });

    let response = await Promise.all(promises);
    cache.clear();
    return response;
}
770085
def merge_peptides(fns, ns):
    """"""
    peptides_to_map = reader.generate_peptides_multiple_fractions(fns, ns)
    psmmap = create_merge_psm_map(peptides_to_map, ns)
    peptides = reader.generate_peptides_multiple_fractions(fns, ns)
    for peptide in peptides:
        seq = reader.get_peptide_seq(peptide, ns)
        psm_ids = reader.get_psm_ids_from_peptide(peptide, ns)
        # remove current psm ids, repopulate with stored ones
        psm_ids.clear()
        for new_psm_id in psmmap[seq]:
            etree.SubElement(psm_ids, 'psm_id').text = new_psm_id
        yield formatting.string_and_clear(peptide, ns)
868219
def plotter(molecules, ensemble_lookup, options):
    

    try:
        import matplotlib

        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
    except ImportError:
        print("\n Plotting requires matplotlib to be installed\n")
        sys.exit(1)

    for ensemble in ensemble_lookup.keys():

        # create figure
        fig = plt.figure()

        # create the queries subplot, the left subplot
        # create the left hand subplot
        ax1 = fig.add_subplot(121)

        for query in sorted(ensemble_lookup[ensemble]):
            query_list = []
            query_list.append(query)
            score_structure = classification.make_score_structure(molecules, query_list)
            auc_structure = classification.make_auc_structure(score_structure)
            tpf = []
            fpf = []
            for mol in auc_structure:
                fpf.append(mol[4])
                tpf.append(mol[5])

            # add axis-labels and a title
            ax1.set_xlabel('FPF')
            ax1.set_ylabel('TPF')
            title = 'query performance'
            ax1.set_title(title)

            # add plot data and labels for the legend
            lbl = query
            ax1.plot(fpf, tpf, lw=3, label=lbl)

        # get legend handles and labels, then reverse their order
        handles, labels = ax1.get_legend_handles_labels()
        ax1.legend(handles[::-1], labels[::-1])

        # add the legend
        ax1.legend(handles, labels, loc='best')

        # create the ensemble subplot, the right subplot
        score_structure = classification.make_score_structure(molecules, ensemble_lookup[ensemble])
        auc_structure = classification.make_auc_structure(score_structure)
        tpf = []
        fpf = []
        for mol in auc_structure:
            fpf.append(mol[4])
            tpf.append(mol[5])

        # create right hand subplot
        ax2 = fig.add_subplot(122)

        # add axis-labels and a title
        ax2.set_xlabel('FPF')
        ax2.set_ylabel('TPF')
        title = 'ensemble performance'
        ax2.set_title(title)

        # add plot data and a label for the legend
        lbl = 'ensemble'
        ax2.plot(fpf, tpf, lw=3, label=lbl)

        # get legend handles and labels, then reverse their order
        handles, labels = ax2.get_legend_handles_labels()
        ax2.legend(handles[::-1], labels[::-1])

        # add the legend
        ax2.legend(handles, labels, loc='best')

        # save figure
        figurename = options.outname + '_' + ensemble.replace('.csv', '') + '.pdf'
        filename = os.path.join(os.getcwd(), figurename)
        plt.savefig(filename, bbobx='tight', format='pdf')
889831
def get_peptable_headerfields(headertypes, lookup=False, poolnames=False):
    """"""
    field_defs = {'isoquant': get_isoquant_fields,
                  'precursorquant': get_precursorquant_fields,
                  'peptidefdr': get_peptidefdr_fields,
                  'peptidepep': get_peptidepep_fields,
                  'proteindata': get_proteininfo_fields,
                  }
    return generate_headerfields(headertypes, field_defs, poolnames, lookup)
120641
def read_dependent_peptides(filename):
    
    df = (pd.read_perseus(filename, usecols=_cols)
            .dropna(subset=['DP Ratio mod/base']))
    df['DP Ratio mod/base'] = df['DP Ratio mod/base'].astype(float)
    dep = df.pivot_table('DP Ratio mod/base', index=_index_columns,
                columns='Raw file', aggfunc=np.median)
    localization = _count_localizations(df)
    return dep, localization
547150
def map_psmnrcol_to_quantcol(quantcols, psmcols, tablefn_map):
    """"""
    if not psmcols:
        for fn in quantcols:
            for qcol in quantcols[fn]:
                yield (tablefn_map[fn], qcol)
    else:
        for fn in quantcols:
            for qcol, psmcol in zip(quantcols[fn], psmcols[fn]):
                yield (tablefn_map[fn], qcol, psmcol)
861267
def _run_qc_tools(bam_file, data):
    
    from bcbio.qc import (atropos, contamination, coverage, damage, fastqc, kraken,
                          qsignature, qualimap, samtools, picard, srna, umi, variant,
                          viral, preseq, chipseq)
    tools = {"fastqc": fastqc.run,
             "atropos": atropos.run,
             "small-rna": srna.run,
             "samtools": samtools.run,
             "qualimap": qualimap.run,
             "qualimap_rnaseq": qualimap.run_rnaseq,
             "qsignature": qsignature.run,
             "contamination": contamination.run,
             "coverage": coverage.run,
             "damage": damage.run,
             "variants": variant.run,
             "peddy": peddy.run_qc,
             "kraken": kraken.run,
             "picard": picard.run,
             "umi": umi.run,
             "viral": viral.run,
             "preseq": preseq.run,
             "chipqc": chipseq.run
             }
    qc_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "qc", data["description"]))
    metrics = {}
    qc_out = utils.deepish_copy(dd.get_summary_qc(data))
    for program_name in dd.get_algorithm_qc(data):
        if not bam_file and program_name != "kraken":  # kraken doesn't need bam
            continue
        if dd.get_phenotype(data) == "germline" and program_name != "variants":
            continue
        qc_fn = tools[program_name]
        cur_qc_dir = os.path.join(qc_dir, program_name)
        out = qc_fn(bam_file, data, cur_qc_dir)
        qc_files = None
        if out and isinstance(out, dict):
            # Check for metrics output, two cases:
            # 1. output with {"metrics"} and files ("base")
            if "metrics" in out:
                metrics.update(out.pop("metrics"))
            # 2. a dictionary of metrics
            elif "base" not in out:
                metrics.update(out)
            # Check for files only output
            if "base" in out:
                qc_files = out
        elif out and isinstance(out, six.string_types) and os.path.exists(out):
            qc_files = {"base": out, "secondary": []}
        if not qc_files:
            qc_files = _organize_qc_files(program_name, cur_qc_dir)
        if qc_files:
            qc_out[program_name] = qc_files

    metrics["Name"] = dd.get_sample_name(data)
    metrics["Quality format"] = dd.get_quality_format(data).lower()
    return {"qc": qc_out, "metrics": metrics}
737426
def generate_tsv_pep_protein_quants(fns):
    """"""
    for fn in fns:
        header = get_tsv_header(fn)
        for pquant in generate_split_tsv_lines(fn, header):
            yield os.path.basename(fn), header, pquant
189322
def create_searchspace(lookup, fastafn, proline_cut=False,
                       reverse_seqs=True, do_trypsinize=True):
    """"""
    allpeps = []
    for record in SeqIO.parse(fastafn, 'fasta'):
        if do_trypsinize:
            pepseqs = trypsinize(record.seq, proline_cut)
        else:
            pepseqs = [record.seq]
        # Exchange all leucines to isoleucines because MS can't differ
        pepseqs = [(str(pep).replace('L', 'I'),) for pep in pepseqs]
        allpeps.extend(pepseqs)
        if len(allpeps) > 1000000:  # more than x peps, write to SQLite
            lookup.write_peps(allpeps, reverse_seqs)
            allpeps = []
    # write remaining peps to sqlite
    lookup.write_peps(allpeps, reverse_seqs)
    lookup.index_peps(reverse_seqs)
    lookup.close_connection()
427197
def _label_select_column(self, select, column,
                             populate_result_map,
                             asfrom, column_clause_args,
                             name=None,
                             within_columns_clause=True):
        """"""

        if column.type._has_column_expression and \
                populate_result_map:
            col_expr = column.type.column_expression(column)
            add_to_result_map = lambda keyname, name, objects, type_: \
                self._add_to_result_map(
                    keyname, name,
                    (column,) + objects, type_)
        else:
            col_expr = column
            if populate_result_map:
                add_to_result_map = self._add_to_result_map
            else:
                add_to_result_map = None

        if not within_columns_clause:
            result_expr = col_expr
        elif isinstance(column, elements.Label):
            if col_expr is not column:
                result_expr = compiler._CompileLabel(
                    col_expr,
                    column.name,
                    alt_names=(column.element,)
                )
            else:
                result_expr = col_expr

        elif select is not None and name:
            result_expr = compiler._CompileLabel(
                col_expr,
                name,
                alt_names=(column._key_label,)
            )

        elif \
                                                asfrom and \
                                                isinstance(column,
                                                           elements.ColumnClause) and \
                                        not column.is_literal and \
                                        column.table is not None and \
                        not isinstance(column.table, selectable.Select):
            result_expr = compiler._CompileLabel(col_expr,
                                                 elements._as_truncated(
                                                     column.name),
                                                 alt_names=(column.key,))
        elif (
                        not isinstance(column, elements.TextClause) and
                        (
                                    not isinstance(column,
                                                   elements.UnaryExpression) or
                                    column.wraps_column_expression
                        ) and
                    (
                                not hasattr(column, 'name') or
                                isinstance(column, functions.Function)
                    )
        ):
            result_expr = compiler._CompileLabel(col_expr, column.anon_label)
        elif col_expr is not column:
            # TODO: are we sure "column" has a .name and .key here ?
            # assert isinstance(column, elements.ColumnClause)
            result_expr = compiler._CompileLabel(col_expr,
                                                 elements._as_truncated(
                                                     column.name),
                                                 alt_names=(column.key,))
        else:
            result_expr = col_expr

        column_clause_args.update(
            within_columns_clause=within_columns_clause,
            add_to_result_map=add_to_result_map
        )
        return result_expr._compiler_dispatch(
            self,
            **column_clause_args
        )
6574
def compare(molecules, ensemble_lookup, options):
    

    print(" Analyzing differences ... ")
    print('')
    sort_order = classification.get_sort_order(molecules)

    ensemble1 = sorted(ensemble_lookup.keys())[0]
    ensemble2 = sorted(ensemble_lookup.keys())[1]

    stats = {}
    stats['header'] = [' ']
    name = os.path.basename(ensemble1).replace('.csv', '')
    stats['header'].append(name)
    name = os.path.basename(ensemble2).replace('.csv', '')
    stats['header'].append(name)
    stats['header'].append('Difference')
    stats['header'].append('95% CI')
    stats['header'].append('p-value')

    molecules1 = copy.deepcopy(molecules)
    molecules2 = copy.deepcopy(molecules)

    score_structure1 = classification.make_score_structure(molecules1, ensemble_lookup[ensemble1])
    score_structure2 = classification.make_score_structure(molecules2, ensemble_lookup[ensemble2])

    auc_structure_1 = classification.make_auc_structure(score_structure1)
    auc_structure_2 = classification.make_auc_structure(score_structure2)

    # calculate auc value differences
    auc_diff = classification.calculate_auc_diff(auc_structure_1, auc_structure_2, sort_order)

    stats['AUC'] = auc_diff

    # calculate enrichment factor differences
    fpfList = make_fpfList(options)
    for fpf in fpfList:
        fpf = float(fpf)
        ef_structure1 = classification.make_ef_structure(score_structure1, fpf, sort_order)
        ef_structure2 = classification.make_ef_structure(score_structure2, fpf, sort_order)

        if ef_structure1 and ef_structure2:
            ef_diff = classification.calculate_ef_diff(ef_structure1, ef_structure2, fpf)
            title = 'E%s' % fpf
            stats[title] = ef_diff

    # write results summary
    output.write_diff_summary(stats, options)

    # write roc curves
    if options.write_roc:
        print(" Writing ROC data ... ")
        print('')
        output.write_roc(auc_structure_1, ensemble1, options)
        output.write_roc(auc_structure_2, ensemble2, options)

    # plot
    if options.plot:
        print(" Making plots ... ")
        print('')
        plotter(molecules, ensemble_lookup, options)
941997
function(columns, desc) {
    var self = this.chain()

    if (columns) {
      if (typeof desc === 'boolean') {
        if (!Array.isArray(columns)) columns = [columns]
      } else {
        columns = [].concat.apply(
          [],
          this.definition.store.utils.args(arguments)
        )
        desc = false
      }

      for (var i in columns) {
        self.addInternal('order', {
          column: self.store.toInternalAttributeName(columns[i]),
          order: desc ? 'DESC' : 'ASC'
        })
      }
    } else {
      self.clearInternal('order')
    }

    return self
  }
1183403
function categorizeBGDatums(opts) {
    var treatments = opts.treatments;
    // this sorts the treatments collection in order.
    treatments.sort(function (a, b) {
        var aDate = new Date(tz(a.timestamp));
        var bDate = new Date(tz(b.timestamp));
        //console.error(aDate);
        return bDate.getTime() - aDate.getTime();
    });
    var profileData = opts.profile;
    if (typeof(opts.glucose) !== 'undefined') {
        //var glucoseData = opts.glucose;
        var glucoseData = opts.glucose.map(function prepGlucose (obj) {
            //Support the NS sgv field to avoid having to convert in a custom way
            obj.glucose = obj.glucose || obj.sgv;
            return obj;
        });
    }
    // this sorts the glucose collection in order.
    glucoseData.sort(function (a, b) {
        var aDate = new Date(tz(a.date));
        var bDate = new Date(tz(b.date));
        //console.error(aDate);
        return bDate.getTime() - aDate.getTime();
    });
    // if (typeof(opts.preppedGlucose) !== 'undefined') {
        // var preppedGlucoseData = opts.preppedGlucose;
    // }
    //starting variable at 0
    var boluses = 0;
    var maxCarbs = 0;
    //console.error(treatments);
    if (!treatments) return {};

    //console.error(glucoseData);
    var IOBInputs = {
        profile: profileData
    ,   history: opts.pumpHistory
    };
    var CSFGlucoseData = [];
    var ISFGlucoseData = [];
    var basalGlucoseData = [];
    var UAMGlucoseData = [];
    var CRData = [];

    var bucketedData = [];
    bucketedData[0] = glucoseData[0];
    j=0;
    //for loop to validate and bucket the data
    for (var i=1; i < glucoseData.length; ++i) {
        var BGTime;
        var lastBGTime;
        if (glucoseData[i].date) {
            BGTime = new Date(glucoseData[i].date);
        } else if (glucoseData[i].displayTime) {
            BGTime = new Date(glucoseData[i].displayTime.replace('T', ' '));
        } else if (glucoseData[i].dateString) {
            BGTime = new Date(glucoseData[i].dateString);
        } else { console.error("Could not determine BG time"); }
        if (glucoseData[i-1].date) {
            lastBGTime = new Date(glucoseData[i-1].date);
        } else if (glucoseData[i-1].displayTime) {
            lastBGTime = new Date(glucoseData[i-1].displayTime.replace('T', ' '));
        } else if (glucoseData[i-1].dateString) {
            lastBGTime = new Date(glucoseData[i-1].dateString);
        } else { console.error("Could not determine last BG time"); }
        if (glucoseData[i].glucose < 39 || glucoseData[i-1].glucose < 39) {
            continue;
        }
        var elapsedMinutes = (BGTime - lastBGTime)/(60*1000);
        if(Math.abs(elapsedMinutes) > 2) {
            j++;
            bucketedData[j]=glucoseData[i];
            bucketedData[j].date = BGTime.getTime();
            if (! bucketedData[j].dateString) {
                bucketedData[j].dateString = BGTime.toISOString();
            }
        } else {
            // if duplicate, average the two
            bucketedData[j].glucose = (bucketedData[j].glucose + glucoseData[i].glucose)/2;
        }
    }
    //console.error(bucketedData);
    //console.error(bucketedData[bucketedData.length-1]);
    // go through the treatments and remove any that are older than the oldest glucose value
    //console.error(treatments);
    for (var i=treatments.length-1; i>0; --i) {
        var treatment = treatments[i];
        //console.error(treatment);
        if (treatment) {
            var treatmentDate = new Date(tz(treatment.timestamp));
            var treatmentTime = treatmentDate.getTime();
            var glucoseDatum = bucketedData[bucketedData.length-1];
            //console.error(glucoseDatum);
            var BGDate = new Date(glucoseDatum.date);
            var BGTime = BGDate.getTime();
            if ( treatmentTime < BGTime ) {
                treatments.splice(i,1);
            }
        }
    }
    //console.error(treatments);
    var calculatingCR = false;
    var absorbing = 0;
    var uam = 0; // unannounced meal
    var mealCOB = 0;
    var mealCarbs = 0;
    var CRCarbs = 0;
    var type="";
    // main for loop
    var fullHistory = IOBInputs.history;
    for (var i=bucketedData.length-5; i > 0; --i) {
        var glucoseDatum = bucketedData[i];
        //console.error(glucoseDatum);
        var BGDate = new Date(glucoseDatum.date);
        var BGTime = BGDate.getTime();
        // As we're processing each data point, go through the treatment.carbs and see if any of them are older than
        // the current BG data point.  If so, add those carbs to COB.
        var treatment = treatments[treatments.length-1];
        var myCarbs = 0;
        if (treatment) {
            var treatmentDate = new Date(tz(treatment.timestamp));
            var treatmentTime = treatmentDate.getTime();
            //console.error(treatmentDate);
            if ( treatmentTime < BGTime ) {
                if (treatment.carbs >= 1) {
                    mealCOB += parseFloat(treatment.carbs);
                    mealCarbs += parseFloat(treatment.carbs);
                    myCarbs = treatment.carbs;
                }
                treatments.pop();
            }
        }

        var BG;
        var avgDelta;
        var delta;
        // TODO: re-implement interpolation to avoid issues here with gaps
        // calculate avgDelta as last 4 datapoints to better catch more rises after COB hits zero
        if (typeof(bucketedData[i].glucose) != 'undefined' && typeof(bucketedData[i+4].glucose) != 'undefined') {
            //console.error(bucketedData[i]);
            BG = bucketedData[i].glucose;
            if ( BG < 40 || bucketedData[i+4].glucose < 40) {
                //process.stderr.write("!");
                continue;
            }
            avgDelta = (BG - bucketedData[i+4].glucose)/4;
            delta = (BG - bucketedData[i+1].glucose);
        } else { console.error("Could not find glucose data"); }

        avgDelta = avgDelta.toFixed(2);
        glucoseDatum.avgDelta = avgDelta;

        //sens = ISF
        var sens = ISF.isfLookup(IOBInputs.profile.isfProfile,BGDate);
        IOBInputs.clock=BGDate.toISOString();
        // trim down IOBInputs.history to just the data for 6h prior to BGDate
        //console.error(IOBInputs.history[0].created_at);
        var newHistory = [];
        for (h=0; h<fullHistory.length; h++) {
            var hDate = new Date(fullHistory[h].created_at)
            //console.error(fullHistory[i].created_at, hDate, BGDate, BGDate-hDate);
            //if (h == 0 || h == fullHistory.length - 1) {
                //console.error(hDate, BGDate, hDate-BGDate)
            //}
            if (BGDate-hDate < 6*60*60*1000 && BGDate-hDate > 0) {
                //process.stderr.write("i");
                //console.error(hDate);
                newHistory.push(fullHistory[h]);
            }
        }
        IOBInputs.history = newHistory;
        // process.stderr.write("" + newHistory.length + " ");
        //console.error(newHistory[0].created_at,newHistory[newHistory.length-1].created_at,newHistory.length);


        // for IOB calculations, use the average of the last 4 hours' basals to help convergence;
        // this helps since the basal this hour could be different from previous, especially if with autotune they start to diverge.
        // use the pumpbasalprofile to properly calculate IOB during periods where no temp basal is set
        currentPumpBasal = basal.basalLookup(opts.pumpbasalprofile, BGDate);
        BGDate1hAgo = new Date(BGTime-1*60*60*1000);
        BGDate2hAgo = new Date(BGTime-2*60*60*1000);
        BGDate3hAgo = new Date(BGTime-3*60*60*1000);
        basal1hAgo = basal.basalLookup(opts.pumpbasalprofile, BGDate1hAgo);
        basal2hAgo = basal.basalLookup(opts.pumpbasalprofile, BGDate2hAgo);
        basal3hAgo = basal.basalLookup(opts.pumpbasalprofile, BGDate3hAgo);
        var sum = [currentPumpBasal,basal1hAgo,basal2hAgo,basal3hAgo].reduce(function(a, b) { return a + b; });
        IOBInputs.profile.currentBasal = Math.round((sum/4)*1000)/1000;

        // this is the current autotuned basal, used for everything else besides IOB calculations
        currentBasal = basal.basalLookup(opts.basalprofile, BGDate);

        //console.error(currentBasal,basal1hAgo,basal2hAgo,basal3hAgo,IOBInputs.profile.currentBasal);
        // basalBGI is BGI of basal insulin activity.
        basalBGI = Math.round(( currentBasal * sens / 60 * 5 )*100)/100; // U/hr * mg/dL/U * 1 hr / 60 minutes * 5 = mg/dL/5m
        //console.log(JSON.stringify(IOBInputs.profile));
        // call iob since calculated elsewhere
        var iob = getIOB(IOBInputs)[0];
        //console.error(JSON.stringify(iob));

        // activity times ISF times 5 minutes is BGI
        var BGI = Math.round(( -iob.activity * sens * 5 )*100)/100;
        // datum = one glucose data point (being prepped to store in output)
        glucoseDatum.BGI = BGI;
        // calculating deviation
        deviation = avgDelta-BGI;

        // set positive deviations to zero if BG is below 80
        if ( BG < 80 && deviation > 0 ) {
            deviation = 0;
        }

        // rounding and storing deviation
        deviation = deviation.toFixed(2);
        glucoseDatum.deviation = deviation;


        // Then, calculate carb absorption for that 5m interval using the deviation.
        if ( mealCOB > 0 ) {
            var profile = profileData;
            ci = Math.max(deviation, profile.min_5m_carbimpact);
            absorbed = ci * profile.carb_ratio / sens;
            // Store the COB, and use it as the starting point for the next data point.
            mealCOB = Math.max(0, mealCOB-absorbed);
        }


        // Calculate carb ratio (CR) independently of CSF and ISF
        // Use the time period from meal bolus/carbs until COB is zero and IOB is < currentBasal/2
        // For now, if another meal IOB/COB stacks on top of it, consider them together
        // Compare beginning and ending BGs, and calculate how much more/less insulin is needed to neutralize
        // Use entered carbs vs. starting IOB + delivered insulin + needed-at-end insulin to directly calculate CR.

        if (mealCOB > 0 || calculatingCR ) {
            // set initial values when we first see COB
            CRCarbs += myCarbs;
            if (!calculatingCR) {
                CRInitialIOB = iob.iob;
                CRInitialBG = glucoseDatum.glucose;
                CRInitialCarbTime = new Date(glucoseDatum.date);
                console.error("CRInitialIOB:",CRInitialIOB,"CRInitialBG:",CRInitialBG,"CRInitialCarbTime:",CRInitialCarbTime);
            }
            // keep calculatingCR as long as we have COB or enough IOB
            if ( mealCOB > 0 && i>1 ) {
                calculatingCR = true;
            } else if ( iob.iob > currentBasal/2 && i>1 ) {
                calculatingCR = true;
            // when COB=0 and IOB drops low enough, record end values and be done calculatingCR
            } else {
                CREndIOB = iob.iob;
                CREndBG = glucoseDatum.glucose;
                CREndTime = new Date(glucoseDatum.date);
                console.error("CREndIOB:",CREndIOB,"CREndBG:",CREndBG,"CREndTime:",CREndTime);
                var CRDatum = {
                    CRInitialIOB: CRInitialIOB
                ,   CRInitialBG: CRInitialBG
                ,   CRInitialCarbTime: CRInitialCarbTime
                ,   CREndIOB: CREndIOB
                ,   CREndBG: CREndBG
                ,   CREndTime: CREndTime
                ,   CRCarbs: CRCarbs
                };
                //console.error(CRDatum);

                var CRElapsedMinutes = Math.round((CREndTime - CRInitialCarbTime) / 1000 / 60);
                //console.error(CREndTime - CRInitialCarbTime, CRElapsedMinutes);
                if ( CRElapsedMinutes < 60 || ( i==1 && mealCOB > 0 ) ) {
                    console.error("Ignoring",CRElapsedMinutes,"m CR period.");
                } else {
                    CRData.push(CRDatum);
                }

                CRCarbs = 0;
                calculatingCR = false;
            }
        }


        // If mealCOB is zero but all deviations since hitting COB=0 are positive, assign those data points to CSFGlucoseData
        // Once deviations go negative for at least one data point after COB=0, we can use the rest of the data to tune ISF or basals
        if (mealCOB > 0 || absorbing || mealCarbs > 0) {
            // if meal IOB has decayed, then end absorption after this data point unless COB > 0
            if ( iob.iob < currentBasal/2 ) {
                absorbing = 0;
            // otherwise, as long as deviations are positive, keep tracking carb deviations
            } else if (deviation > 0) {
                absorbing = 1;
            } else {
                absorbing = 0;
            }
            if ( ! absorbing && ! mealCOB ) {
                mealCarbs = 0;
            }
            // check previous "type" value, and if it wasn't csf, set a mealAbsorption start flag
            //console.error(type);
            if ( type != "csf" ) {
                glucoseDatum.mealAbsorption = "start";
                console.error(glucoseDatum.mealAbsorption,"carb absorption");
            }
            type="csf";
            glucoseDatum.mealCarbs = mealCarbs;
            //if (i == 0) { glucoseDatum.mealAbsorption = "end"; }
            CSFGlucoseData.push(glucoseDatum);
        } else {
          // check previous "type" value, and if it was csf, set a mealAbsorption end flag
          if ( type === "csf" ) {
            CSFGlucoseData[CSFGlucoseData.length-1].mealAbsorption = "end";
            console.error(CSFGlucoseData[CSFGlucoseData.length-1].mealAbsorption,"carb absorption");
          }

          if ((iob.iob > currentBasal || deviation > 6 || uam) ) {
            if (deviation > 0) {
                uam = 1;
            } else {
                uam = 0;
            }
            if ( type != "uam" ) {
                glucoseDatum.uamAbsorption = "start";
                console.error(glucoseDatum.uamAbsorption,"uannnounced meal absorption");
            }
            type="uam";
            UAMGlucoseData.push(glucoseDatum);
          } else {
            if ( type === "uam" ) {
                console.error("end unannounced meal absorption");
            }


            // Go through the remaining time periods and divide them into periods where scheduled basal insulin activity dominates. This would be determined by calculating the BG impact of scheduled basal insulin (for example 1U/hr * 48 mg/dL/U ISF = 48 mg/dL/hr = 5 mg/dL/5m), and comparing that to BGI from bolus and net basal insulin activity.
            // When BGI is positive (insulin activity is negative), we want to use that data to tune basals
            // When BGI is smaller than about 1/4 of basalBGI, we want to use that data to tune basals
            // When BGI is negative and more than about 1/4 of basalBGI, we can use that data to tune ISF,
            // unless avgDelta is positive: then that's some sort of unexplained rise we don't want to use for ISF, so that means basals
            if (basalBGI > -4 * BGI) {
                type="basal";
                basalGlucoseData.push(glucoseDatum);
            } else {
                if ( avgDelta > 0 && avgDelta > -2*BGI ) {
                    //type="unknown"
                    type="basal"
                    basalGlucoseData.push(glucoseDatum);
                } else {
                    type="ISF";
                    ISFGlucoseData.push(glucoseDatum);
                }
            }
          }
        }
        // debug line to print out all the things
        BGDateArray = BGDate.toString().split(" ");
        BGTime = BGDateArray[4];
        console.error(absorbing.toString(),"mealCOB:",mealCOB.toFixed(1),"mealCarbs:",mealCarbs,"basalBGI:",basalBGI.toFixed(1),"BGI:",BGI.toFixed(1),"IOB:",iob.iob.toFixed(1),"at",BGTime,"dev:",deviation,"avgDelta:",avgDelta,type);
    }

    var IOBInputs = {
        profile: profileData
    ,   history: opts.pumpHistory
    };
    var treatments = find_insulin(IOBInputs);
    CRData.forEach(function(CRDatum) {
        var dosedOpts = {
            treatments: treatments
            , profile: opts.profile
            , start: CRDatum.CRInitialCarbTime
            , end: CRDatum.CREndTime
        };
        var insulinDosed = dosed(dosedOpts);
        CRDatum.CRInsulin = insulinDosed.insulin;
        //console.error(CRDatum);
    });

    var CSFLength = CSFGlucoseData.length;
    var ISFLength = ISFGlucoseData.length;
    var UAMLength = UAMGlucoseData.length;
    var basalLength = basalGlucoseData.length;

    if (opts.categorize_uam_as_basal) {
        console.error("Categorizing all UAM data as basal.");
        basalGlucoseData = basalGlucoseData.concat(UAMGlucoseData);
    } else if (2*basalLength < UAMLength) {
        //console.error(basalGlucoseData, UAMGlucoseData);
        console.error("Warning: too many deviations categorized as UnAnnounced Meals");
        console.error("Adding",UAMLength,"UAM deviations to",basalLength,"basal ones");
        basalGlucoseData = basalGlucoseData.concat(UAMGlucoseData);
        //console.error(basalGlucoseData);
        // if too much data is excluded as UAM, add in the UAM deviations, but then discard the highest 50%
        basalGlucoseData.sort(function (a, b) {
            return a.deviation - b.deviation;
        });
        var newBasalGlucose = basalGlucoseData.slice(0,basalGlucoseData.length/2);
        //console.error(newBasalGlucose);
        basalGlucoseData = newBasalGlucose;
        console.error("and selecting the lowest 50%, leaving", basalGlucoseData.length, "basal+UAM ones");

        console.error("Adding",UAMLength,"UAM deviations to",ISFLength,"ISF ones");
        ISFGlucoseData = ISFGlucoseData.concat(UAMGlucoseData);
        //console.error(ISFGlucoseData.length, UAMLength);
    }
    var basalLength = basalGlucoseData.length;
    var ISFLength = ISFGlucoseData.length;
    if ( 4*basalLength + ISFLength < CSFLength && ISFLength < 10 ) {
        console.error("Warning: too many deviations categorized as meals");
        //console.error("Adding",CSFLength,"CSF deviations to",basalLength,"basal ones");
        //var basalGlucoseData = basalGlucoseData.concat(CSFGlucoseData);
        console.error("Adding",CSFLength,"CSF deviations to",ISFLength,"ISF ones");
        ISFGlucoseData = ISFGlucoseData.concat(CSFGlucoseData);
        CSFGlucoseData = [];
    }


    return {
        CRData: CRData,
        CSFGlucoseData: CSFGlucoseData,
        ISFGlucoseData: ISFGlucoseData,
        basalGlucoseData: basalGlucoseData
    };
}