Rhs2116 - Search through possible step sizes to find global step size

[bparks13]

Instead of choosing the smallest step size for a new amplitude, search through all existing amplitudes and all possible step sizes to see if one fits all amplitudes simultaneously.
Modified the existing error message to be more clear why the current amplitudes are incompatible with the new step size.

• Fixes Instead of using the smallest step size, search for a step size that works for all stimuli #345

[bparks13]

Commenting here the algorithm used to define the step-size.

1. User types in an amplitude value [0.01 - 2550 µA]
   1. These bounds are checked to ensure it is a valid amplitude
2. A list of possible step sizes are found that can be used to generate this amplitude
   1. If the list contains the current step size that other channels use, use the current step size
      1. Note: This currently does not evaluate error for the requested amplitude. Should the newly requested amplitude take priority over existing amplitudes? Meaning, should we change the step size if the new amplitude is slightly off, or should we prioritize the existing amplitudes and round the new amplitude to match?
   2. Otherwise, create a new list of step sizes that are valid for all existing amplitudes
      1. Using this new list, iterate over all existing and new amplitudes to calculate the error, and find a step size that minimizes the maximum error of each amplitude
   3. If no step size exists that fits all amplitudes, choose the largest step size that is valid for the new amplitude, but give a warning that the step size is out of range
      1. Note: This does not change the step size for the existing amplitudes yet, the step size for the existing amplitudes is only updated if the user tries to apply the stimulus settings to a new channel

[jonnew]

Commenting here the algorithm used to define the step-size.

1. User types in an amplitude value [0.01 - 2550 µA]

   1. These bounds are checked to ensure it is a valid amplitude

good. Can be done by using numeric up down instead of text box.

Remove this step [

2. A list of possible step sizes are found that can be used to generate this amplitude

   1. If the list contains the current step size that other channels use, use the current step size

      1. Note: This currently does not evaluate error for the requested amplitude. Should the newly requested amplitude take priority over existing amplitudes? Meaning, should we change the step size if the new amplitude is slightly off, or should we prioritize the existing amplitudes and round the new amplitude to match?

]

Start here:

2. Otherwise, create a new list of step sizes that are valid for all existing amplitudes

   1. Using this new list, iterate over all existing and new amplitudes to calculate the error, and find a step size that minimizes the maximum error of each amplitude

3. If no step size exists that fits all amplitudes, choose the largest step size that is valid for the new amplitude, but give a warning that the step size is out of range

   1. Note: This does not change the step size for the existing amplitudes yet, the step size for the existing amplitudes is only updated if the user tries to apply the stimulus settings to a new channel

The last point is why this process should be uniform. There are no special cases: everything is filtered through this Minimax scheme.

[jonnew]

Unfortunately, this is still incorrect:

Here I'm trying to change one amplitude to 100 uA while the rest of the channels are at 1000 uA. 0.5 uA is selected and I get the following warning:

But both 5 uA and 10 uA steps would perfectly fit 100 and 1000 uA amplitudes.

[jonnew]

I'm attaching a starting point for correcting this, potentially. I would completely forget the GUI and think about how to modify the definitions of Rhs2116Stimulus, Rhs2116StimulusSequence, and Rhs2116StimulusSequencePair in order to verify a requested stimulus and optimize the step size.

• I would make Rhs2116StimulusSequence.StepSize have a private setter to ensure its only going to be manipulated via internal logic
• I would make Rhs2116StimulusSequence.Stimuli have private setter as well.
• I would then provide this class with a public function that can be called to try to add stimuli to the Stimuli while doing verification, error optimization, etc.

** This is a starting point, the logic is not tested **

using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Xml.Serialization;
using System.Runtime.CompilerServices;

[assembly: InternalsVisibleTo("OpenEphys.Onix1.Design")]

namespace OpenEphys.Onix1
{
    /// <summary>
    /// A stimulus sequence for a Rhs2116 device.
    /// </summary>
    public class Rhs2116StimulusSequence
    {
        internal const int NumberOfChannels = 16;
        const int NumberOfStepSizes = 10;

        /// <summary>
        /// Initializes a new instance of the <see cref="Rhs2116StimulusSequence"/> class with 16 default
        /// stimuli.
        /// </summary>
        public Rhs2116StimulusSequence()
        {
            Stimuli = new Rhs2116Stimulus[NumberOfChannels];

            for (var i = 0; i < Stimuli.Length; i++)
            {
                Stimuli[i] = new Rhs2116Stimulus();
            }
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Rhs2116StimulusSequence"/> by performing a
        /// shallow copy of a reference sequence.
        /// </summary>
        /// <param name="sequence">Existing <see cref="Rhs2116StimulusSequence"/> to copy.</param>
        public Rhs2116StimulusSequence(Rhs2116StimulusSequence sequence)
        {
            Stimuli = Array.ConvertAll(sequence.Stimuli, stimulus => stimulus.Clone());
            CurrentStepSize = sequence.CurrentStepSize;
        }


        public bool DefineStimulus(int channel, uint numberOfStimuli, bool anodicFirst, double delayMillis, double dwellMillis, 
            double anodicMillis, double anodicMicroAmps, double cathodicMillis, double cathodicMicroAmps, double interStimMillis,
            bool force)
        {
            // TODO: sanity checks other than amplitude checking
            var error = new double[NumberOfChannels * 2, NumberOfStepSizes];
            var maxError = new List<double>(NumberOfStepSizes);
            var stepSizes = Enum.GetValues(typeof(Rhs2116StepSize)).Cast<Rhs2116StepSize>();
            var stepSizesuA = stepSizes.Select(s => GetStepSizeuA(s)).ToArray();
            var currentStepSizeuA = GetStepSizeuA(CurrentStepSize);

            for (int c = 0; c < NumberOfChannels; c += 2)
            {
                var anodicAmp = c == channel ? anodicMicroAmps : Stimuli[c].AnodicAmplitudeSteps * currentStepSizeuA;
                var cathodicAmp = c == channel ? cathodicMicroAmps : Stimuli[c].CathodicAmplitudeSteps * currentStepSizeuA;

                for (int s = 0; s < NumberOfStepSizes; s++)
                {
                    error[c, s] = anodicAmp < stepSizesuA[s] || anodicAmp > s * 255 ?
                        double.PositiveInfinity :
                        Math.Abs((anodicAmp - (stepSizesuA[s] * Math.Round(anodicAmp / stepSizesuA[s]))) / anodicAmp);

                    error[c + 1, s] = cathodicAmp < stepSizesuA[s] || cathodicAmp > s * 255 ?
                        double.PositiveInfinity :
                        Math.Abs((cathodicAmp - (stepSizesuA[s] * Math.Round(cathodicAmp / stepSizesuA[s]))) / cathodicAmp);

                    maxError[s] = Math.Max(maxError[s], Math.Max(error[c, s], error[c + 1, s]));

                }
            }

            var optimalStep = maxError.ToList().IndexOf(maxError.Min());
            var optimalError = maxError[optimalStep];

            if (optimalError == double.PositiveInfinity && !force)
                return false;

            CurrentStepSize = stepSizes.ElementAt(optimalStep);
            // TODO: regenerate Stimuli by assigning Stimuli[channel] and updated the others' parameters with the new step size

            return true;
        }

        /// <summary>
        /// Gets or sets the array of <see cref="Rhs2116Stimulus"/> sequences.
        /// </summary>
        public Rhs2116Stimulus[] Stimuli { get; private set; }

        /// <summary>
        /// Gets or sets the <see cref="Rhs2116StepSize"/>. 
        /// </summary>
        public Rhs2116StepSize CurrentStepSize { get; private set; } = Rhs2116StepSize.Step5000nA;

        /// <summary>
        /// Maximum length of the sequence across all channels.
        /// </summary>
        [XmlIgnore]
        public uint SequenceLengthSamples
        {
            get
            {
                uint max = 0;

                foreach (var stim in Stimuli)
                {
                    var len = stim.DelaySamples + stim.NumberOfStimuli * (stim.AnodicWidthSamples + stim.CathodicWidthSamples + stim.DwellSamples + stim.InterStimulusIntervalSamples);
                    max = len > max ? len : max;

                }

                return max;
            }
        }

        /// <summary>
        /// Maximum peak to peak amplitude of the sequence across all channels.
        /// </summary>
        [XmlIgnore]
        public int MaximumPeakToPeakAmplitudeSteps
        {
            get
            {
                int max = 0;

                foreach (var stim in Stimuli)
                {
                    var p2p = stim.CathodicAmplitudeSteps + stim.AnodicAmplitudeSteps;
                    max = p2p > max ? p2p : max;

                }

                return max;
            }
        }

        /// <summary>
        /// Is the stimulus sequence well defined.
        /// </summary>
        [XmlIgnore]
        public bool Valid => Stimuli.ToList().All(s => s.Valid);

        /// <summary>
        /// Does the sequence fit in hardware.
        /// </summary>
        [XmlIgnore]
        public bool FitsInHardware => StimulusSlotsRequired <= Rhs2116.StimMemorySlotsAvailable;

        /// <summary>
        /// The maximum number of memory slots available.
        /// </summary>
        [XmlIgnore]
        public int MaxMemorySlotsAvailable => Rhs2116.StimMemorySlotsAvailable;

        /// <summary>
        /// Number of hardware memory slots required by the sequence.
        /// </summary>
        [XmlIgnore]
        public int StimulusSlotsRequired => DeltaTable.Count;

        /// <summary>
        /// Gets the current step size in µA.
        /// </summary>
        [XmlIgnore]
        public double CurrentStepSizeuA
        {
            get => GetStepSizeuA(CurrentStepSize);
        }

        /// <summary>
        /// Gets the step size in microamps of a given <see cref="Rhs2116StepSize"/>.
        /// </summary>
        /// <param name="stepSize">The <see cref="Rhs2116StepSize"/> value to convert to microamps.</param>
        /// <returns>Returns a double corresponding to the given step size in microamps.</returns>
        /// <exception cref="ArgumentException">Invalid stimulus step size selection.</exception>
        public static double GetStepSizeuA(Rhs2116StepSize stepSize)
        {
            return stepSize switch
            {
                Rhs2116StepSize.Step10nA => 0.01,
                Rhs2116StepSize.Step20nA => 0.02,
                Rhs2116StepSize.Step50nA => 0.05,
                Rhs2116StepSize.Step100nA => 0.1,
                Rhs2116StepSize.Step200nA => 0.2,
                Rhs2116StepSize.Step500nA => 0.5,
                Rhs2116StepSize.Step1000nA => 1.0,
                Rhs2116StepSize.Step2000nA => 2.0,
                Rhs2116StepSize.Step5000nA => 5.0,
                Rhs2116StepSize.Step10000nA => 10.0,
                _ => throw new ArgumentException("Invalid stimulus step size selection."),
            };
        }

        /// <summary>
        /// Gets the maximum possible amplitude for a single pulse, in µA.
        /// </summary>
        [XmlIgnore]
        public double MaxPossibleAmplitudePerPhaseMicroAmps => CurrentStepSizeuA * 255;

        internal IEnumerable<byte> AnodicAmplitudes => Stimuli.ToList().Select(x => x.AnodicAmplitudeSteps);

        internal IEnumerable<byte> CathodicAmplitudes => Stimuli.ToList().Select(x => x.CathodicAmplitudeSteps);

        /// <summary>
        /// Generate the delta-table representation of this stimulus sequence that can be uploaded to the
        /// RHS2116 device. The resultant dictionary has a time, in samples, as the key and a combined
        /// [polarity, enable] bit field as the value.
        /// </summary>
        [XmlIgnore]
        internal Dictionary<uint, uint> DeltaTable
        {
            get => GetDeltaTable();
        }

        private Dictionary<uint, uint> GetDeltaTable()
        {
            var table = new Dictionary<uint, BitArray>();

            for (int i = 0; i < NumberOfChannels; i++)
            {
                var s = Stimuli[i];

                var e0 = s.AnodicFirst ? s.AnodicAmplitudeSteps > 0 : s.CathodicAmplitudeSteps > 0;
                var e1 = s.AnodicFirst ? s.CathodicAmplitudeSteps > 0 : s.AnodicAmplitudeSteps > 0;
                var d0 = s.AnodicFirst ? s.AnodicWidthSamples : s.CathodicWidthSamples;
                var d1 = d0 + s.DwellSamples;
                var d2 = d1 + (s.AnodicFirst ? s.CathodicWidthSamples : s.AnodicWidthSamples);

                var t0 = s.DelaySamples;

                for (int j = 0; j < s.NumberOfStimuli; j++)
                {
                    AddOrInsert(ref table, i, t0, s.AnodicFirst, e0);
                    AddOrInsert(ref table, i, t0 + d0, s.AnodicFirst, false);
                    AddOrInsert(ref table, i, t0 + d1, !s.AnodicFirst, e1);
                    AddOrInsert(ref table, i, t0 + d2, !s.AnodicFirst, false);

                    t0 += d2 + s.InterStimulusIntervalSamples;
                }
            }

            return table.ToDictionary(d => d.Key, d =>
            {
                int[] i = new int[1];
                d.Value.CopyTo(i, 0);
                return (uint)i[0];
            });
        }

        private static void AddOrInsert(ref Dictionary<uint, BitArray> table, int channel, uint key, bool polarity, bool enable)
        {
            if (table.ContainsKey(key))
            {
                table[key][channel] = enable;
                table[key][channel + 16] = polarity;
            }
            else
            {
                table.Add(key, new BitArray(32, false));
                table[key][channel] = enable;
                table[key][channel + 16] = polarity;
            }
        }
    }
}

[jonnew]

I terms of behavior, this looks very good. I played around with several extreme combinations of amplitudes

• At margins of min and max for an step size i had in mind
• Looking at the rounding behavior versus actual applied amplitude in grey box
• Looking at display while i was doing this.

All of that looked good. However, the underlying code does need work.

[bparks13]

@jonnew I have addressed the review comments. Let me know if this PR is approved, and if it is I will go ahead and rebase to clean up the history before merging

[jonnew]

Looks good. Can we add the following functionality before the merge please: