C++ TimeCompression::npts方法代码示例

void CumulativeLostHydrogens ( float **model_trace, float **incorporation_rise,
                               TimeCompression &time_c, RegionTracker &my_regions,
                               bool use_vectorization,
                               float *vec_tau_top )
{
  // Now solve the top of the well cumulative lost hydrogen ions
  // happen faster, hence tau_top
#ifdef __INTEL_COMPILER
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
    {
      RedSolveHydrogenFlowInWell ( model_trace[fnum],incorporation_rise[fnum],time_c.npts(),my_regions.cache_step.i_start_coarse_step[fnum],&time_c.deltaFrame[0],vec_tau_top[fnum] ); // we lose hydrogen ions fast!

    }
  }
#else
  if ( use_vectorization )
  {
    RedSolveHydrogenFlowInWell_Vec ( NUMFB,model_trace,incorporation_rise,time_c.npts(),&time_c.deltaFrame[0],vec_tau_top ); // we lose hydrogen ions fast!

  }
  else
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
    {
      RedSolveHydrogenFlowInWell ( model_trace[fnum],incorporation_rise[fnum],time_c.npts(),my_regions.cache_step.i_start_coarse_step[fnum],&time_c.deltaFrame[0],vec_tau_top[fnum] ); // we lose hydrogen ions fast!

    }
  }
#endif

  for ( int fnum=0; fnum<NUMFB; fnum++ )
    DiminishVector ( incorporation_rise[fnum],model_trace[fnum],time_c.npts() );  // cumulative lost hydrogen ions instead of retained hydrogen ions

}

void single_flow_optimizer::AllocLevMar (TimeCompression &time_c, PoissonCDFApproxMemo *_math_poiss, float damp_kmult, bool _var_kmult_only, float kmult_low_limit, float kmult_hi_limit, float AmplLowerLimit)
{
  //@TODO:  All these fitters share the curried calc_trace object, so we can do something useful by sharing it
  math_poiss = _math_poiss;
  // one parameter fit
  oneFlowFit = new BkgModSingleFlowFit (time_c.npts(),&time_c.frameNumber[0],&time_c.deltaFrame[0],&time_c.deltaFrameSeconds[0],math_poiss,1);
  // two parameter fit
  oneFlowFitKrate = new BkgModSingleFlowFit (time_c.npts(),&time_c.frameNumber[0],&time_c.deltaFrame[0],&time_c.deltaFrameSeconds[0], math_poiss, 2);

  // set up levmar state before going off on a wild adventure
  oneFlowFit->SetLambdaThreshold (10.0);
  oneFlowFitKrate->SetLambdaThreshold (1.0);

  SetLowerLimitAmplFit (AmplLowerLimit,kmult_low_limit);
  SetUpperLimitAmplFit (LAST_POISSON_TABLE_COL,kmult_hi_limit);

  float my_prior[2] = {0.0, 1.0}; // amplitude, kmult
  float my_damper[2] = {0.0,damp_kmult*time_c.npts()}; // must make a difference to change kmult
  oneFlowFitKrate->SetPrior (&my_prior[0]);
  oneFlowFitKrate->SetDampers (&my_damper[0]);
  var_kmult_only = _var_kmult_only;

  // boot up the new projection fit >which is not an instance of LevMar<
  ProjectionFit = new ProjectionSearchOneFlow (time_c.npts(), &time_c.deltaFrame[0], &time_c.deltaFrameSeconds[0], math_poiss);
  ProjectionFit->max_paramA = LAST_POISSON_TABLE_COL;
  ProjectionFit->min_paramA = AmplLowerLimit;

  AltFit = new AlternatingDirectionOneFlow (time_c.npts(), &time_c.deltaFrame[0], &time_c.deltaFrameSeconds[0], math_poiss);
}

void MathModel::AccumulateSingleNeighborExcessHydrogenOneParameter ( float *my_xtflux, float *neighbor_signal,
    BeadParams *p, reg_params *reg_p,
    BeadScratchSpace &my_scratch, 
    buffer_params_block_flows &my_cur_buffer_block,
    TimeCompression &time_c,
    RegionTracker &my_regions, const FlowBufferInfo & my_flow,
    bool use_vectorization,
    float multiplier, bool rescale_flag,
    int flow_block_size, int flow_block_start )
{

  FillBufferParamsBlockFlows ( &my_cur_buffer_block,p,reg_p,my_flow.flow_ndx_map, flow_block_start,
    flow_block_size );

  float block_model_trace[my_scratch.bead_flow_t], block_incorporation_rise[my_scratch.bead_flow_t];

  // "In parallel, across flows"

  float* model_trace[flow_block_size];
  float* scratch_trace[flow_block_size];
  float* incorporation_rise[flow_block_size];
  float* lost_hydrogens[flow_block_size];
  float* neighbor_local[flow_block_size];

  // should this be using cur_buffer_block as usual?
  float vec_tau_well[flow_block_size];
  float vec_tau_empty[flow_block_size];

  for ( int fnum=0; fnum<flow_block_size; fnum++ )
  {
    vec_tau_well[fnum] = my_cur_buffer_block.tauB[fnum];
    vec_tau_empty[fnum] = my_cur_buffer_block.etbR[fnum]*vec_tau_well[fnum];
  }

  // setup parallel pointers into the structure
  for ( int fnum=0; fnum<flow_block_size; fnum++ )
  {

    scratch_trace[fnum] = model_trace[fnum] = &block_model_trace[fnum*time_c.npts() ];
    lost_hydrogens[fnum] = incorporation_rise[fnum] = &block_incorporation_rise[fnum*time_c.npts() ];   // set up each flow information
    neighbor_local[fnum] = &neighbor_signal[fnum*time_c.npts() ];
  }

  IncorporationRiseFromNeighborSignal ( incorporation_rise,neighbor_local, time_c,my_regions, my_scratch, my_cur_buffer_block, flow_block_size );
  // uses a scratch buffer here [recycles model_trace as we don't need it yet], 
  // turns incorporation_rise into lost_hydrogens
  CumulativeLostHydrogens ( incorporation_rise, scratch_trace, time_c, my_regions, use_vectorization, vec_tau_well, flow_block_size );
  // lost_hydrogens=incorporation_rise returned as lost_hydrogens above
  // now we generate the real model_trace we're accumulating
  IonsFromBulk ( model_trace,lost_hydrogens, time_c, my_regions, my_flow, use_vectorization, vec_tau_empty, flow_block_size );
  
  if (rescale_flag)
    RescaleTraceByBuffering(model_trace, vec_tau_well, vec_tau_empty,time_c.npts(), flow_block_size);

  // universal
  MultiplyVectorByScalar ( block_model_trace,multiplier,my_scratch.bead_flow_t ); // scale down the quantity of ions

  // add to the bulk cross-talk we're creating
  AccumulateVector ( my_xtflux, block_model_trace,my_scratch.bead_flow_t );
}

void IonsFromBulk ( float **model_trace, float **incorporation_rise,
                    TimeCompression &time_c, RegionTracker &my_regions, flow_buffer_info my_flow,
                    bool use_vectorization,
                    float *vec_tau_bulk )
{
  // finally solve the way hydrogen ions diffuse out of the bulk
#ifdef __INTEL_COMPILER
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
    {
      // Now solve the bulk
      RedSolveHydrogenFlowInWell ( model_trace[fnum],incorporation_rise[fnum],time_c.npts(),
                                   my_regions.cache_step.i_start_coarse_step[my_flow.flow_ndx_map[fnum]],&time_c.deltaFrame[0],vec_tau_bulk[fnum] ); // we retain hydrogen ions variably in the bulk depending on direction
    }
  }
#else
  if ( use_vectorization )
  {
    // Now solve the bulk
    RedSolveHydrogenFlowInWell_Vec ( NUMFB,model_trace,incorporation_rise,time_c.npts(),&time_c.deltaFrame[0],vec_tau_bulk );

  }
  else
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
    {
      // Now solve the bulk
      RedSolveHydrogenFlowInWell ( model_trace[fnum],incorporation_rise[fnum],time_c.npts(),
                                   my_regions.cache_step.i_start_coarse_step[fnum],&time_c.deltaFrame[0],vec_tau_bulk[fnum] ); // we retain hydrogen ions variably in the bulk depending on direction

    }
  }
#endif
}

void EmptyTrace::ShiftMe (
    float tshift, const TimeCompression &time_cp, const float *my_buff, float *out_buff,
    int flow_block_size ) const
{
  for (int fnum=0;fnum<flow_block_size;fnum++)
    {
      float *fbkg = out_buff + fnum*time_cp.npts();
      const float *bg = &my_buff[fnum*imgFrames];         // get ptr to start of neighbor background
      memset (fbkg,0,sizeof (float[time_cp.npts()])); // on general principles
      //fprintf(stdout, "tshift %f\n", tshift);
      
      for (int i=0;i < time_cp.npts();i++){
        // get the frame number of this data point (might be fractional because this point could be
        // the average of several frames of data.  This number is the average time of all the averaged
        // data points
        float t=time_cp.frameNumber[i];
        float fn=t-tshift;
        if (fn < 0.0f) fn = 0.0f;
        if (fn > (imgFrames-2)) fn = imgFrames-2;
        int ifn= (int) fn;
          float frac = fn - ifn;
          fbkg[i] = ( (1-frac) *bg[ifn] + frac*bg[ifn+1]);
          assert ( !isnan(fbkg[i]) );
      }
    }
}

// Multiflow: Solve incorporation
// Solve lost hydrogens to bulk
// solve bulk resistance to lost hydrogens
// this function too large, should be componentized
void AccumulateSingleNeighborXtalkTrace ( float *my_xtflux, bead_params *p, reg_params *reg_p,
    BeadScratchSpace &my_scratch, TimeCompression &time_c, RegionTracker &my_regions, flow_buffer_info my_flow,
    PoissonCDFApproxMemo *math_poiss, bool use_vectorization,
    float tau_top, float tau_bulk, float multiplier )
{

  // Compute the hydrogen signal of xtalk in the bulk fluid above the affected well
  // 1) Xtalk happens fast -> we compute the cumulative lost hydrogen ions at the top of the bead instead of at the bottom (tau_top)
  // 1a) xtalk happens as the well loses hydrogen ions so the cumulative lost = cumulative total generated - number in well currently
  // 2) hydrogen ions are "retained" by the bulk for a finite amount of time as they are swept along
  // 2a) this 'retainer' may be asymmetric (tau_bulk) models the decay rate
  // 3) multiplier: hydrogen ions spread out, so are of different proportion to the original signal

  // this is pre-calculated outside for the current region parameters
  //my_regions.cache_step.CalculateNucRiseCoarseStep(reg_p,time_c);


  // over-ride buffering parameters for bead
  // use the same incorporation parameters, though
  FillIncorporationParamsBlockFlows ( &my_scratch.cur_bead_block, p,reg_p,my_flow.flow_ndx_map,my_flow.buff_flow );
  //params_IncrementHits(p);

  float block_model_trace[my_scratch.bead_flow_t], block_incorporation_rise[my_scratch.bead_flow_t];

  // "In parallel, across flows"
  float* nuc_rise_ptr[NUMFB];
  float* model_trace[NUMFB];
  float* incorporation_rise[NUMFB];
  // should this be using cur_buffer_block as usual?
  float vec_tau_top[NUMFB];
  float vec_tau_bulk[NUMFB];

  for ( int fnum=0; fnum<NUMFB; fnum++ )
  {
    vec_tau_top[fnum] = tau_top;
    vec_tau_bulk[fnum] = tau_bulk;
  }

  // setup parallel pointers into the structure
  for ( int fnum=0; fnum<NUMFB; fnum++ )
  {
    nuc_rise_ptr[fnum] = my_regions.cache_step.NucCoarseStep ( fnum );
    model_trace[fnum] = &block_model_trace[fnum*time_c.npts() ];
    incorporation_rise[fnum] = &block_incorporation_rise[fnum*time_c.npts() ];   // set up each flow information
  }

  IncorporationRiseFromNeighborParameters ( incorporation_rise, nuc_rise_ptr, p, time_c, my_regions, my_scratch, math_poiss );

  CumulativeLostHydrogens ( model_trace,incorporation_rise, time_c, my_regions, use_vectorization, vec_tau_top );

  IonsFromBulk ( model_trace,incorporation_rise, time_c, my_regions, my_flow, use_vectorization, vec_tau_bulk );

  // universal
  MultiplyVectorByScalar ( block_model_trace,multiplier,my_scratch.bead_flow_t ); // scale down the quantity of ions

  // add to the bulk cross-talk we're creating
  AccumulateVector ( my_xtflux, block_model_trace,my_scratch.bead_flow_t );
}

void MathModel::AccumulateSingleNeighborExcessHydrogen ( float *my_xtflux, float *neighbor_signal, BeadParams *p, reg_params *reg_p,
    BeadScratchSpace &my_scratch, 
    buffer_params_block_flows &my_cur_buffer_block,
    TimeCompression &time_c,
    RegionTracker &my_regions, const FlowBufferInfo & my_flow,
    bool use_vectorization,
    float tau_top, float tau_bulk, float multiplier,
    int flow_block_size, int flow_block_start )
{
  // over-ride buffering parameters for bead

  FillBufferParamsBlockFlows ( &my_cur_buffer_block,p,reg_p,my_flow.flow_ndx_map,flow_block_start,
    flow_block_size );

  float block_model_trace[my_scratch.bead_flow_t], block_incorporation_rise[my_scratch.bead_flow_t];

  // "In parallel, across flows"

  float* model_trace[flow_block_size];
  float* scratch_trace[flow_block_size];
  float* incorporation_rise[flow_block_size];
  float* lost_hydrogens[flow_block_size];
  float* neighbor_local[flow_block_size];
  // should this be using cur_buffer_block as usual?
  float vec_tau_top[flow_block_size];
  float vec_tau_bulk[flow_block_size];

  for ( int fnum=0; fnum<flow_block_size; fnum++ )
  {
    vec_tau_top[fnum] = tau_top;
    vec_tau_bulk[fnum] = tau_bulk;
  }

  // setup parallel pointers into the structure
  for ( int fnum=0; fnum<flow_block_size; fnum++ )
  {

    scratch_trace[fnum]=model_trace[fnum] = &block_model_trace[fnum*time_c.npts() ];
    lost_hydrogens[fnum]=incorporation_rise[fnum] = &block_incorporation_rise[fnum*time_c.npts() ];   // set up each flow information
    neighbor_local[fnum] = &neighbor_signal[fnum*time_c.npts() ];
  }
  // make incorporation_rise
  IncorporationRiseFromNeighborSignal ( incorporation_rise,neighbor_local, time_c,my_regions, my_scratch, my_cur_buffer_block, flow_block_size );
  // use scratch_trace to hold temporary trace - same memory as model_trace because we don't need it yet
  // turn incorporation_rise into lost_hydrogens
  CumulativeLostHydrogens ( incorporation_rise, scratch_trace,time_c, my_regions, use_vectorization, vec_tau_top, flow_block_size );
  // lost_hydrogens = incorporation_rise 
  // now get model_trace for real in cross-talk and overwrite any temporary uses of that space
  IonsFromBulk ( model_trace,lost_hydrogens, time_c, my_regions, my_flow, use_vectorization, vec_tau_bulk, flow_block_size );
  

  // universal
  MultiplyVectorByScalar ( block_model_trace,multiplier,my_scratch.bead_flow_t ); // scale down the quantity of ions

  // add to the bulk cross-talk we're creating
  AccumulateVector ( my_xtflux, block_model_trace,my_scratch.bead_flow_t );
}

// 2nd-order background function with non-uniform bead well
void MultiFlowComputeTraceGivenIncorporationAndBackground ( float *fval,struct bead_params *p,struct reg_params *reg_p, float *ival, float *sbg,
    RegionTracker &my_regions, buffer_params_block_flows &cur_buffer_block,
    TimeCompression &time_c, flow_buffer_info &my_flow, bool use_vectorization, int bead_flow_t )
{
  float *vb_out[NUMFB];
  float *bkg_for_flow[NUMFB];
  float *new_hydrogen_for_flow[NUMFB];

  //@TODO: the natural place for vectorization is here at the flow level
  // flows are logically independent: apply "compute trace" to all flows
  // this makes for an obvious fit to solving 4 at once using the processor vectorization routines

  // parallel fill one bead parameter for block of flows
  FillBufferParamsBlockFlows ( &cur_buffer_block,p,reg_p,my_flow.flow_ndx_map,my_flow.buff_flow );

  // parallel compute across flows
  for ( int fnum=0;fnum<NUMFB;fnum++ )
  {
    vb_out[fnum] = fval + fnum*time_c.npts();        // get ptr to start of the function evaluation for the current flow
    bkg_for_flow[fnum] = &sbg[fnum*time_c.npts() ];            // get ptr to pre-shifted background
    new_hydrogen_for_flow[fnum] = &ival[fnum*time_c.npts() ];
  }

  // do the actual computation
#ifdef __INTEL_COMPILER
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
      PurpleSolveTotalTrace ( vb_out[fnum],bkg_for_flow[fnum], new_hydrogen_for_flow[fnum], time_c.npts(),
                              &time_c.deltaFrame[0], cur_buffer_block.tauB[fnum], cur_buffer_block.etbR[fnum] );
  }
#else // assumed to be GCC
  if ( use_vectorization )
  {
    PurpleSolveTotalTrace_Vec ( NUMFB, vb_out, bkg_for_flow, new_hydrogen_for_flow, time_c.npts(),
                                &time_c.deltaFrame[0], cur_buffer_block.tauB, cur_buffer_block.etbR, p->gain );
  }
  else
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
      PurpleSolveTotalTrace ( vb_out[fnum],bkg_for_flow[fnum], new_hydrogen_for_flow[fnum],time_c.npts(),
                              &time_c.deltaFrame[0], cur_buffer_block.tauB[fnum], cur_buffer_block.etbR[fnum] );
  }
#endif

  // adjust for well sensitivity, unexplained systematic effects
  // gain naturally parallel across flows
  MultiplyVectorByScalar ( fval,p->gain,bead_flow_t );

  // Dark Matter is extra background term of unexplained origin
  // Possibly should be applied directly to the observed signal rather than synthesized here inside a loop.
  if (my_regions.missing_mass.mytype == PerNucAverage)
     ApplyDarkMatter ( fval,reg_p,my_regions.missing_mass.dark_matter_compensator,my_flow.flow_ndx_map,time_c.npts() );
  else
     ApplyPCADarkMatter ( fval,p,my_regions.missing_mass.dark_matter_compensator,time_c.npts() );
}

void MultiCorrectBeadBkg ( float *block_signal_corrected, bead_params *p,
                           BeadScratchSpace &my_scratch, flow_buffer_info &my_flow, TimeCompression &time_c, RegionTracker &my_regions, float *sbg, bool use_vectorization )
{
  float vb[my_scratch.bead_flow_t];
  float* vb_out[my_flow.numfb];
  float* sbgPtr[my_flow.numfb];
  float block_bkg_plus_xtalk[my_scratch.bead_flow_t]; // set up instead of shifted background
  memset ( vb,0,sizeof ( float[my_scratch.bead_flow_t] ) );

  // add cross-talk for this bead to the empty-trace
  CopyVector ( block_bkg_plus_xtalk,sbg,my_scratch.bead_flow_t );
  AccumulateVector ( block_bkg_plus_xtalk,my_scratch.cur_xtflux_block,my_scratch.bead_flow_t );

  // compute the zeromer
  // setup pointers into the arrays
  for ( int fnum=0; fnum<my_flow.numfb; fnum++ )
  {
    // remove zeromer background - just like oneFlowFit.
    // should include xtalk (0) so I can reuse this routine
    sbgPtr[fnum] = &block_bkg_plus_xtalk[fnum*time_c.npts() ];
    vb_out[fnum] = &vb[fnum*time_c.npts() ];
  }

  // do the actual calculation in parallel or not
#ifdef __INTEL_COMPILER
  {
    for ( int fnum=0; fnum<my_flow.numfb; fnum++ )
      BlueSolveBackgroundTrace ( vb_out[fnum],sbgPtr[fnum],time_c.npts(),&time_c.deltaFrame[0],
                                 my_scratch.cur_buffer_block.tauB[fnum],my_scratch.cur_buffer_block.etbR[fnum] );
  }
#else
  if ( use_vectorization )
  {
    BlueSolveBackgroundTrace_Vec ( my_flow.numfb, vb_out, sbgPtr, time_c.npts(), &time_c.deltaFrame[0],
                                   my_scratch.cur_buffer_block.tauB, my_scratch.cur_buffer_block.etbR );
  }
  else
  {
    for ( int fnum=0; fnum<my_flow.numfb; fnum++ )
      BlueSolveBackgroundTrace ( vb_out[fnum],sbgPtr[fnum],time_c.npts(),&time_c.deltaFrame[0],
                                 my_scratch.cur_buffer_block.tauB[fnum],my_scratch.cur_buffer_block.etbR[fnum] );
  }
#endif

  MultiplyVectorByScalar ( vb,p->gain,my_scratch.bead_flow_t );

  ApplyDarkMatter ( vb,&my_regions.rp, my_regions.missing_mass.dark_matter_compensator,my_flow.flow_ndx_map,time_c.npts() );

  // zeromer computed, now remove from observed
  DiminishVector ( block_signal_corrected,vb,my_scratch.bead_flow_t ); // remove calculated background to produce corrected signal

}

void AccumulateSingleNeighborExcessHydrogenOneParameter ( float *my_xtflux, float *neighbor_signal,
    bead_params *p, reg_params *reg_p,
    BeadScratchSpace &my_scratch, TimeCompression &time_c,
    RegionTracker &my_regions, flow_buffer_info my_flow,
    bool use_vectorization,
    float multiplier, bool rescale_flag )
{

  FillBufferParamsBlockFlows ( &my_scratch.cur_buffer_block,p,reg_p,my_flow.flow_ndx_map,my_flow.buff_flow );

  float block_model_trace[my_scratch.bead_flow_t], block_incorporation_rise[my_scratch.bead_flow_t];

  // "In parallel, across flows"

  float* model_trace[NUMFB];
  float* incorporation_rise[NUMFB];
  float* neighbor_local[NUMFB];

  // should this be using cur_buffer_block as usual?
  float vec_tau_well[NUMFB];
  float vec_tau_empty[NUMFB];

  for ( int fnum=0; fnum<NUMFB; fnum++ )
  {
    vec_tau_well[fnum] = my_scratch.cur_buffer_block.tauB[fnum];
    vec_tau_empty[fnum] = my_scratch.cur_buffer_block.etbR[fnum]*vec_tau_well[fnum];
  }

  // setup parallel pointers into the structure
  for ( int fnum=0; fnum<NUMFB; fnum++ )
  {

    model_trace[fnum] = &block_model_trace[fnum*time_c.npts() ];
    incorporation_rise[fnum] = &block_incorporation_rise[fnum*time_c.npts() ];   // set up each flow information
    neighbor_local[fnum] = &neighbor_signal[fnum*time_c.npts() ];
  }

  IncorporationRiseFromNeighborSignal ( incorporation_rise,neighbor_local, time_c,my_regions, my_scratch );

  CumulativeLostHydrogens ( model_trace,incorporation_rise, time_c, my_regions, use_vectorization, vec_tau_well );

  IonsFromBulk ( model_trace,incorporation_rise, time_c, my_regions, my_flow, use_vectorization, vec_tau_empty );
  if (rescale_flag)
    RescaleTraceByBuffering(model_trace, vec_tau_well, vec_tau_empty,time_c.npts());

  // universal
  MultiplyVectorByScalar ( block_model_trace,multiplier,my_scratch.bead_flow_t ); // scale down the quantity of ions

  // add to the bulk cross-talk we're creating
  AccumulateVector ( my_xtflux, block_model_trace,my_scratch.bead_flow_t );
}

void MultiFlowComputeCumulativeIncorporationSignal ( struct bead_params *p,struct reg_params *reg_p, float *ivalPtr,
    NucStep &cache_step, incorporation_params_block_flows &cur_bead_block,
    TimeCompression &time_c, flow_buffer_info &my_flow,  PoissonCDFApproxMemo *math_poiss )
{
  // only side effect should be new values in ivalPtr

  // this is region wide
  //This will short-circuit if has been computed
  cache_step.CalculateNucRiseCoarseStep ( reg_p,time_c,my_flow );
  // pretend I'm making a parallel process
  FillIncorporationParamsBlockFlows ( &cur_bead_block, p,reg_p,my_flow.flow_ndx_map,my_flow.buff_flow );
  // "In parallel, across flows"
  float* nuc_rise_ptr[NUMFB];
  float* incorporation_rise[NUMFB];
  int my_start[NUMFB];
  for ( int fnum=0; fnum<NUMFB; fnum++ )
  {
    nuc_rise_ptr[fnum]      = cache_step.NucCoarseStep ( fnum );
    incorporation_rise[fnum]= &ivalPtr[fnum*time_c.npts() ];
    my_start[fnum] = cache_step.i_start_coarse_step[fnum];
  }
  // this is >almost< a parallel operation by flows now
  bool use_my_parallel=true;
  if ( use_my_parallel )
  {
    //@TODO handle cases of fewer than 4 flows remaining
    for ( int fnum=0; fnum<NUMFB; fnum+=4 )
    {

      ParallelSimpleComputeCumulativeIncorporationHydrogens ( &incorporation_rise[fnum], time_c.npts(), &time_c.deltaFrameSeconds[0],&nuc_rise_ptr[fnum],
          ISIG_SUB_STEPS_MULTI_FLOW,  &my_start[fnum], &p->Ampl[fnum],
          &cur_bead_block.SP[fnum],&cur_bead_block.kr[fnum], &cur_bead_block.kmax[fnum], &cur_bead_block.d[fnum], &cur_bead_block.molecules_to_micromolar_conversion[fnum], math_poiss );
      for ( int q=0; q<4; q++ )
        MultiplyVectorByScalar ( incorporation_rise[fnum+q], cur_bead_block.sens[fnum+q],time_c.npts() );  // transform hydrogens to signal
    }
  }
  else
  {
    for ( int fnum=0;fnum<NUMFB;fnum++ )
    {
      ComputeCumulativeIncorporationHydrogens ( incorporation_rise[fnum], time_c.npts(),
          &time_c.deltaFrameSeconds[0], nuc_rise_ptr[fnum], ISIG_SUB_STEPS_MULTI_FLOW, my_start[fnum],
          cur_bead_block.C[fnum], p->Ampl[fnum], cur_bead_block.SP[fnum], cur_bead_block.kr[fnum], cur_bead_block.kmax[fnum], cur_bead_block.d[fnum],cur_bead_block.molecules_to_micromolar_conversion[fnum], math_poiss );

      MultiplyVectorByScalar ( incorporation_rise[fnum], cur_bead_block.sens[fnum],time_c.npts() );  // transform hydrogens to signal
    }
  }
}

void extern_links::DumpTimeAndEmphasisByRegionH5 ( int reg, TimeCompression &time_c, EmphasisClass &emphasis_data )
{
  if ( mPtrs->mEmphasisParam!=NULL )
  {

    ION_ASSERT ( emphasis_data.numEv <= MAX_HPLEN+1, "emphasis_data.numEv > MAX_HPLEN+1" );
    //ION_ASSERT(time_c.npts <= MAX_COMPRESSED_FRAMES, "time_c.npts > MAX_COMPRESSED_FRAMES");
    int npts = std::min ( time_c.npts(), MAX_COMPRESSED_FRAMES );

    if ( mPtrs->mEmphasisParam!=NULL )
    {
      // use copyCube_element to copy DataCube element to mEmphasisParam in BkgModel::DumpTimeAndEmphasisByRegionH5
      for ( int hp=0; hp<emphasis_data.numEv; hp++ )
      {
        for ( int t=0; t< npts; t++ )
          mPtrs->copyCube_element ( mPtrs->mEmphasisParam,reg,hp,t,emphasis_data.EmphasisVectorByHomopolymer[hp][t] );
        for ( int t=npts; t< MAX_COMPRESSED_FRAMES; t++ )
          mPtrs->copyCube_element ( mPtrs->mEmphasisParam,reg,hp,t,0 ); // pad 0's, memset faster here?
      }

      for ( int hp=emphasis_data.numEv; hp<MAX_HPLEN+1; hp++ )
      {
        for ( int t=0; t< MAX_COMPRESSED_FRAMES; t++ )
        {
          mPtrs->copyCube_element ( mPtrs->mEmphasisParam,reg,hp,t,0 ); // pad 0's, memset faster here?
        }
      }
    }
  }
}

// note: input is incorporation_rise, returns lost_hydrogens in the same buffer, recycling the memory
void CumulativeLostHydrogens ( float **incorporation_rise_to_lost_hydrogens, float **scratch_trace, 
                               TimeCompression &time_c, RegionTracker &my_regions,
                               bool use_vectorization,
                               float *vec_tau_top )
{
  // Put the model trace from buffering the incorporation_rise into scratch_trace
#ifdef __INTEL_COMPILER
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
    {
      RedSolveHydrogenFlowInWell ( scratch_trace[fnum],incorporation_rise_to_lost_hydrogens[fnum],time_c.npts(),my_regions.cache_step.i_start_coarse_step[fnum],&time_c.deltaFrame[0],vec_tau_top[fnum] ); // we lose hydrogen ions fast!

    }
  }
#else
  if ( use_vectorization )
  {
    RedSolveHydrogenFlowInWell_Vec ( NUMFB,scratch_trace,incorporation_rise_to_lost_hydrogens,time_c.npts(),&time_c.deltaFrame[0],vec_tau_top ); // we lose hydrogen ions fast!

  }
  else
  {
    for ( int fnum=0; fnum<NUMFB; fnum++ )
    {
      RedSolveHydrogenFlowInWell ( scratch_trace[fnum],incorporation_rise_to_lost_hydrogens[fnum],time_c.npts(),my_regions.cache_step.i_start_coarse_step[fnum],&time_c.deltaFrame[0],vec_tau_top[fnum] ); // we lose hydrogen ions fast!

    }
  }
#endif

  // return lost_hydrogens in the incorporation_rise variables by subtracting the trace from the cumulative
  for ( int fnum=0; fnum<NUMFB; fnum++ )
    DiminishVector ( incorporation_rise_to_lost_hydrogens[fnum],scratch_trace[fnum],time_c.npts() );  // cumulative lost hydrogen ions instead of retained hydrogen ions

}

static void IncorporationRiseFromNeighborParameters ( 
    float **incorporation_rise, const float * const *nuc_rise_ptr,
    BeadParams *p,
    TimeCompression &time_c, RegionTracker &my_regions,
    BeadScratchSpace &my_scratch, 
    incorporation_params_block_flows &my_cur_bead_block,
    PoissonCDFApproxMemo *math_poiss, int hydrogenModelType,
    int flow_block_size)
{
  // fill in each flow incorporation
  for ( int fnum=0; fnum<flow_block_size; fnum++ )
  {
    // compute newly generated ions for the amplitude of each flow
    MathModel::ComputeCumulativeIncorporationHydrogens ( incorporation_rise[fnum],
        time_c.npts(), &time_c.deltaFrameSeconds[0], nuc_rise_ptr[fnum], ISIG_SUB_STEPS_MULTI_FLOW, my_regions.cache_step.i_start_coarse_step[fnum],
        my_cur_bead_block.C[fnum], p->Ampl[fnum], my_cur_bead_block.SP[fnum], my_cur_bead_block.kr[fnum], my_cur_bead_block.kmax[fnum], my_cur_bead_block.d[fnum],my_cur_bead_block.molecules_to_micromolar_conversion[fnum], math_poiss, hydrogenModelType );
    MultiplyVectorByScalar ( incorporation_rise[fnum], my_cur_bead_block.sens[fnum],time_c.npts() );  // transform hydrogens to signal       // variables used for solving background signal shape
  }
}

void IncorporationRiseFromNeighborSignal ( float **incorporation_rise, float **neighbor_local,
    TimeCompression &time_c, RegionTracker &my_regions,
    BeadScratchSpace &my_scratch
                                         )
{
  // fill in each flow incorporation
  for ( int fnum=0; fnum<NUMFB; fnum++ )
  {
    // grab my excess hydrogen ions from the observed signal
    IntegrateRedFromObservedTotalTrace ( incorporation_rise[fnum], neighbor_local[fnum], &my_scratch.shifted_bkg[fnum*time_c.npts() ], time_c.npts(),&time_c.deltaFrame[0], my_scratch.cur_buffer_block.tauB[fnum], my_scratch.cur_buffer_block.etbR[fnum] );
    // variables used for solving background signal shape
  }

  // construct the incorporation in the neighbor well just for research purposes
  for ( int fnum=0; fnum<NUMFB; fnum++ )
  {
    // bad! put trace back in neighbor signal
    RedSolveHydrogenFlowInWell ( neighbor_local[fnum],incorporation_rise[fnum],time_c.npts(),my_regions.cache_step.i_start_coarse_step[fnum],&time_c.deltaFrame[0],my_scratch.cur_buffer_block.tauB[fnum] ); // we lose hydrogen ions fast!
  }
}