Created
August 3, 2024 15:23
-
-
Save acaly/127af291d0f5747fea1cc29aa4457cac to your computer and use it in GitHub Desktop.
Orthogonal fractional factorial experiment design
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using System.Diagnostics; | |
using System.Runtime.CompilerServices; | |
using System.Runtime.InteropServices; | |
namespace ExperimentDesign | |
{ | |
internal class Program | |
{ | |
static void Main(string[] args) | |
{ | |
Console.WriteLine("DOE orthogonal plan generator."); | |
Console.WriteLine("Input number of levels for each factors (one per line):"); | |
List<uint> levels = []; | |
while (true) | |
{ | |
var line = Console.ReadLine(); | |
if (!uint.TryParse(line, out var level) || level == 0) | |
{ | |
if (!string.IsNullOrWhiteSpace(line)) | |
{ | |
Console.WriteLine("Invalid number, assuming end of input"); | |
} | |
break; | |
} | |
levels.Add(level); | |
} | |
if (levels.Count == 0) | |
{ | |
Console.WriteLine("No input provided."); | |
return; | |
} | |
var plan = new Plan(levels); | |
plan.Generate(); | |
var result = new uint[levels.Count]; | |
for (uint i = 0; i < plan.SampleCount; ++i) | |
{ | |
plan.GetSamples(i, result); | |
foreach (var v in result) | |
{ | |
Console.Write(v); | |
Console.Write('\t'); | |
} | |
Console.WriteLine(); | |
} | |
} | |
internal abstract class OrthogonalArray | |
{ | |
protected readonly uint _base, _power, _size; | |
public abstract void SetColumns(uint n); | |
public abstract uint GetCell(uint row, uint column); | |
public static OrthogonalArray Create(uint nbase, uint power, uint size) | |
{ | |
Debug.Assert(Math.Pow(nbase, power) == size); | |
ArgumentOutOfRangeException.ThrowIfGreaterThan(nbase, 255u); | |
if (power == 1) | |
{ | |
return new SingleColumnArray(nbase, power, size); | |
} | |
else | |
{ | |
return new GeneralArray(nbase, power, size); | |
} | |
} | |
private OrthogonalArray(uint nbase, uint power, uint size) | |
{ | |
_base = nbase; | |
_power = power; | |
_size = size; | |
} | |
private sealed class SingleColumnArray : OrthogonalArray | |
{ | |
public SingleColumnArray(uint nbase, uint power, uint size) | |
: base(nbase, power, size) | |
{ | |
Debug.Assert(power == 1); | |
} | |
public override void SetColumns(uint n) | |
{ | |
ArgumentOutOfRangeException.ThrowIfNotEqual(n, 1u); | |
} | |
public override uint GetCell(uint row, uint column) | |
{ | |
ArgumentOutOfRangeException.ThrowIfNotEqual(column, 0u); | |
ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(row, _base); | |
return row; | |
} | |
} | |
/* | |
* Example for 3^3=27 trials experiments, which can holds 13 factors, as 13*(3-1)=27-1. | |
* Three elementary factors are called A, B, and C. | |
* Other factors are combinations of them. | |
* | |
* C | |
* B, B+C, B+2C | |
* A, A+B, A+2B, A+C, A+B+C, A+2B+C, A+2C, A+B+2C, A+2B+2C | |
* | |
* One can confirm any pair among these factors are orthogonal. | |
* The columns are reordered vertically: C, B, A, B+C, A+B, B+2C, A+2B, A+C, ... | |
* This can avoid adjacent factors to be correlated. | |
* | |
* In the implementation below, A+B+2C is represented as [2,1,1]. | |
* Representations for all columns are concatenated and saved in the list for later use. | |
*/ | |
private sealed class GeneralArray : OrthogonalArray | |
{ | |
private readonly List<byte> _columns = []; | |
private readonly uint _maxColumns; | |
private uint _columnCount = 0; | |
public GeneralArray(uint nbase, uint power, uint size) | |
: base(nbase, power, size) | |
{ | |
_maxColumns = (size - 1) / (nbase - 1); | |
} | |
public override void SetColumns(uint n) | |
{ | |
ArgumentOutOfRangeException.ThrowIfGreaterThan(n, _maxColumns); | |
_ = checked((int)(n * _power)); | |
_columns.Clear(); | |
//Implement the reordering of columns. | |
var power = _power; | |
var nbase = _base; | |
var coeff = new byte[power + 1]; | |
var startingLine = 0; | |
while ((uint)_columns.Count < n * power) | |
{ | |
for (int line = startingLine; line < power; ++line) | |
{ | |
for (int i = startingLine; i < line; ++i) | |
{ | |
_columns.Add(0); | |
} | |
for (int i = 0; i < startingLine; ++i) | |
{ | |
_columns.Add(coeff[i]); | |
} | |
_columns.Add(1); | |
for (int i = line + 1; i < power; ++i) | |
{ | |
_columns.Add(0); | |
} | |
} | |
int inc = 0; | |
while (true) | |
{ | |
if (++coeff[inc] < nbase) | |
{ | |
break; | |
} | |
coeff[inc] = 0; | |
inc += 1; | |
} | |
if (inc >= startingLine) | |
{ | |
startingLine = inc + 1; | |
} | |
} | |
_columnCount = n; | |
} | |
public override uint GetCell(uint row, uint column) | |
{ | |
ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(row, _size); | |
ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(column, _columnCount); | |
var coeff = CollectionsMarshal.AsSpan(_columns)[(int)(column * _power)..]; | |
uint sum = 0; | |
for (uint i = 0; i < _power; ++i) | |
{ | |
(row, var r) = Math.DivRem(row, _base); | |
sum = (sum + r * coeff[(int)i]) % _base; | |
} | |
return sum; | |
} | |
} | |
} | |
private sealed class PrimeSeqence | |
{ | |
private static readonly uint[] _values = | |
[ | |
2, 3, 5, 7, 11, 13, 17, 19, | |
23, 29, 31, 37, 41, 43, 47, | |
]; | |
public uint Get(int index) | |
{ | |
return _values[index]; | |
} | |
} | |
private sealed class Plan | |
{ | |
private sealed class PseudoFactorGroup | |
{ | |
private uint _count = 0; | |
private uint _samples = 0; | |
private OrthogonalArray? _oa; | |
public uint Level { get; } | |
public uint SampleCount => _samples == 0 ? throw new InvalidOperationException() : _samples; | |
public PseudoFactorGroup(uint level) | |
{ | |
Debug.Assert(level >= 2); | |
Level = level; | |
} | |
public uint AddFactor() | |
{ | |
Debug.Assert(_samples == 0); | |
return _count++; | |
} | |
public void Generate() | |
{ | |
Debug.Assert(_samples == 0); | |
var level = Level; | |
var dof = (level - 1) * _count; | |
uint power = 0; | |
uint exp = 1; | |
while (exp <= dof) | |
{ | |
power += 1; | |
exp = checked(exp * level); | |
} | |
_samples = exp; | |
_oa = OrthogonalArray.Create(level, power, exp); | |
_oa.SetColumns(_count); | |
} | |
public uint Get(uint index, uint column) | |
{ | |
var oa = _oa ?? throw new InvalidOperationException(); | |
return oa.GetCell(index, column); | |
} | |
} | |
private readonly struct FactorMapping | |
{ | |
public int PseudoFactorIndex { get; init; } | |
public uint ColumnIndex { get; init; } | |
} | |
private readonly PrimeSeqence _primes = new(); | |
private readonly List<uint> _levels; | |
private readonly List<PseudoFactorGroup> _pseudoFactors = []; | |
private readonly Dictionary<uint, int> _pseudoFactorsByLevel = []; | |
private readonly List<List<FactorMapping>> _factorMappings = []; | |
private uint[] _pseudoFactorIndicesBuffer = null!; | |
private uint _samples = 0; | |
public uint SampleCount => _samples == 0 ? throw new InvalidOperationException() : _samples; | |
public Plan(List<uint> levels) | |
{ | |
_levels = levels; | |
} | |
public void Generate() | |
{ | |
SplitFactors(); | |
PreparePseudoFactors(); | |
} | |
private void SplitFactors() | |
{ | |
foreach (uint l in _levels) | |
{ | |
uint level = l; | |
var mapping = new List<FactorMapping>(); | |
_factorMappings.Add(mapping); | |
if (level > 1) | |
{ | |
FindAndAddPseudoFactor(ref level, 2, mapping); | |
} | |
if (level > 1) | |
{ | |
FindAndAddPseudoFactor(ref level, 3, mapping); | |
} | |
if (level > 1) | |
{ | |
FindAndAddPseudoFactor(ref level, 5, mapping); | |
} | |
if (level > 1) | |
{ | |
FindAndAddPseudoFactor(ref level, 7, mapping); | |
} | |
int pindex = 4; | |
while (level > 1) | |
{ | |
FindAndAddPseudoFactor(ref level, _primes.Get(pindex++), mapping); | |
} | |
} | |
} | |
[MethodImpl(MethodImplOptions.AggressiveInlining)] | |
private void FindAndAddPseudoFactor(ref uint level, uint pf, List<FactorMapping> mapping) | |
{ | |
while (level > 1) | |
{ | |
var (q, r) = Math.DivRem(level, pf); | |
if (r != 0) | |
{ | |
return; | |
} | |
if (!_pseudoFactorsByLevel.TryGetValue(pf, out var index)) | |
{ | |
index = _pseudoFactors.Count; | |
_pseudoFactorsByLevel[pf] = index; | |
_pseudoFactors.Add(new PseudoFactorGroup(pf)); | |
} | |
var pfo = _pseudoFactors[index]; | |
var col = pfo.AddFactor(); | |
mapping.Add(new() | |
{ | |
PseudoFactorIndex = index, | |
ColumnIndex = col, | |
}); | |
level = q; | |
} | |
} | |
private void PreparePseudoFactors() | |
{ | |
uint samples = 1; | |
foreach (var pf in _pseudoFactors) | |
{ | |
pf.Generate(); | |
samples = checked(samples * pf.SampleCount); | |
} | |
_samples = samples; | |
_pseudoFactorIndicesBuffer = new uint[_pseudoFactors.Count]; | |
} | |
public void GetSamples(uint index, Span<uint> result) | |
{ | |
if (result.Length != _levels.Count) | |
{ | |
throw new ArgumentException("length of result Span must be equal to factor count"); | |
} | |
for (int i = 0; i < _pseudoFactors.Count; ++i) | |
{ | |
var pf = _pseudoFactors[i]; | |
(index, _pseudoFactorIndicesBuffer[i]) = Math.DivRem(index, pf.SampleCount); | |
} | |
for (int i = 0; i < _levels.Count; ++i) | |
{ | |
var mapping = _factorMappings[i]; | |
uint multiplier = 1; | |
uint sum = 0; | |
foreach (var m in mapping) | |
{ | |
var pf = _pseudoFactors[m.PseudoFactorIndex]; | |
var pfv = pf.Get(_pseudoFactorIndicesBuffer[m.PseudoFactorIndex], m.ColumnIndex); | |
sum += pfv * multiplier; | |
multiplier *= pf.Level; | |
} | |
result[i] = sum; | |
} | |
} | |
} | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment