tiny.simd.CostDistribution

1pub const CostDistribution = struct {2    num_values: usize = 0,3    pointer_padding: [8 - @sizeOf(usize)]u8 = @splat(0),4    online_n: f64 = 0,5    values: [max_values]f64 = @splat(0),67    pub const max_values: usize = 14;89    const m1_index: usize = 0;10    const m2_index: usize = 1;11    const mean_index: usize = 2;12    const standard_deviation_index: usize = 3;13    const lower_index: usize = 4;14    const upper_index: usize = 5;1516    pub fn notify(self: *@This(), value: f64) void {17        if (value < 0) return;18        if (self.isOnline()) {19            self.onlineNotify(value);20            return;21        }22        std.debug.assert(self.num_values < max_values);23        self.values[self.num_values] = value;24        self.num_values += 1;25        if (self.num_values == max_values) self.warmUpOnline();26    }2728    pub fn estimateCost(self: *@This()) f64 {29        std.debug.assert(self.num_values != 0);30        if (!self.isOnline()) self.warmUpOnline();31        return self.mean();32    }3334    pub fn isOnline(self: *const @This()) bool {35        return self.online_n > 0;36    }3738    pub fn bufferedCount(self: *const @This()) usize {39        return self.num_values;40    }4142    pub fn onlineCount(self: *const @This()) f64 {43        return self.online_n;44    }4546    pub fn m1(self: *const @This()) f64 {47        return self.values[m1_index];48    }4950    pub fn m2(self: *const @This()) f64 {51        return self.values[m2_index];52    }5354    pub fn mean(self: *const @This()) f64 {55        return self.values[mean_index];56    }5758    pub fn standardDeviation(self: *const @This()) f64 {59        return self.values[standard_deviation_index];60    }6162    pub fn lower(self: *const @This()) f64 {63        return self.values[lower_index];64    }6566    pub fn upper(self: *const @This()) f64 {67        return self.values[upper_index];68    }6970    fn median(to_sort: []f64) f64 {71        std.debug.assert(to_sort.len >= 2);72        std.mem.sort(f64, to_sort, {}, std.sort.asc(f64));73        const middle = to_sort.len / 2;74        if (to_sort.len % 2 != 0) return to_sort[middle];75        return (to_sort[middle] + to_sort[middle - 1]) * 0.5;76    }7778    fn medianAbsoluteDeviation(values: []const f64, median_value: f64) f64 {79        var deviations: [max_values]f64 = undefined;80        for (values, deviations[0..values.len]) |value, *deviation| {81            deviation.* = @abs(value - median_value);82        }83        return median(deviations[0..values.len]);84    }8586    fn removeOutliers(self: *@This()) void {87        if (self.num_values < 3) return;88        std.debug.assert(self.num_values <= max_values);89        const active = self.values[0..self.num_values];90        const median_value = median(active);91        const mad = medianAbsoluteDeviation(active, median_value);92        if (mad == 0) {93            const skewness = (active[active.len - 1] - median_value) -94                (median_value - active[0]);95            const trim = @max(self.num_values / 2, 2);96            const left = @max(if (skewness < 0) trim * 3 / 4 else trim / 4, 1);97            self.num_values -= trim;98            std.debug.assert(self.num_values >= 1);99            std.mem.copyForwards(100                f64,101                self.values[0..self.num_values],102                self.values[left .. left + self.num_values],103            );104            return;105        }106107        const upper_bound = median_value + 5 * mad;108        const lower_bound = median_value - 5 * mad;109        var right = self.num_values - 1;110        while (self.values[right] > upper_bound) right -= 1;111        std.debug.assert(right >= self.num_values / 2);112        var left: usize = 0;113        while (left < right and self.values[left] < lower_bound) left += 1;114        std.debug.assert(left <= self.num_values / 2);115        self.num_values = right - left + 1;116        std.mem.copyForwards(117            f64,118            self.values[0..self.num_values],119            self.values[left .. left + self.num_values],120        );121    }122123    fn sampleMean(self: *const @This()) f64 {124        std.debug.assert(!self.isOnline());125        std.debug.assert(self.num_values != 0);126        std.debug.assert(self.num_values <= max_values);127        var sum: f64 = 0;128        for (self.values[0..self.num_values]) |value| sum += value;129        return sum / @as(f64, @floatFromInt(self.num_values));130    }131132    fn sampleVariance(self: *const @This(), sample_mean: f64) f64 {133        std.debug.assert(sample_mean >= 0);134        std.debug.assert(!self.isOnline());135        std.debug.assert(self.num_values != 0);136        std.debug.assert(self.num_values <= max_values);137        if (self.num_values == 1) return 0;138        var squared_sum: f64 = 0;139        for (self.values[0..self.num_values]) |value| {140            const difference = value - sample_mean;141            squared_sum += difference * difference;142        }143        return squared_sum / @as(f64, @floatFromInt(self.num_values - 1));144    }145146    fn onlineNotify(self: *@This(), unbounded_value: f64) void {147        const value = @min(@max(self.lower(), unbounded_value), self.upper());148        const n_minus_one = self.online_n;149        self.online_n += 1;150        const difference = value - self.values[m1_index];151        const difference_div_n = difference / self.online_n;152        self.values[m1_index] += difference_div_n;153        std.debug.assert(self.values[m1_index] >= self.lower());154        self.values[m2_index] += difference * n_minus_one * difference_div_n;155        const standard_deviation = @sqrt(self.values[m2_index] / @max(1, n_minus_one));156        self.values[mean_index] = self.values[m1_index] * 0.2 + self.values[mean_index] * 0.8;157        self.values[standard_deviation_index] = standard_deviation * 0.2 +158            self.values[standard_deviation_index] * 0.8;159        self.values[lower_index] = self.values[mean_index] -160            3.5 * self.values[standard_deviation_index];161        self.values[upper_index] = self.values[mean_index] +162            3.5 * self.values[standard_deviation_index];163    }164165    fn warmUpOnline(self: *@This()) void {166        self.removeOutliers();167        const sample_mean = self.sampleMean();168        const sample_variance = self.sampleVariance(sample_mean);169        var copy: [max_values]f64 = undefined;170        @memcpy(copy[0..self.num_values], self.values[0..self.num_values]);171        self.values[m1_index] = 0;172        self.values[m2_index] = 0;173        self.values[mean_index] = sample_mean;174        self.values[standard_deviation_index] = @sqrt(sample_variance);175        if (self.values[standard_deviation_index] == 0) {176            self.values[standard_deviation_index] = self.values[mean_index] / 2;177        }178        self.values[lower_index] = self.values[mean_index] -179            4 * self.values[standard_deviation_index];180        self.values[upper_index] = self.values[mean_index] +181            4 * self.values[standard_deviation_index];182        for (copy[0..self.num_values]) |value| self.onlineNotify(value);183        std.debug.assert(self.isOnline());184    }185};

1pub const CostDistribution = autotune.CostDistribution;