// TODO: declare audio ended so we can ignore any incomplete bass frequencies

sclass ContinuousOscillators_v2 extends Meta {
  macro dprint { if (debug) printVars } bool debug;
    
  settable bool bassFirst;
  settable bool printSubtractions;
  int channels = 1;
  settable double sampleRate = 48000;
  double currentSample;
  IAudioSample audio;              // integrated audio sample
  SubtractedAudio subtractedAudio; // audio with deleted frequencies
  settable bool keepRendered;
  transient double[] rendered;
  
  bool renderUsingSine; // or square wave
  bool realOnly, imagOnly; // for testing
  double maxRelFreq = 0.5; // maximum oscillator frequency relative to sample rate
  
  event escalating(Oscillator o);

  long haarFeaturesChecked;

  persistable class Oscillator extends HasKey<Frequency> {
    double currentPeriodStart;
    // all recorded intensities
    new TreeMap<DoubleRange, Channels<Complex>> intensities;
    double interval; // length of period in samples
    
    *(Frequency f) {
      super(f);
      interval = frequency().interval()*sampleRate;
    }
    
    Frequency frequency() { ret key(); }
    double interval() { ret interval; }
    
    scaffolded void record_impl(DoubleRange r, Channels<Complex> c) {
      intensities.put(r, c);
      subtractMeasured(r, c);
      currentPeriodStart = phase360();
      if (!bassFirst) escalate();
    }
    
    void escalate {
      escalating(this);
      rcallF(oscillators.higher(this), o -> { o?.process(currentPeriodStart); });
    }
    
    void subtractMeasured(DoubleRange r, Channels<Complex> c) {
      renderer from method;
      if (printSubtractions) {
        print("subtract " + r + " = " + c);
        enableScaffolding(renderer);
      }
      renderer.subtractPeriodFrom(subtractedAudio, r, frequency(), c.getSingleton());
      renderer().addPeriodTo(rendered, r, frequency(), c.getSingleton());
    }
    
    scaffolded void process(double upTo) {
      //printVars(+phase360(), +upTo);
      while (phase360() <= upTo)
        measure();
      subtractedAudioSanityCheck();
      escalate();
    }
    
    scaffolded void measure() {
      DoubleRange p1 = doubleRangeWithLength(currentPeriodStart, interval);
      AudioHaarFeature haar = new(relevantAudio(), p1);
      haarFeaturesChecked += 2;
      record(this, p1, haar.getComplex());
    }
    
    double phase360() { ret currentPeriodStart + interval(); }
  } // end of Oscillator
  
  IAudioSample relevantAudio() { ret subtractedAudio; }
  
  TreeHasKeyMap<Frequency, Oscillator> oscillators;

  *() {}
  *(IAudioSample audio, double *currentSample) { setAudio(audio); }
  *(IAudioSample audio) { setAudio(audio); } // start at time 0

  void setFrequencies(Iterable<Frequency> frequencies) {
    oscillators = new TreeHasKeyMap<Frequency, Oscillator>(map(f -> new Oscillator(f), filter(frequencies,
      f -> between(f!, 1, sampleRate*maxRelFreq))));
  }
  
  void setAudio(IAudioSample audio) {
    this.audio = audio;
    subtractedAudio = new SubtractedAudio(audio);
    if (keepRendered) rendered = new double[iceil(audio.length())];
    assertEquals("bounds", audio.bounds(), subtractedAudio.bounds());
  }
  
  int nFrequencies() { ret oscillators.size(); }
  
  scaffolded void stepTo(DoubleRange timeRange) {
    var o = first(oscillators);
    o?.process(timeRange.end);
  }
  
  void record(Oscillator o, DoubleRange r, Channels<Complex> c) {
    if (o == null) ret;
    o.record_impl(r, c);
  }
  
  Frequency lowestFrequency() { ret oscillators.firstKey(); }
  Frequency highestFrequency() { ret oscillators.lastKey(); }
  
  int minWindowSize() {
    ret iceil(sampleRate*lowestFrequency().interval());
  }
  
  S stats() { ret renderVars(+audio, +haarFeaturesChecked); }
  
  // how much horizontal area have we covered
  L<Double> coverageByFreq() {
    ret mapReversed(oscillators, o -> totalLengthOfDoubleRanges(keys(o.intensities)));
  }
  
  IIntegralImage imageColumn(DoubleRange timeRange) {
    int n = nFrequencies();
    double[] col = rawImageColumn(timeRange);
    double[] col2 = normalizeDoubles(col, 255);
    ret bwIntegralImageFromFunction(1, n, (x, y) -> ifloor(col2[y]));
  }
  
  double[] rawImageColumn(DoubleRange timeRange) {
    int n = nFrequencies();
    var l = reversedList(oscillators);
    double[] col = new[n];
    for y to n: {
      TreeMap<DoubleRange, Channels<Complex>> l2 = l.get(y).intensities;
      
      // find right period
      DoubleRange period = l2.floorKey(timeRange);
      if (period == null) continue;
      
      double amplitude = l2.get(period).first().abs();
      col[y] = amplitude;
    }
    
    ret col;
  }
  
  MakesBufferedImage image(DoubleRange r, double pixelsPerSecond) {
    new L<IIntegralImage> images;
    double len = r.length()/sampleRate;
    int pixels = iround(len*pixelsPerSecond);
    for i to pixels: {
      DoubleRange range = doubleRange(
        i/pixelsPerSecond*sampleRate+r.start(),
        (i+1)/pixelsPerSecond*sampleRate+r.start());
      //print(+range);
      images.add(imageColumn(range));
    }
    
    ret mergeBWImagesHorizontally(map toBWImage(images), spacing := 0);
  }
  
  double[] toAudio(IntRange sampleRange, Iterable<Oscillator> oscillatorsToRender default oscillators) {
    double[] samples = new[l(sampleRange)];
    
    class PerFreq {
      Oscillator o;
      DoubleRange period;
      Channels<Complex> intensity;

      *(Oscillator *o) {
        setPeriod(o.intensities.floorKey(toDoubleRange(sampleRange)));
        //printVars("PerFreq", f := o.frequency(), +sampleRange, +period, firstPeriod := firstKey(o.intensities));
      }
      
      void setPeriod(DoubleRange period) {
        this.period = period;
        intensity = period == null ?: o.intensities.get(period);
      }
      
      void nextPeriod(double t) {
        while (period != null && t >= period.end)
          setPeriod(o.intensities.higherKey(period));
      }
    }
    
    L<PerFreq> perFreqs = map(oscillatorsToRender, o -> new PerFreq(o));
    
    for i over samples: {
      double t = sampleRange.start+i; // t is time in samples
      double sum = 0;
      for (pf : perFreqs) {
        pf.nextPeriod(t);
        bool shouldPrint = toAudio_shouldPrint(i, l(samples));
        
        if (pf.intensity != null) {
          sum += renderer().renderPeriodSample(pf.period, pf.o.frequency(), pf.intensity.first(), t, shouldPrint);
        } else
          if (shouldPrint)
            printVars(f := pf.o.frequency(), +i, interval := pf.o.interval, +t, c := 0);
        
      }
      
      samples[i] = sum;
    }
  
    ret samples;
  }
  
  RenderFrequencySample renderer() {
    ret renderUsingSine
      ? new RenderFrequencySample_Sine
      : new RenderFrequencySample_SquareWave;
  }

  swappable bool toAudio_shouldPrint(int i, int n) { false; }
  
  void subtractedAudioSanityCheck() {
    if (rendered == null) ret;
    int n = l(rendered);
    for i to n: {
      reimmutable i;
      assertDoublesVeryClose_warnOnly(-> "Sample " + i,
        subtractedAudio.sampleSum(0, i, i+1),
        audio.sampleSum(0, i, i+1)-rendered[i],
        1e-5);
    }
  }
}

Began life as a copy of #1033245