import { euclidean, chebyshev } from "../metrics/index.js";
import { MDS } from "../dimred/MDS.js";
import { Randomizer } from "../util/index.js";
import { BallTree } from "../knn/BallTree.js";
import { Matrix } from "../matrix/index.js";
import { neumair_sum } from "../numerical/index.js";
/**
*
*/
export class OAP {
constructor(X, depth_field_lag, step_size, depth_weight, d = 2, metric = euclidean, seed = 1212) {
this._X = X;
this._d = d;
[this._N, this._D] = X.shape;
this._depth_field_lag = depth_field_lag;
this._step_size = step_size;
this._depth_weight = depth_weight;
this._J = 3;
this._max_iter = 1;
this._metric = metric;
this._seed = seed;
this._randomizer = new Randomizer(seed);
}
_data_depth(technique = "chebyshev") {
const X = this._X;
const N = this._N;
const h = new Float32Array(N);
let deepest_point = 0;
if (technique === "mdb") {
h.fill(1);
/*
// Modified Band Depth
// https://www.tandfonline.com/doi/pdf/10.1198/jasa.2009.0108?casa_token=E1Uzntgs-5AAAAAA:Eo8mUpJDhpLQ5RHBkCB3Mdz0tbGM3Q0v78bwyCIAv7-peLGwfG3TcXLqShIaYuJLEqKc7GvaKlgvUg
const randomizer = this._randomizer;
const h = new Float32Array(this._N);
const J = this._J;
const N = this._N;
const D = this._D;
const X = this._X;
const one_div_by_n_choose_j = 1;
for (let row = 0; row < N; ++row) {
const x = X.row(row);
const B_min = new Float32Array(D).fill(Infinity);
const B_max = new Float32Array(D).fill(-Infinity);
let r = Math.floor(randomizer.random * N);
for (let i = 0; i < J; ++i) {
const x_j = X.row(r);
for (let d = 0; d < D; ++d) {
const x_jd = x_j[d]
B_min[d] = Math.min(B_min[d], x_jd);
B_max[d] = Math.max(B_max[d], x_jd);
}
r += Math.floor(randomizer.random * (N - 1));
r = r % N;
}
for (let d = 0; d < D; ++d) {
const x_d = x[d];
if (x_d >= B_min[d] && x_d <= B_max[d]) {
++h[row]
}
}
}
this._h = h;*/
} else if (technique === "chebyshev") {
// L∞ Depth
// https://arxiv.org/pdf/1506.01332.pdf
for (let i = 0; i < N; ++i) {
let x = X.row(i);
let sum = 0;
for (let j = 0; j < N; ++j) {
if (i !== j) {
sum += chebyshev(x, X.row(j));
}
}
h[i] = 1 / (1 + sum / N);
if (h[deepest_point] < h[i]) {
deepest_point = i;
}
}
}
this._h = h;
this._deepest_point = deepest_point;
}
init() {
this._iter = 0;
// init with MDS
const init_MDS = new MDS(this._X, this._d, this._metric);
//console.log(init_MDS)
this._Y = init_MDS.transform();
// try häääh?
this._X_distances = init_MDS._d_X;
/*let max = -Infinity
init_MDS._d_X._data.forEach(dx => max = Math.max(dx, max));
this._X_distances = init_MDS._d_X.divide(max);*/
// end try hääääh?
// compute order statistics
this._data_depth();
this._M = this._monotonic_field(this._Y);
//
return this;
}
set depth_field_lag(value) {
this._depth_field_lag = value;
}
get depth_field_lag() {
return this._depth_field_lag;
}
set step_size(value) {
this._step_size = value;
}
get step_size() {
return this._step_size;
}
set depth_weight(value) {
this._depth_weight = value;
}
get depth_weight() {
return this._depth_weight;
}
transform(iterations = this._max_iter) {
for (let i = 0; i < iterations; ++i) {
this.next();
}
return this._Y;
}
*transform_iter() {
while (true) {
this.next();
yield this._Y;
}
}
_monotonic_field(Y) {
const h = this._h;
const Y_ = this._Y_;
Y, h, Y_;
const nn = new BallTree();
nn.add(Y.to2dArray);
const N = 5;
let M = (x) => {
let neighbors = nn.search(x, N).toArray();
let d_sum = 0; //neighbors.reduce((a, b) => a + b.value, 0);
let m = 0;
for (let i = 0; i < N; ++i) {
d_sum += neighbors[i].value;
m += h[neighbors[i].element.index] * neighbors[i].value;
}
//console.log(m, d_sum)
m /= d_sum;
return m;
};
return M;
}
next() {
const iter = ++this._iter;
const l = this._depth_field_lag;
const step_size = this._step_size;
const w = this._depth_weight;
const N = this._N;
const dim = this._d;
const d_X = this._X_distances;
const h = this._h;
let Y = this._Y;
if (iter % l === 1) {
// compute monotonic field
this._Y_ = this._Y.clone();
this._M = this._monotonic_field(Y);
}
const M = this._M;
// perform gradient step
// MDS stress step
/*for (let i = 0; i < N; ++i) {
const d_x = d_X.row(i);
const y_i = Y.row(i)
const delta_mds_stress = new Float32Array(dim);
for (let j = 0; j < N; ++j) {
if (i !== j) {
const y_j = Y.row(j)
const d_y = metric(y_i, y_j);
const d_x_j = d_x[j] === 0 ? 1e-2 : d_x[j]
const mult = 1 - (d_x_j / d_y)
for (let d = 0; d < dim; ++d) {
delta_mds_stress[d] += (mult * (y_i[d] - y_j[d]));
}
}
}
for (let d = 0; d < dim; ++d) {
Y.set_entry(i, d, Y.entry(i, d) - step_size * delta_mds_stress[d] / N)
}
}*/
// MDS stress step
const d_Y = new Matrix();
d_Y.shape = [
N,
N,
(i, j) => {
return i < j ? euclidean(Y.row(i), Y.row(j)) : d_Y.entry(j, i);
},
];
const ratio = new Matrix(); //d_X.divide(d_Y).mult(-1);
ratio.shape = [
N,
N,
(i, j) => {
if (i === j) return 1e-8;
return i < j ? -d_X.entry(i, j) / d_Y.entry(i, j) : ratio.entry(j, i);
},
];
for (let i = 0; i < N; ++i) {
ratio.sub_entry(i, i, neumair_sum(ratio.row(i)));
}
const mds_Y = ratio.dot(Y).divide(N);
// Data depth step
const diff_Y = new Matrix(N, dim, (i, j) => mds_Y.entry(i, j) - Y.entry(i, j));
for (let i = 0; i < N; ++i) {
const m = M(Y.row(i));
const dm = M(mds_Y.row(i));
const h_i = h[i];
for (let d = 0; d < dim; ++d) {
Y.add_entry(i, d, step_size * (diff_Y.entry(i, d) + w * 2 * (m - h_i) * dm));
}
}
this._Y = Y;
return this._Y;
}
get projection() {
return this._Y;
}
}