diff --git a/js/cnvpytor/CombinedCaller.js b/js/cnvpytor/CombinedCaller.js index 996da916f..f9d7d519d 100644 --- a/js/cnvpytor/CombinedCaller.js +++ b/js/cnvpytor/CombinedCaller.js @@ -1,25 +1,33 @@ import t_dist from './t_dist.js' -import g_utils from './GeneralUtil.js' - - -class CombinedCaller{ - constructor(wigFeatures, binSize) { - this.wigFeatures = wigFeatures - this.binSize = binSize +import gUtils from './GeneralUtil.js' +import baseCNVpytorVCF from './baseCNVpytorVCF.js' + + +class CombinedCaller extends baseCNVpytorVCF { + /** + * Creates an instance of CombinedCaller. + * + * @param {Array} wigFeatures - An array of arrays containing wig formatted data for each chromosome and bin. + * @param {number} binSize - The size of the bins used in the wig data. + * @param {string} refGenome - GC content data indexed by chromosome and bin. + */ + constructor(wigFeatures, binSize, refGenome) { + super(wigFeatures, binSize, refGenome) } - get_fit(){ - var fit_info = new g_utils.GetFit(this.wigFeatures) - var [globalMean, globalStd] = fit_info.fit_data() - return {globalMean:globalMean, globalStd:globalStd} - - } async call_2d(omin=null, mcount=null, event_type="both", max_distance=0.1, baf_threshold=0, max_copy_number=10, min_cell_fraction=0.0){ - let fit_obj = this.get_fit() - this.globalMean = fit_obj.globalMean - this.globalStd = fit_obj.globalStd + // let fit_obj = this.get_fit_v2() + // this.globalMean = fit_obj.globalMean + // this.globalStd = fit_obj.globalStd + // this.getGcCorrectionSignal(this.globalMean) + + // applying gc correction + await this.apply_gcCorrection() + + let binScoreField = this.gcFlag ? "gcCorrectedBinScore": "binScore" ; + let overlap_min = omin==null ? 0.05 * this.binSize / 3e9: omin ; let min_count = mcount == null ? parseInt(this.binSize / 10000) : mcount ; @@ -41,11 +49,14 @@ class CombinedCaller{ if (bin.hets_count > 4 ){ if( bin.dp_count > min_count ){ - segments.push([bin_idx]) - levels.push(bin.binScore) - likelihoods.push(bin.likelihood_score) - delete bin.likelihood_score - + if(bin[binScoreField]){ + segments.push([bin_idx]) + levels.push(bin[binScoreField]) + likelihoods.push(bin.likelihood_score) + delete bin.likelihood_score + + } + } } }) @@ -234,7 +245,7 @@ class CombinedCaller{ if(points == 0){ points = 1 } - let x = g_utils.linspace(min_cell_fraction, 1, points) + let x = gUtils.linspace(min_cell_fraction, 1, points) let master_lh = {} let germline_lh = {} for(let cn=10; cn > -1; cn--){ @@ -330,12 +341,17 @@ class CombinedCaller{ } var rawbinScore = this.formatDataStructure(this.wigFeatures, 'binScore', this.globalMean) + + let gcCorrectedBinScore = []; + if (this.gcFlag) { + gcCorrectedBinScore = this.formatDataStructure(this.wigFeatures, 'gcCorrectedBinScore', this.globalMean); + } var callScore = this.formatDataStructure(this.wigFeatures, 'segment_score', this.globalMean) - return {binScore: rawbinScore, segment_score: callScore} + return {binScore: rawbinScore, gcCorrectedBinScore: gcCorrectedBinScore, segmentScore: callScore} } - + formatDataStructure(wigFeatures, feature_column, scaling_factor = 1) { const results = [] for (const [chr, wig] of Object.entries(wigFeatures)) { @@ -351,7 +367,7 @@ class CombinedCaller{ return results } - + /* formatDataStructure_BAF(feature_column, scaling_factor = -1) { const baf1 = [] const baf2 = [] @@ -375,7 +391,7 @@ class CombinedCaller{ return [baf1, baf2] - } + }*/ } function arrayMax(arr) { @@ -424,6 +440,7 @@ function normal_overlap_approx(m1, s1, m2, s2){ * @returns {float} - Overlap area. */ function likelihood_overlap(likelihood_1, likelihood_2){ + // console.log(likelihood_1, likelihood_2) let sum; try{ sum = likelihood_1.reduce((accumulator, currentValue, currentIndex) => {return accumulator + Math.min(currentValue, likelihood_2[currentIndex])}); diff --git a/js/cnvpytor/GeneralUtil.js b/js/cnvpytor/GeneralUtil.js index 953d69ed3..30e44cadd 100644 --- a/js/cnvpytor/GeneralUtil.js +++ b/js/cnvpytor/GeneralUtil.js @@ -1,99 +1,141 @@ class GetFit { - constructor(allBins) { - this.allBins = allBins - } - getValues() { - const bins = Object.values(this.allBins).reduce( - (binResult, bin) => { return binResult.concat(bin.filter(a => a.binScore > 0).map(a => a.binScore)) }, []) - return bins - } - getMean(data) { - return (data.reduce(function (a, b) { return a + b; }) / data.length); - } - fit_data() { - let rd_list = this.getValues() - let distParmas = getDistParams(rd_list) - return distParmas - } - - histogram(data, bins) { - const step = bins[1] - bins[0]; - const hist_bins = []; - - data.forEach((value, index) => { - bins.forEach((bin_value, bin_index) => { - if (!hist_bins[bin_value]) { - hist_bins[bin_value] = { count: 0 }; - } - if (bin_value <= value && value < bin_value + step) { - hist_bins[bin_value].count++; - return false; - } - }); - }); - const dist_p = [] - hist_bins.forEach((bin, index) => { dist_p.push(bin.count); }); - return dist_p - } + /** + * Creates an instance of GetFit. + * @param {Object} allBins - An object containing all the bins with their respective data. + */ + constructor(allBins) { + this.allBins = allBins // Stores all bins data + } + + /** + * Extracts bin scores greater than zero from all bins. + * @returns {Array} An array of bin scores. + */ + + getValues() { + const bins = Object.values(this.allBins).reduce( + (binResult, bin) => { return binResult.concat(bin.filter(a => a.binScore > 0).map(a => a.binScore)) }, []) + return bins + } + + /** + * Calculates the mean of the given data. + * @param {Array} data - The data array to calculate the mean from. + * @returns {number} The mean value of the data. + */ + getMean(data) { + return (data.reduce(function (a, b) { return a + b; }) / data.length); + } + + fit_data() { + let rd_list = this.getValues() + let distParmas = getDistParams(rd_list) + return distParmas + } + + + histogram(data, bins) { + const step = bins[1] - bins[0]; + const hist_bins = []; + + data.forEach((value, index) => { + bins.forEach((bin_value, bin_index) => { + if (!hist_bins[bin_value]) { + hist_bins[bin_value] = { count: 0 }; + } + if (bin_value <= value && value < bin_value + step) { + hist_bins[bin_value].count++; + return false; + } + }); + }); + const dist_p = [] + hist_bins.forEach((bin, index) => { dist_p.push(bin.count); }); + return dist_p + } } + function range_function(start, stop, step) { - const data_array = Array(Math.ceil((stop - start) / step)) - .fill(start) - .map((x, y) => x + y * step); - return data_array; + const data_array = Array(Math.ceil((stop - start) / step)) + .fill(start) + .map((x, y) => x + y * step); + return data_array; } + function Gaussian([a, x0, sigma]) { - return x => - (a * Math.exp(-Math.pow(x - x0, 2) / (2 * Math.pow(sigma, 2)))) / (Math.sqrt(2 * Math.PI) * sigma); + return x => + (a * Math.exp(-Math.pow(x - x0, 2) / (2 * Math.pow(sigma, 2)))) / (Math.sqrt(2 * Math.PI) * sigma); } function filterOutliers(someArray) { - if (someArray.length < 4) - return someArray; + if (someArray.length < 4) + return someArray; - let values, q1, q3, iqr, maxValue, minValue; + let values, q1, q3, iqr, maxValue, minValue; - values = someArray.slice().sort((a, b) => a - b); //copy array fast and sort + values = someArray.slice().sort((a, b) => a - b); //copy array fast and sort - if ((values.length / 4) % 1 === 0) { //find quartiles - q1 = 1 / 2 * (values[(values.length / 4)] + values[(values.length / 4) + 1]); - q3 = 1 / 2 * (values[(values.length * (3 / 4))] + values[(values.length * (3 / 4)) + 1]); - } else { - q1 = values[Math.floor(values.length / 4 + 1)]; - q3 = values[Math.ceil(values.length * (3 / 4) + 1)]; - } + if ((values.length / 4) % 1 === 0) { //find quartiles + q1 = 1 / 2 * (values[(values.length / 4)] + values[(values.length / 4) + 1]); + q3 = 1 / 2 * (values[(values.length * (3 / 4))] + values[(values.length * (3 / 4)) + 1]); + } else { + q1 = values[Math.floor(values.length / 4 + 1)]; + q3 = values[Math.ceil(values.length * (3 / 4) + 1)]; + } - iqr = q3 - q1; - maxValue = q3 + iqr * 1.5; - minValue = q1 - iqr * 1.5; + iqr = q3 - q1; + maxValue = q3 + iqr * 1.5; + minValue = q1 - iqr * 1.5; - return values.filter((x) => (x >= minValue) && (x <= maxValue)); + return values.filter((x) => (x >= minValue) && (x <= maxValue)); } function getDistParams(bins) { - let filteredBins = filterOutliers(bins) - const n = filteredBins.length - const mean = filteredBins.reduce((a, b) => a + b) / n - const std = Math.sqrt(filteredBins.map(x => Math.pow(x - mean, 2)).reduce((a, b) => a + b) / n) - return [mean, std] + let filteredBins = filterOutliers(bins) + const n = filteredBins.length + const mean = filteredBins.reduce((a, b) => a + b) / n + const std = Math.sqrt(filteredBins.map(x => Math.pow(x - mean, 2)).reduce((a, b) => a + b) / n) + return [mean, std] } function linspace(a, b, n) { - if (typeof n === "undefined") n = Math.max(Math.round(b - a) + 1, 1); - if (n < 2) { - return n === 1 ? [a] : []; - } - var ret = Array(n); - n--; - for (let i = n; i >= 0; i--) { - ret[i] = (i * b + (n - i) * a) / n; - } - return ret; + if (typeof n === "undefined") n = Math.max(Math.round(b - a) + 1, 1); + if (n < 2) { + return n === 1 ? [a] : []; + } + var ret = Array(n); + n--; + for (let i = n; i >= 0; i--) { + ret[i] = (i * b + (n - i) * a) / n; + } + return ret; } -export default { range_function, getDistParams, linspace, GetFit}; +export function histogram2d(data1, data2, binsX, binsY) { + // Calculate bin sizes + const minX = math.min(data1); + const maxX = math.max(data1); + const minY = math.min(data2); + const maxY = math.max(data2); + const binSizeX = (maxX - minX) / binsX; + const binSizeY = (maxY - minY) / binsY; + + // Create the histogram array + const histogram = math.zeros(binsX, binsY); + + // Populate the histogram + for (let i = 0; i < data1.length; i++) { + const xBin = Math.floor((data1[i] - minX) / binSizeX); + const yBin = Math.floor((data2[i] - minY) / binSizeY); + histogram.set([xBin, yBin], histogram.get([xBin, yBin]) + 1); + } + + return histogram; + } + +export default { range_function, getDistParams, linspace, GetFit, filterOutliers }; diff --git a/js/cnvpytor/HDF5IndexedReader.js b/js/cnvpytor/HDF5IndexedReader.js index 4fd31d066..5b7b374ec 100644 --- a/js/cnvpytor/HDF5IndexedReader.js +++ b/js/cnvpytor/HDF5IndexedReader.js @@ -36,6 +36,8 @@ class HDF5Reader { this.h5_file = h5_file; this.bin_size = bin_size; this.h5_obj = undefined + this.pytorKeys = []; + this.availableBins = []; } async fetch(){ @@ -59,100 +61,76 @@ class HDF5Reader { return h5_obj.keys } - async get_rd_signal(bin_size=this.bin_size){ + async get_rd_signal(bin_size = this.bin_size, chrom=undefined){ + // Fetch the pytor file and get keys + const h5Obj = await this.fetch(); + this.pytorKeys = h5Obj.keys; - let h5_obj = await this.fetch(); + // get available bin sizes + const signalBin = new ParseSignals(this.pytorKeys); + this.availableBins = signalBin.getAllBins(); + + // check if the user provided bin is available, else set the last bin_size + if(! this.availableBins.includes(bin_size)){ + bin_size = this.availableBins.at(-1); + } - this.h5_obj = h5_obj - this.pytor_keys = h5_obj.keys + // get rd chromosomes and rd stat + const rdChromosomes = await this.getChromosomes(chrom); - // get available bin sizes - let signal_bin = new ParseSignals(this.pytor_keys); - let rd_bins = signal_bin.get_rd_bins() - let snp_bins = signal_bin.get_snp_bins() + let rd_stat = await this.rd_stat(bin_size) - // merge the bins get from two type of signals - this.available_bins = [...new Set(rd_bins, snp_bins)] + // prepare wig formatted file for all chromosome + const wigFeatures = await this.getWigFeatures(rdChromosomes, bin_size, rd_stat); - // check if the user provided bin is available, else set the last bin_size - if(! this.available_bins.includes(bin_size)){ - bin_size = this.available_bins.at(-1); - } + this.setCallers(wigFeatures); + return { [bin_size]: wigFeatures }; + } - const chr_ds = await h5_obj.get("rd_chromosomes") - const type = await chr_ds.dtype - const t0 = Date.now() - let rd_chromosomes = await chr_ds.value - let rd_flag = "" + async getWigFeatures(rdChromosomes, binSize, rdStat) { + const wigFeatures = { + RD_Raw: [], + RD_Raw_gc_coor: [], + ReadDepth: [], + "2D": [], + BAF1: [], + BAF2: [] + }; - // get rd stat - let rd_stat = await this.rd_stat(bin_size) - - var wigFeatures = [] - var wigFeatures_gc = [] - var wigFeatures_rd_call_meanshift = [] - var wigFeatures_rd_call_combined = [] - var wigFeatures_baf1 = [] - var wigFeatures_baf2 = [] - - for (let chrom of rd_chromosomes) { - let signal_name_obj = new SignalNames(chrom, bin_size) - - // for normal rd signal - // var signal_rd = `his_rd_p_${chrom}_${bin_size}${rd_flag}` - var signal_rd = signal_name_obj.signals['raw_RD'] - let chr_wig = await this.get_chr_signal( chrom, bin_size, signal_rd, rd_stat) - wigFeatures = wigFeatures.concat(chr_wig) - - // rd gc corrected - // var signal_rd_gc = `his_rd_p_${chrom}_${bin_size}_GC` - var signal_rd_gc = signal_name_obj.signals['gc_RD'] - let chr_wig_gc = await this.get_chr_signal(chrom, bin_size, signal_rd_gc, rd_stat) - wigFeatures_gc = wigFeatures_gc.concat(chr_wig_gc) + for (const chrom of rdChromosomes) { + const signalNameObj = new SignalNames(chrom, binSize); - // rd call MeanShift - - // let signal_rd_call = `his_rd_p_${chrom}_${bin_size}_partition_GC_merge` - let signal_rd_call = signal_name_obj.signals['gc_partition'] - let chr_wig_rd_call_meanshift = await this.get_chr_signal(chrom, bin_size, signal_rd_call, rd_stat) - wigFeatures_rd_call_meanshift = wigFeatures_rd_call_meanshift.concat(chr_wig_rd_call_meanshift) - - let chr_wig_rd_call = await this.rd_call_combined(chrom, bin_size, rd_stat, signal_name_obj) - wigFeatures_rd_call_combined = wigFeatures_rd_call_combined.concat(chr_wig_rd_call) + wigFeatures.RD_Raw.push(...await this.get_chr_signal(chrom, binSize, signalNameObj.signals.raw_RD, rdStat)); + wigFeatures.RD_Raw_gc_coor.push(...await this.get_chr_signal(chrom, binSize, signalNameObj.signals.gc_RD, rdStat)); + wigFeatures.ReadDepth.push(...await this.get_chr_signal(chrom, binSize, signalNameObj.signals.gc_partition, rdStat)); - // baf likelihood - // let signal_baf_1 = `snp_likelihood_${chrom}_${bin_size}_mask` - // let signal_baf_1 = signal_name_obj.signals['baf'] - // let chr_wig_bafs = await this.get_baf_signals(chrom, bin_size, signal_baf_1) - - // let signal_baf_1 = `snp_i1_${chrom}_${bin_size}_mask` - let signal_baf_1 = signal_name_obj.signals['baf_i1'] - let chr_wig_bafs = await this.get_baf_signals_v2(chrom, bin_size, signal_baf_1) + wigFeatures["2D"].push(...await this.rd_call_combined(chrom, binSize, rdStat, signalNameObj)); + const [baf1, baf2] = await this.getBafSignals(chrom, binSize, signalNameObj.signals.baf_i1); + wigFeatures.BAF1.push(...baf1); + wigFeatures.BAF2.push(...baf2); + } - wigFeatures_baf1 = wigFeatures_baf1.concat(chr_wig_bafs[0]) - wigFeatures_baf2 = wigFeatures_baf2.concat(chr_wig_bafs[1]) + return wigFeatures; + } + + async getChromosomes(refChroms) { + // return chromosome names if they exists in the rd_chromosomes + const rdChroms_obj = await this.h5_obj.get("rd_chromosomes"); + const rdChroms = await rdChroms_obj.value + if(!refChroms){ + return rdChroms + }else{ + let refChromsSet = new Set(refChroms) + return rdChroms.filter(item => refChromsSet.has(item)); } - this.callers = [] - if (wigFeatures_rd_call_combined.length != 0){ - this.callers.push('ReadDepth') - } - if (wigFeatures_rd_call_combined.length != 0){ - this.callers.push('2D') - } + } - var obj = {} - var signal_obj = { - "RD_Raw": wigFeatures, - "RD_Raw_gc_coor" : wigFeatures_gc, - "ReadDepth": wigFeatures_rd_call_meanshift, - "2D": wigFeatures_rd_call_combined, - "BAF1": wigFeatures_baf1, - "BAF2": wigFeatures_baf2 - } - obj[bin_size] = signal_obj - return obj + setCallers(wigFeatures) { + this.callers = []; + if (wigFeatures.ReadDepth.length) this.callers.push('ReadDepth'); + if (wigFeatures["2D"].length) this.callers.push('2D'); } decode_segments(segments_arr){ @@ -174,9 +152,8 @@ class HDF5Reader { let chr_wig = []; let segments - // let mosaic_call_segments = `his_rd_p_${chrom}_${bin_size}_partition_GC_mosaic_segments_2d` let mosaic_call_segments = signal_name_obj.signals['Mosaic_segments'] - if (this.pytor_keys.includes(mosaic_call_segments)){ + if (this.pytorKeys.includes(mosaic_call_segments)){ const chrom_dataset = await this.h5_obj.get(mosaic_call_segments) const t0 = Date.now() let chrom_data = await chrom_dataset.value @@ -184,9 +161,8 @@ class HDF5Reader { } - // let mosaic_calls = `his_rd_p_${chrom}_${bin_size}_partition_GC_mosaic_call_2d` let mosaic_calls = signal_name_obj.signals['Mosaic_calls'] - if (this.pytor_keys.includes(mosaic_calls)){ + if (this.pytorKeys.includes(mosaic_calls)){ const segments_call_dataset = await this.h5_obj.get(mosaic_calls) let segments_call = await segments_call_dataset.to_array() //create_nested_array(value, shape) segments.forEach((ind_segment, segment_idx) => { @@ -209,12 +185,9 @@ class HDF5Reader { let rd_stat_signal = `rd_stat_${bin_size}_auto` let rd_stat; - - if (this.pytor_keys.includes(rd_stat_signal)){ + if (this.pytorKeys.includes(rd_stat_signal)){ const rd_stat_dataset = await this.h5_obj.get(rd_stat_signal) - const t0 = Date.now() rd_stat = await rd_stat_dataset.value - //console.log(`rd_stat_signal ${rd_stat_signal} ${Date.now() - t0}`) } return rd_stat } @@ -224,11 +197,10 @@ class HDF5Reader { /* return a list of dictionary for a chromosome */ let chr_wig = []; - if (this.pytor_keys.includes(signal_name)){ + if (this.pytorKeys.includes(signal_name)){ const chrom_dataset = await this.h5_obj.get(signal_name) let chrom_data = await chrom_dataset.value - //console.log(`chr_signal ${signal_name} ${Date.now() - t0}`) chrom_data.forEach((bin_value, bin_idx) => { chr_wig.push({chr:chrom, start: bin_idx*bin_size, end: (bin_idx+1) * bin_size, value: (bin_value/rd_stat[4]) *2}) }); @@ -236,89 +208,68 @@ class HDF5Reader { return chr_wig } - async get_baf_signals(chrom, bin_size, signal_name, scaling_factor = -1){ - /* return two list of dictionary*/ - let chr_wig_1 = []; - let chr_wig_2 = []; - if (this.pytor_keys.includes(signal_name)){ - let chrom_dataset = await this.h5_obj.get(signal_name) - let chrom_data = await chrom_dataset.to_array() //create_nested_array(value, shape) - chrom_data.forEach((bin_value, bin_idx) => { - let max_value = Math.max(...bin_value); - const res = bin_value.indexOf(max_value); - let lh = Math.max(res / 200, 1 - res / 200); - chr_wig_1.push({chr:chrom, start: bin_idx*bin_size, end: (bin_idx+1) * bin_size, value: scaling_factor * lh}) - if(lh != 0.5){ - chr_wig_2.push({chr:chrom, start: bin_idx*bin_size, end: (bin_idx+1) * bin_size, value: scaling_factor *(1-lh)}) - } - }); - } - return [chr_wig_1, chr_wig_2] - } - async get_baf_signals_v2(chrom, bin_size, signal_name, scaling_factor = -1){ + async getBafSignals(chrom, binSize, signalName, scalingFactor = -1) { + const chrWig1 = []; + const chrWig2 = []; - /* return two list of dictionary*/ - let chr_wig_1 = []; - let chr_wig_2 = []; - if (this.pytor_keys.includes(signal_name)){ - let chrom_dataset = await this.h5_obj.get(signal_name) - let chrom_data = await chrom_dataset.to_array() //create_nested_array(value, shape) - chrom_data.forEach((lh, bin_idx) => { - if (!isNaN(lh)){ - chr_wig_1.push({chr:chrom, start: bin_idx*bin_size, end: (bin_idx+1) * bin_size, value: scaling_factor * ( 0.5 - lh )}) - if(lh != 0.5){ - chr_wig_2.push({chr:chrom, start: bin_idx*bin_size, end: (bin_idx+1) * bin_size, value: scaling_factor * ( 0.5 + lh )}) + if (this.pytorKeys.includes(signalName)) { + const chromDataset = await this.h5_obj.get(signalName); + const chromData = await chromDataset.to_array(); + + chromData.forEach((lh, binIdx) => { + if (!isNaN(lh)) { + chrWig1.push({ + chr: chrom, + start: binIdx * binSize, + end: (binIdx + 1) * binSize, + value: scalingFactor * (0.5 - lh) + }); + if (lh !== 0.5) { + chrWig2.push({ + chr: chrom, + start: binIdx * binSize, + end: (binIdx + 1) * binSize, + value: scalingFactor * (0.5 + lh) + }); } } }); } - //console.log(chrom, chr_wig_1, chr_wig_2) - return [chr_wig_1, chr_wig_2] - + return [chrWig1, chrWig2]; } -} -class ParseSignals{ +} +class ParseSignals { /** - * Parse a signal names - * - * @param {*} signals - List of keys in pytor files + * @param {string[]} signals - List of keys in pytor files. */ - constructor(signals){ - this.signals = signals + constructor(signals) { + this.signals = signals; } - /** - * - * @returns - return list of rd bins - */ - get_rd_bins(){ - let rd_keys = []; - this.signals.forEach( val => { - let match = val.match(/^his_rd_p_(.*)_(\d+)$/); - if(match){ - rd_keys.push({chr:match[1], bin_size:match[2]}) - }}); - const rd_bins = [...new Set(rd_keys.map(item => Number(item.bin_size)))]; - return rd_bins + getAllBins() { + const rdBins = this.getRdBins(); + const snpBins = this.getSnpBins(); + return [...new Set([...rdBins, ...snpBins])].sort((a, b) => a - b);; } - /** - * - * @returns - return list of snp bins - */ - get_snp_bins(){ - - let slected_signals = []; - this.signals.forEach( val => { - let match = val.match(/^snp_likelihood_(.*)_(\d+)_mask$/); - if(match){ - slected_signals.push({chr:match[1], bin_size:match[2]}) - }}); - const bins = [...new Set(slected_signals.map(item => Number(item.bin_size)))]; - return bins + getRdBins() { + return this.extractBins(/^his_rd_p_(.*)_(\d+)$/); + } + + getSnpBins() { + return this.extractBins(/^snp_likelihood_(.*)_(\d+)_mask$/); + } + + extractBins(regex) { + return [...new Set( + this.signals + .map(val => val.match(regex)) + .filter(match => match !== null) + .map(match => Number(match[2])) + )]; } } diff --git a/js/cnvpytor/MeanShiftUtil.js b/js/cnvpytor/MeanShiftUtil.js index 20fb877bd..0421cfc9f 100644 --- a/js/cnvpytor/MeanShiftUtil.js +++ b/js/cnvpytor/MeanShiftUtil.js @@ -1,10 +1,438 @@ -import g_utils from './GeneralUtil.js' +// import g_utils from './GeneralUtil.js' import t_dist from './t_dist.js' +import baseCNVpytorVCF from './baseCNVpytorVCF.js' // TODO -- remove this hardcoded value const genome_size = 2871000000; + +class MeanShiftCaller extends baseCNVpytorVCF{ + /** + * Creates an instance of CombinedCaller. + * + * @param {Array} wigFeatures - An array of arrays containing wig formatted data for each chromosome and bin. + * @param {number} binSize - The size of the bins used in the wig data. + * @param {string} refGenome - reference genome name + */ + constructor(wigFeatures, binSize, refGenome) { + super(wigFeatures, binSize, refGenome) + this.binBands = [2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 64, 80, 96, 112, 128] + } + + async callMeanshift(repeats = 3){ + // applying gc correction + await this.apply_gcCorrection() + + let partitionLevels = this.partition() + // console.log("Partition: ", partitionLevels) + + // console.log("WigFeatures: ", this.wigFeatures) + + // let cnvs = this.cnvCalling(partitionLevels) + // console.log("cnvs: ", cnvs) + + Object.entries(this.wigFeatures).forEach(([chr, chrRD]) => { + chrRD.forEach((bin, index) => { + if (partitionLevels[chr]){ + bin.partitionLevel = parseInt(partitionLevels[chr][index]) + } + // bin.cnvCall = parseInt(caller_array[0][chr][index]) + }); + }) + + var rawbinScore = this.formatDataStructure('binScore', this.globalMean) + var gcCorrectedBinScore = this.formatDataStructure('gcCorrectedBinScore', this.globalMean) + var partitionBinScore = this.formatDataStructure('partitionLevel', this.globalMean) + // var cnvLevels = this.formatDataStructure('cnvCall', this.globalMean) + + const fetchedData = {binScore: rawbinScore, gcCorrectedBinScore: gcCorrectedBinScore, segmentsCNV: partitionBinScore} + return fetchedData + + } + + getRDSignalBandWidth(data_array) { + const threshold = this.globalMean / 4; + + // The values are reversed and squared as they will be used to calculate a gradient function. + // This optimization is done to speed up the calculations. + + const constantValue = 4 / this.globalStd ** 2; + return data_array.map(value => { + return value > threshold ? this.globalMean / (this.globalStd ** 2 * value) : constantValue; + }); + } + + partition(repeats = 3){ + + // sort the dictionary based on chromosome names; + let sortedDictionary = {}; + Object.keys(this.wigFeatures).sort((a, b) => a.localeCompare(b, undefined, {numeric: true})).forEach(key => { + sortedDictionary[key] = this.wigFeatures[key]; + }); + + let binScoreField = this.gcFlag ? "gcCorrectedBinScore": "binScore" ; + + var chrLevels = {} + // Object.entries(this.wigFeatures).forEach(([chr, chr_rd]) => { + for (const [chr, chrWig] of Object.entries(sortedDictionary)) { + + // console.log("chr: ", chr, chrWig.length, chrWig) + + // boolean array; Initiate with all false values + var masked = new Array(chrWig.length).fill(false) + + // set the level; score from either RAW or GC corrected bin score + var levels = chrWig.map((item, index) => !masked[index] ? item[binScoreField] : undefined); + // console.log("Levels: ", chr, levels) + // var levels = chrWig.map((item, index) => !masked[index] ? item : undefined); + + + this.binBands.forEach((bin_band, bin_band_index) => { + + // console.log("BinBand: ", bin_band) + + // not masked levels at current bin + // get boolean values + var not_masked = masked.map((value, index) => { return !value; }) + // console.log(not_masked) + + // not masked level at current bin + // Object.entries(chrWig).forEach(([k, v]) => { nm_levels.push(v.binScore) }) + // var nm_levels = Object.keys(chrWig).map(k => chrWig[k].binScore); + // var nm_levels = levels + + // not mask level + var nm_levels = levels.filter((_, index) => !masked[index]); + // console.log("nm_levels : ( Bin Band: " , bin_band, " )", nm_levels) + // const nm_levels = Object.keys(chrWig).map(k => chrWig[k].binScore).filter(score => !isNaN(score)); + + // console.log(nm_levels) + + // set the mask border + var mask_borders = [0] + var count = 0 + + // the masked array was declared previously + masked.forEach(item => { + if (item) { + if (count > 0) { + mask_borders.push(mask_borders[mask_borders.length - 1] + count - 1); + count = 0; + } + } else if (!item) { count++; } + }); + // console.log(mask_borders) + // console.log("Mask Borders: ", mask_borders) + mask_borders.shift() + + // repeating steps + for (let step = 0; step < repeats; step++) { + var isig = this.getRDSignalBandWidth(nm_levels) + + // get the direction of meanshift vector for a bin + // it compares near by bins to get the direction of the vector; bin band defines the range of comparison + + var grad = new Array(nm_levels.length).fill(0); + + for (let i = 0; i < nm_levels.length; i++) { + const start_bin = Math.max(0, i - 3 * bin_band); + const end_bin = Math.min(nm_levels.length - 1, i + 3 * bin_band + 1); + // let bin_length = end_bin - start_bin + + for (let j = start_bin ; j <= end_bin; j++) { + // for (let j = i - 3 * bin_band; j <= i + 3 * bin_band + 1; j++) { + // if (j < 0 || j >= nm_levels.length) continue; + // if (Math.abs(i - j) >= nm_levels.length) continue; + + var g_value = (j - i) * Math.exp((-0.5 * (j - i) ** 2) / bin_band ** 2) * Math.exp(-0.5 * (nm_levels[i] - nm_levels[j]) ** 2 * isig[i]); + + grad[i] += g_value + } + } + // console.log("grad: ", grad) + // get the border; if there is a change of gradient, it is a border + var border = new Array(); + for (var i = 0; i < grad.length - 1; i++) { + if ((grad[i] < 0) & (grad[i + 1] >= 0)) border.push(i); + } + + border.push(grad.length - 1) + border = border.concat(mask_borders).sort((a, b) => a - b) + border = Array.from(new Set(border)) + + var pb = 0; + for (var i = 0; i < border.length; i++) { + var range_array = nm_levels.slice(pb, border[i] + 1) + var range_mean = range_array.reduce((acc, n) => acc + n) / range_array.length + + nm_levels.fill(range_mean, pb, border[i] + 1) + pb = border[i] + 1 + } + } + + for (var i = 0, j = 0; i < levels.length; i++) { + if (not_masked[i]) { + levels[i] = nm_levels[j] + j++ + } + } + + //get the border + var border = new Array(); + for (var i = 0; i < levels.length - 1; i++) { + //if(i == levels.length -1) continue; + var diff = Math.abs(levels[i + 1] - levels[i]); + + if (diff > 0.01) border.push(i + 1); + } + + border.unshift(0); + border.push(levels.length); + + // reset the mask + masked = new Array(this.wigFeatures.length).fill(false); + + // check the borders + for (var i = 1; i < border.length; i++) { + var seg = [border[i - 1], border[i]] + var seg_left = [border[i - 1], border[i - 1]] + if (i > 1) { seg_left[0] = border[i - 2] } else continue; + + var seg_right = [border[i], border[i]]; + if (i < border.length - 1) { seg_right[1] = border[i + 1] } else continue; + + var n = seg[1] - seg[0]; + var n_left = seg_left[1] - seg_left[0]; + var n_right = seg_right[1] - seg_right[0]; + if (n <= 1) continue; + var seg_array = new DataStat(levels.slice(seg[0], seg[1])); + + if (n_right <= 15 || n_left <= 15 || n <= 15) { + var ns = 1.8 * Math.sqrt(levels[seg_left[0]] / this.globalMean) * this.globalStd; + if (Math.abs(levels[seg_left[0]] - levels[seg[0]]) < ns) { continue } + + ns = 1.8 * Math.sqrt(levels[seg_right[0]] / this.globalMean) * this.globalStd; + if (Math.abs(levels[seg_right[0]] - levels[seg[0]]) < ns) { continue } + } else { + var seg_left_array = levels.slice(seg_left[0], seg_left[1]) + var seg_left_1 = new DataStat(seg_left_array); + + var seg_right_array = levels.slice(seg_right[0], seg_right[1]) + var seg_right_1 = new DataStat(seg_right_array); + + var ttest_2sample_1 = t_test_2_samples(seg_array.mean, seg_array.std, seg_array.data.length, + seg_left_1.mean, seg_left_1.std, seg_left_1.data.length); + if (ttest_2sample_1 > (0.01 / genome_size) * this.binSize * (n + n_left)) { continue } + + var ttest_2sample_2 = t_test_2_samples(seg_array.mean, seg_array.std, seg_array.data.length, + seg_right_1.mean, seg_right_1.std, seg_right_1.data.length); + if (ttest_2sample_2 > (0.01 / genome_size) * this.binSize * (n + n_right)) { continue } + } + + var ttest_1sample_1 = t_test_1_sample(this.globalMean, seg_array.mean, seg_array.std, seg_array.data.length) + if (ttest_1sample_1 > 0.05) { continue } + let segments_rd = nm_levels.slice(seg[0], seg[1]) + // console.log("segments_rd: ", segments_rd) + var raw_seg_data = new DataStat(segments_rd); + + masked.fill(true, seg[0], seg[1]); + levels.fill(raw_seg_data.mean, seg[0], seg[1]); + } + + }); + // console.log("after applying partition: ", levels) + chrLevels[chr] = levels + // break + } + + return chrLevels + } + + cnvCalling(levels) { + + // console.log("levels: ", levels) + var delta = 0.25 * this.globalMean + var min = this.globalMean - delta, max = this.globalMean + delta; + var normal_genome_size = 2971000000 + + // var levels = this.meanShiftCaller(bin_size) + // var levels = this.MeanShiftCallerV2(bin_size) + + var merged_level = {} + var cnv_levels = [] + + Object.entries(levels).forEach(([chr, chr_levels]) => { + + var done = false + while (!done) { + done = true + + // + // get the borders + // + var borders = new Array(1).fill(0); + for (let i = 0; i < chr_levels.length - 1; i++) { + var diff = Math.abs(chr_levels[i + 1] - chr_levels[i]); + if (diff > 0.01) borders.push(i + 1); + } + borders.push(chr_levels.length); + + for (let ix = 0; ix < borders.length - 2; ix++) { + var v1 = Math.abs(chr_levels[borders[ix]] - chr_levels[borders[ix + 1]]) + // console.log(ix, v1); + if (v1 < delta) { + var v2 = v1 + 1, v3 = v1 + 1; + + if (ix > 0) { v2 = Math.abs(chr_levels[borders[ix]] - chr_levels[borders[ix - 1]]) } + if (ix < borders.length - 3) { v3 = Math.abs(levels[borders[ix + 1]] - chr_levels[borders[ix + 2]]) } + + if (v1 < v2 && v1 < v3) { + done = false + + var tmp_array = new DataStat(chr_levels.slice(borders[ix], borders[ix + 2])) + chr_levels.fill(tmp_array.mean, borders[ix], borders[ix + 2]); + borders.splice(ix + 1, ix + 1); + } + } + } + } + // console.log('updated levels', chr_levels) + + // var chr_rd = this.rd[chr] + var chr_rd = [] + Object.entries(this.wigFeatures[chr]).forEach(([bin, binDict]) => { chr_rd.push(binDict.binScore) }); + // console.log('cnv_calling', chr_rd) + + // + // Calling Segments + // + + //var flags = [""] * levels.length; + var flags = new Array(chr_levels.length).fill("") + var segments = [] + + // console.log('default levels', chr, chr_levels) + var b = 0 + var pval = (0.05 * this.binSize) / normal_genome_size + while (b < chr_levels.length) { + var b0 = b, border_start = b; + while ((b < chr_levels.length) & (chr_levels[b] < min)) { + b += 1 + } + var border_end = b + + if (border_end > border_start + 1) { + // console.log('border min', border_start, border_end) + var adj = adjustToEvalue(this.globalMean, this.std, chr_rd, border_start, border_end, pval) + // console.log(adj) + if (adj) { + var border_start, border_end = adj; + segments.push([border_start, border_end + 1]) + flags.fill("D", border_start, border_end) + } + } + border_start = b; + while ((b < chr_levels.length) & (chr_levels[b] > max)) { b += 1 } + border_end = b + + if (border_end > border_start + 1) { + adj = adjustToEvalue(this.globalMean, this.std, chr_rd, border_start, border_end, pval) + // console.log(adj) + if (adj) { + border_start, (border_end = adj); + segments.push([border_start, border_end, +1]); + // flags[bs:be] = ["A"] * (be - bs) + flags.fill("A", border_start, border_end); + } + } + if (b == b0) b += 1; + } + + //console.log(chr, 'segments', segments) + // + // Calling additional deletions + // + b = 0; + while (b < chr_levels.length) { + while ((b < chr_levels.length) & (flags[b] != "")) b += 1; + border_start = b; + while ((b < chr_levels.length) & (chr_levels[b] < min)) b += 1; + border_end = b; + if (border_end > border_start + 1) { + if (gaussianEValue(this.globalMean, this.std, chr_rd, border_start, border_end) < 0.05 / normal_genome_size) { + segments.push([border_start, border_end, -1]) + flags.fill(["d"] * (border_end - border_start), border_start, border_end) + } + b -= 1; + } + b += 1; + } + + b = 0; + var cf; + if (b < chr_levels.length) cf = flags[b]; + + border_start = 0; + + //var merge = [...this.rd] + var merge = [...chr_rd] + // console.log('initial merge', merge) + while (b < chr_levels.length) { + while (flags[b] == cf) { + b += 1; + if (b >= flags.length) break; + } + if (b > border_start) { + var merge_arr = new DataStat(merge.slice(border_start, b)); + merge.fill(merge_arr.mean, border_start, b); + } + if (b < chr_levels.length) cf = flags[b]; + border_start = b; + } + + merged_level[chr] = merge + + + // + // calculate calls + // + b = 0 + while (b < chr_levels.length) { + cf = flags[b] + if (cf == "") { + b += 1 + continue + } + border_start = b + while (b < chr_levels.length & cf == flags[b]) { b += 1 } + var border_arr = new DataStat(merge.slice(border_start, b)) + let cnv = border_arr.mean / this.globalMean + let event_type; + if (cf == "D") { + event_type = "deletion" + } else { + event_type = "duplication" + } + let cnv_dict = { + chr: chr, + start: this.binSize * border_start + 1, + end: this.binSize * b, + size: this.binSize * (b - border_start + 1), + value: cnv * 2, + event_type: event_type + } + cnv_levels.push(cnv_dict) + + } + }); + + return [merged_level, cnv_levels] + } + + +} + function erf(x) { var m = 1.0, s = 1.0, sum = x * 1.0; for (var i = 1; i < 50; i++) { @@ -658,4 +1086,4 @@ export class Partition { -export default { Partition }; \ No newline at end of file +export default { Partition, MeanShiftCaller }; \ No newline at end of file diff --git a/js/cnvpytor/baseCNVpytorVCF.js b/js/cnvpytor/baseCNVpytorVCF.js new file mode 100644 index 000000000..92708b76e --- /dev/null +++ b/js/cnvpytor/baseCNVpytorVCF.js @@ -0,0 +1,384 @@ + +import {BGZip, igvxhr, StringUtils} from "../../node_modules/igv-utils/src/index.js" +import lm from '../vendor/lm-esm.js' + + +class FitingMethod { + /** + * + * @param {*} binScores ; array containing list of rd values + */ + constructor(binScores) { + this.binScores = binScores + } + + get_histogram(){ + /** + * Get the 2D histogram of the array. + */ + + let max_bin = Math.max(...this.binScores); + let rd_bin_size = Math.floor(max_bin / 1000); + let rd_bins = this.range(0, Math.floor(max_bin / rd_bin_size) * rd_bin_size + rd_bin_size, rd_bin_size); + const { counts, bins } = this.histogram(this.binScores, rd_bins); + return { counts, bins }; + } + + range(start, end, step) { + const result = []; + for (let i = start; i < end; i += step) { + result.push(i); + } + return result; + } + + histogram(data, bins) { + const counts = Array(bins.length).fill(0); + + data.forEach(value => { + for (let i = 0; i < bins.length - 1; i++) { + if (value >= bins[i] && value < bins[i + 1]) { + counts[i]++; + break; + } + } + }); + + return { counts, bins }; + } + + normal_distribution([a, x0, sigma]) { + /** + * Normal distribution. + * + * @param {number} a - Area. + * @param {number} x0 - Mean value. + * @param {number} sigma - Sigma. + * @returns {number} - Value of distribution at x. + */ + // console.log('Model Parameters', a, x0, sigma) + return (x) => a * Math.exp(-Math.pow(x - x0, 2) / (2 * Math.pow(sigma, 2))) / (Math.sqrt(2 * Math.PI) * sigma); + } + + get_initial_model_values(x, y){ + /** + * calculate the initial values for the normal distribution + * x: The x values of the data. + * y: The y values of the data. + * return [area, mean, sigma]: The initial values for the normal distribution. + */ + + const sumY = y.slice(1, -1).reduce((acc, val) => acc + val, 0); + if (sumY === 0) { + console.debug("Problem with fit: all data points have zero value. Return zeros instead fit parameters!"); + return [0, 0, 0], null; + } + // calculate area + const bin_width = x[1] - x[0]; + const area = bin_width * sumY; + + // Calculate the mean + const weightedSum = x.reduce((acc, val, idx) => acc + val * y[idx], 0); + const mean = weightedSum / sumY; + + // Calculate the sigma + const weightedSquareSum = x.reduce((acc, val, idx) => acc + y[idx] * Math.pow(val - mean, 2), 0); + const sigma = Math.sqrt(weightedSquareSum / sumY); + // console.log('Initial Model Parameters', area, mean, sigma) + return [area, mean, sigma]; + } + + normal_fit(x, y) { + /** + * Fit a normal distribution to the histogram using levenberg-marquardt algorithm. + * x: The x values of the data. + * y: The y values of the data. + * return model: The model of the normal distribution. + */ + + // Initial values were calculated for arithmetic mean and standard deviation. + let initialValues = this.get_initial_model_values(x, y); + + // paramters for the levenberg-marquardt algorithm + const options = { + initialValues: initialValues, + maxIterations: 100, + } + + // Fit the model + const model = lm({x, y}, this.normal_distribution, options); + + return model; + } + +} + +class FetchGCInfo { + /** + * Creates an instance of gcCorrection. + * @param {number} binSize - binSize + * @param {string} refGenome - Reference genome name. + */ + constructor(binSize, refGenome){ + this.binSize = binSize + this.refGenome = refGenome + + // Determine the appropriate GC bin size + this.gcBin = this.getGCbinSize() + } + + /** + * Determines the appropriate GC bin size. + * @returns {number|boolean} The GC bin size or false if not found. + */ + getGCbinSize(){ + for (let gcBin of [100000, 10000]){ + // Check if the binSize is a multiple of the current gcBin + if (this.binSize % gcBin == 0) return gcBin; + } + // Return false if no appropriate GC bin size is found + return false + } + + /** + * Fetches GC values for a predefined gc bin sizes (i.e., 10k, 100k) from a remote JSON file. + * @returns {Object} The GC values. + */ + async getBinGC(){ + const gcData = {}; + + // If no appropriate GC bin size is found, return an empty object + if (!this.gcBin) { + return gcData; + } + + // URL for the GC info JSON file + const gcInfoURL = 'https://storage.googleapis.com/cnvpytor_data/gcInfoData/GCinfo.json' + + try { + // Load the GC info JSON file + const gcInfoJson = await igvxhr.loadJson(gcInfoURL, {timeout: 5000}) + + // Find the reference genome data within the JSON file + const gcRef = gcInfoJson.find(item => item.id === this.refGenome); + if (!gcRef){ + console.warn("GC data not found for ", this.refGenome); + return gcData; + } + + // Find the matching bin size data + const matchingInfo = gcRef.Bins.find(bin => bin.BinSize === this.gcBin); + if (!matchingInfo) { + console.warn("GC data not found for ", this.refGenome, " Bin : ", this.gcBin); + return gcData; + } + + // Construct the URL to fetch the GC data file + const parentDirectory = gcInfoURL.split('/').slice(0, -1).join('/'); + let gcURL = `${parentDirectory}${matchingInfo.fileURL}`; + + // Fetch the GC data file + const data = await igvxhr.load(gcURL, {}); + + // Parse the GC data file and populate the gcData object + data.split("\n").forEach((row) => { + if (row.trim() !== "") { + const [refName, start, gcContent, gcCount, atCount] = row.split("\t"); + if (!gcData[refName]) { + gcData[refName] = []; + } + gcData[refName].push({ + start: +start, + gcContent: +gcContent, + gcCount: +gcCount, + atCount: +atCount, + }); + } + }); + + } catch(e){ + console.error(e); + console.warn("Errors loading GC correction data."); + } + + return gcData; + + } + +} + + +class baseCNVpytorVCF extends FetchGCInfo { + /** + * Class for gcCorrection + * + * @param {*} wigFeatures + * @param {*} binSize + * @param {*} refGenome + */ + constructor(wigFeatures, binSize, refGenome){ + super(binSize, refGenome) + this.wigFeatures = wigFeatures + this.binSize = binSize + this.refGenome = refGenome + + } + async apply_gcCorrection(){ + // applying fitting method and defineing the variables + if (!this.wigFeatures) { + console.error("BinScore data is not available."); + return null; + } + + // Extract all binScore values into a single array + const binScores = Object.values(this.wigFeatures).reduce( + (binResult, bin) => { return binResult.concat(bin.filter(a => a.binScore > 0).map(a => a.binScore)) }, [] + ) + + let data_fitter = new FitingMethod(binScores) + + const { counts: dist_p, bins } = data_fitter.get_histogram() + let model_parameters = data_fitter.normal_fit(bins, dist_p) + // return { globalMean: model_parameters.parameterValues[1], globalStd: model_parameters.parameterValues[2]} + this.globalMean = model_parameters.parameterValues[1] + this.globalStd = model_parameters.parameterValues[2] + // console.log("globalmean", this.globalMean) + + + this.gcData = await this.getBinGC() + this.gcFlag = Object.keys(this.gcData).length > 0 ? true : false; // Flag indicating whether GC data is available + this.binScoreField = this.gcFlag ? "gcCorrectedBinScore": "binScore" ; + + this.getGcCorrectionSignal(this.globalMean) + + } + + /** + * Applies GC correction to the bin scores. + * @param {number} rdGlobalMean - The global mean read depth. + */ + getGcCorrectionSignal(rdGlobalMean){ + let gcRDMean = this.getGcCorrection(rdGlobalMean) + Object.keys(this.wigFeatures).forEach(chr => { + this.wigFeatures[chr].forEach(bin => { + if (bin.binScore){ + bin.gcCorrectedBinScore = Math.round(gcRDMean[bin.gc] * bin.binScore) + }else{ + bin.gcCorrectedBinScore = null + } + } + ); + }); + } + + /** + * Computes the GC correction values based on the global mean read depth. + * @param {number} rdGlobalMean - The global mean for read depth. + * @returns {Object} An object containing the GC correction values indexed by GC percentage. + */ + getGcCorrection(rdGlobalMean){ + + const gcRDMean = {} + if(this.gcFlag){ + let gcBin = this.getGCbinSize() + const gcRD = {}; + + let gcBinFactor = parseInt(this.binSize/gcBin) + for (let chr in this.wigFeatures){ + for (let k=0; k< this.wigFeatures[chr].length; k++){ + + // collect GC related values for a bin + let baseInfo = {'AT': 0, 'GC': 0} + for (let j = k * gcBinFactor; j < k * gcBinFactor + gcBinFactor; j++){ + if (this.gcData[chr][j]){ + baseInfo['GC'] += this.gcData[chr][j].gcCount + baseInfo['AT'] += this.gcData[chr][j].atCount + } + } + + let gcValue = Math.round((baseInfo['GC'] * 100)/(baseInfo['GC'] + baseInfo['AT'])) + + this.wigFeatures[chr][k].gc = gcValue + if (!gcRD[gcValue]) { + gcRD[gcValue] = []; + } + if (this.wigFeatures[chr][k].binScore){ + gcRD[gcValue].push(this.wigFeatures[chr][k].binScore) + } + } + } + + // get average rd for bins that has same gc content + Object.keys(gcRD).forEach(gc_value => { + + // apply fitting method to get the gc average + let gcRDBins = gcRD[gc_value] + let gcMean + if(gcRDBins.length < 4){ + gcMean = gcRDBins.reduce((sum, val) => sum + val, 0) / gcRDBins.length; + }else{ + let data_fitter = new FitingMethod(gcRD[gc_value]) + + const { counts: dist_p, bins } = data_fitter.get_histogram() + let model_parameters = data_fitter.normal_fit(bins, dist_p) + gcMean = model_parameters.parameterValues[1] + } + gcRDMean[gc_value] = rdGlobalMean/gcMean + + }); + + } + return gcRDMean + } + + /** + * Formats the data structure for output. + * @param {string} feature_column - The column name of the feature to be scaled. + * @param {number} [scaling_factor=1] - The factor by which to scale the feature values. + * @returns {Array} The formatted data structure. + */ + formatDataStructure(feature_column, scaling_factor = 1) { + const results = [] + for (const [chr, wig] of Object.entries(this.wigFeatures)) { + + for(let sample of wig){ + var new_sample = { ...sample } + if (scaling_factor != 1) { + new_sample.value = sample[feature_column] / scaling_factor * 2 + } + results.push(new_sample) + } + } + + return results + } + + formatDataStructure_BAF(feature_column, scaling_factor = -1) { + const baf1 = [] + const baf2 = [] + for (const [chr, wig] of Object.entries(this.wigFeatures)) { + + wig.forEach(sample => { + + var baf1_value = { ...sample } + var baf2_value = { ...sample } + + let value = sample[feature_column] + if (value != 0.5){ + baf2_value.value = scaling_factor * (1 - value) + baf2.push(baf2_value) + } + baf1_value.value = scaling_factor * value + baf1.push(baf1_value) + + }) + } + + return [baf1, baf2] + } + + + +} + +export default baseCNVpytorVCF; diff --git a/js/cnvpytor/cnvpytorTrack.js b/js/cnvpytor/cnvpytorTrack.js index 7cbcc6738..3f282c421 100644 --- a/js/cnvpytor/cnvpytorTrack.js +++ b/js/cnvpytor/cnvpytorTrack.js @@ -33,6 +33,7 @@ class CNVPytorTrack extends TrackBase { this.signal_name = config.signal_name || "rd_snp" this.cnv_caller = config.cnv_caller || '2D' this.colors = config.colors || ['gray', 'black', 'green', 'blue'] + this.hasChroms = {} super.init(config) } @@ -85,7 +86,10 @@ class CNVPytorTrack extends TrackBase { allVariants = this.featureSource.reader.features } - const cnvpytor_obj = new CNVpytorVCF(allVariants, this.bin_size) + const refGenome = this.browser.config.genome + + // Initializing CNVpytorVCF class + const cnvpytor_obj = new CNVpytorVCF(allVariants, this.bin_size, refGenome) let wigFeatures let bafFeatures this.wigFeatures_obj = {} @@ -121,8 +125,21 @@ class CNVPytorTrack extends TrackBase { } else { this.cnvpytor_obj = new HDF5IndexedReader(this.config.url, this.bin_size) - this.wigFeatures_obj = await this.cnvpytor_obj.get_rd_signal(this.bin_size) - this.available_bins = this.cnvpytor_obj.available_bins + // get chrom list that currently user viewing + let chroms = [ ...new Set(this.browser.referenceFrameList.map(val => val.chr))] + + let aliasRecord = this.getAliasChromsList(chroms) + this.wigFeatures_obj = await this.cnvpytor_obj.get_rd_signal(this.bin_size, aliasRecord) + + // Save the processed chroms names to check later for the availability + this.update_hasChroms(this.wigFeatures_obj, chroms) + + this.available_bins = this.cnvpytor_obj.availableBins + // reset the bin size if its not exits + if(! this.available_bins.includes(this.bin_size)){ + this.bin_size = this.available_bins.at(-1); + } + this.available_callers = this.cnvpytor_obj.callers this.set_available_callers() } @@ -135,18 +152,17 @@ class CNVPytorTrack extends TrackBase { for (let bin_size in this.wigFeatures_obj) { let i = 0 - for (const [signal_name, wig] of Object.entries(this.wigFeatures_obj[bin_size])) { if (this.signals.includes(signal_name)) { let tconf = {} tconf.type = "wig" + tconf.isMergedTrack = true tconf.features = wig tconf.name = signal_name tconf.color = this.signal_colors.filter(x => x.singal_name === signal_name).map(x => x.color) const t = await this.browser.createTrack(tconf) if (t) { - t.isMergedTrack = true t.autoscale = false // Scaling done from merged track this.tracks.push(t) } else { @@ -180,6 +196,15 @@ class CNVPytorTrack extends TrackBase { return Promise.all(p) } + getAliasChromsList(chroms){ + let aliasRecord = chroms.map(chr => { + let records = this.browser.genome.chromAlias.aliasRecordCache.get(chr) + return Object.values(records) + }) + aliasRecord = aliasRecord.flat() + return aliasRecord + } + set_available_callers() { if (!this.available_callers.includes(this.cnv_caller)) { if (this.available_callers.length > 0) { @@ -316,7 +341,18 @@ class CNVPytorTrack extends TrackBase { const p = [] if (!(bin_size in this.wigFeatures_obj)) { - this.wigFeatures_obj = {...this.wigFeatures_obj, ...await this.cnvpytor_obj.get_rd_signal(bin_size)} + if(Object.keys(this.hasChroms).length > 0) { + let chroms = [ ...new Set(this.browser.referenceFrameList.map(val => val.chr))] + if(chroms[0] == "all"){ + chroms = this.browser.genome.chromosomeNames + } + + this.wigFeatures_obj = {...this.wigFeatures_obj, ...await this.cnvpytor_obj.get_rd_signal(bin_size, chroms)} + this.update_hasChroms(this.wigFeatures_obj, chroms) + } else{ + this.wigFeatures_obj = {...this.wigFeatures_obj, ...await this.cnvpytor_obj.get_rd_signal(bin_size)} + } + } this.signals = this.get_signals() @@ -328,12 +364,12 @@ class CNVPytorTrack extends TrackBase { if (this.signals.includes(signal_name)) { let tconf = {} tconf.type = "wig" + tconf.isMergedTrack = true tconf.features = wig tconf.name = signal_name tconf.color = this.signal_colors.filter(x => x.singal_name === signal_name).map(x => x.color) const t = await this.browser.createTrack(tconf) if (t) { - t.isMergedTrack = true t.autoscale = false // Scaling done from merged track this.tracks.push(t) } else { @@ -365,7 +401,54 @@ class CNVPytorTrack extends TrackBase { return Promise.all(p) } + update_hasChroms(wigFeatures, chroms){ + for (let binSize in wigFeatures){ + if (!this.hasChroms[binSize]) { + this.hasChroms[binSize] = new Set(); + } + chroms.forEach(item => this.hasChroms[binSize].add(item)) + + } + return this.hasChroms + + } + async getFeatures(chr, bpStart, bpEnd, bpPerPixel) { + + if(Object.keys(this.hasChroms).length > 0) { + + // Need to find the current binSize + if (this.hasChroms[this.bin_size].size != 0){ + let chroms = [ ...new Set(this.browser.referenceFrameList.map(val => val.chr))] + if(chroms[0] == "all"){ + chroms = this.browser.genome.chromosomeNames + } + let newChroms = chroms.filter(val => !this.hasChroms[this.bin_size].has(val)) + + if(newChroms.length != 0){ + + let aliasRecords = this.getAliasChromsList(newChroms) + // update the hasChroms list + let tmp_wig = await this.cnvpytor_obj.get_rd_signal(this.bin_size, aliasRecords) + + this.update_hasChroms(tmp_wig, newChroms) + + // here we need to update the wig + // this part is probaby not required; code improve required + + for (let bin_size in tmp_wig){ + for (const [signal_name, wig] of Object.entries(tmp_wig[bin_size])) { + await this.wigFeatures_obj[bin_size][signal_name].push(...wig) + } + } + + for (let wig of this.tracks){ + await wig.featureSource.updateFeatures(this.wigFeatures_obj[this.bin_size][wig.name] ) + } + } + } + + } if (this.tracks) { const promises = this.tracks.map((t) => t.getFeatures(chr, bpStart, bpEnd, bpPerPixel)) diff --git a/js/cnvpytor/cnvpytorVCF.js b/js/cnvpytor/cnvpytorVCF.js index 4eb24a633..6ec2d9ced 100644 --- a/js/cnvpytor/cnvpytorVCF.js +++ b/js/cnvpytor/cnvpytorVCF.js @@ -14,13 +14,14 @@ class CNVpytorVCF { * * @param {Array} allVariants - An array containing all variants * @param {number} binSize - The bin size for processing variants. + * @param {string} refGenome - Reference genome name */ - constructor(allVariants, binSize) { + constructor(allVariants, binSize, refGenome) { this.allVariants = allVariants this.rowBinSize = 10000 this.binSize = binSize this.binFactor = parseInt(binSize / this.rowBinSize) - + this.refGenome = refGenome } /** @@ -81,45 +82,43 @@ class CNVpytorVCF { // get the chromosome names this.chromosomes = Object.keys(wigFeatures) + let delete_likelihood_scores = caller == 'ReadDepth' ? true : false; + // Step2: Update the binsize according to user provided bin size - var avgbin = this.adjust_bin_size(wigFeatures) + var avgbin = this.adjust_bin_size(wigFeatures, delete_likelihood_scores=delete_likelihood_scores) + // this is to save objects + // console.log("avgbin: ", avgbin) + // Step3: Run the CNV caller var finalFeatureSet if(caller == 'ReadDepth'){ - let caller_obj = new MeanShiftCaller(avgbin, this.binSize) - finalFeatureSet = caller_obj.ReadDepth_caller() + // ------------ new code + // console.log("setting up meanShift CNV calling") + let callerObj = new read_depth_caller.MeanShiftCaller(avgbin, this.binSize, this.refGenome) + + let processedBins = await callerObj.callMeanshift() + + // finalFeatureSet = [processedBins.binScore, processedBins.gcCorrectedBinScore, processedBins.segmentsCNV] + finalFeatureSet = [processedBins.binScore, processedBins.gcCorrectedBinScore, processedBins.segmentsCNV] + // var baf = this.formatDataStructure_BAF(avgbin, 'max_likelihood') + // var baf = callerObj.format_BAF_likelihood(avgbin) + var baf = callerObj.formatDataStructure_BAF('max_likelihood', -1) + - // finalFeatureSet = this.readDepthMeanshift(avgbin) - var baf = this.formatDataStructure_BAF(avgbin, 'max_likelihood') }else if(caller=='2D'){ - let caller_obj = new combined_caller.CombinedCaller(avgbin, this.binSize) + + let caller_obj = new combined_caller.CombinedCaller(avgbin, this.binSize, this.refGenome) let processed_bins = await caller_obj.call_2d() - finalFeatureSet = [processed_bins.binScore, [], processed_bins.segment_score] + finalFeatureSet = [processed_bins.binScore, processed_bins.gcCorrectedBinScore, processed_bins.segmentScore] var baf = caller_obj.formatDataStructure_BAF('max_likelihood', -1) + } return [finalFeatureSet, baf] } - - formatDataStructure(wigFeatures, feature_column, scaling_factor = 1) { - const results = [] - for (const [chr, wig] of Object.entries(wigFeatures)) { - - for(let sample of wig){ - var new_sample = { ...sample } - if (scaling_factor != 1) { - new_sample.value = sample[feature_column] / scaling_factor * 2 - } - results.push(new_sample) - } - } - - return results - } - format_BAF_likelihood(wigFeatures) { const results = [] @@ -149,51 +148,24 @@ class CNVpytorVCF { return sample } - + /* async getAllbins() { const bins = await this.computeDepthFeatures() const fitter = new g_utils.GetFit(bins) const distParams = fitter.fit_data() - return bins - } - - formatDataStructure_BAF(wigFeatures, feature_column, scaling_factor=-1) { - /* Rescale the BAF level from 0:1 to scaling_factpr:0*/ - const baf1 = [] - const baf2 = [] - for (const [chr, wig] of Object.entries(wigFeatures)) { - - for(let sample of wig) { - //delete sample.likelihood_score; - var baf1_value = { ...sample } - var baf2_value = { ...sample } - - let value = sample[feature_column] - if (value != 0.5){ - baf2_value.value = scaling_factor * (1 - value) - baf2.push(baf2_value) - } - baf1_value.value = scaling_factor * value - baf1.push(baf1_value) - - } - } - - - return [baf1, baf2] - } - + }*/ /** * Adjust the bin values to actual bin size * @param {*} wigFeatures - wig features after processing the varaints * @returns */ - adjust_bin_size(wigFeatures){ + adjust_bin_size(wigFeatures, delete_likelihood_scores=false){ var avgbin = {} + const scale = this.binSize/150 for (let chr of this.chromosomes) { if (!avgbin[chr]) { avgbin[chr] = [] } for (let k = 0; k < wigFeatures[chr].length / this.binFactor; k++) { @@ -205,7 +177,7 @@ class CNVpytorVCF { end: (featureBin + 1) * this.binSize, dp_count: 0, hets_count: 0, - binScore: 0, + binScore: null, likelihood_score: [], dp_sum_score: 0 } @@ -246,9 +218,18 @@ class CNVpytorVCF { } } } - avgbin[chr][k].binScore = parseInt(avgbin[chr][k].dp_sum_score / avgbin[chr][k].dp_count) * 100; + if (avgbin[chr][k].dp_count > 0) { + // avgbin[chr][k].binScore = parseInt(avgbin[chr][k].dp_sum_score / avgbin[chr][k].dp_count) * scale; + avgbin[chr][k].binScore = avgbin[chr][k].dp_sum_score / avgbin[chr][k].dp_count * scale; + } const updated_bin = this.get_max_min_score(avgbin[chr][k]) avgbin[chr][k].max_likelihood = updated_bin.value + + // delete the likelihood_score array; otherwise it takes too much memory + // this is used in the combined caller; so it can be deleted depending on the caller i.e., ReadDepth caller + if (delete_likelihood_scores) { + delete avgbin[chr][k].likelihood_score + } } } return avgbin @@ -260,55 +241,5 @@ function beta(a, b, p, phased = true) { return Math.pow(p, a) * Math.pow(1-p, b) + Math.pow(p, b) * Math.pow(1-p, a) } -class MeanShiftCaller{ - constructor(wigFeatures, binSize) { - this.wigFeatures = wigFeatures - this.binSize = binSize - } - - ReadDepth_caller(){ - // Get global mean and standrad deviation - var fit_info = new g_utils.GetFit(this.wigFeatures) - var [globamMean, globalStd] = fit_info.fit_data(); - - // Apply partition method - var partition = new read_depth_caller.Partition(this.wigFeatures, globamMean, globalStd); - var partition_array = partition.meanShiftCaller(this.binSize) - var caller_array = partition.cnv_calling() - - // Assign the partition values to each bin - Object.entries(this.wigFeatures).forEach(([chr, chr_rd]) => { - chr_rd.forEach((bin, index) => { - bin.partition_level = parseInt(partition_array[chr][index]) - bin.partition_call = parseInt(caller_array[0][chr][index]) - }); - }) - - var rawbinScore = this.formatDataStructure('binScore', globamMean) - var partition_levels = this.formatDataStructure('partition_level', globamMean) - var partition_call = this.formatDataStructure('partition_call', globamMean) - - return [rawbinScore, partition_levels, partition_call, caller_array[1]] - } - - formatDataStructure(feature_column, scaling_factor = 1) { - const results = [] - for (const [chr, wig] of Object.entries(this.wigFeatures)) { - - for(let sample of wig){ - var new_sample = { ...sample } - if (scaling_factor != 1) { - new_sample.value = sample[feature_column] / scaling_factor * 2 - } - results.push(new_sample) - } - } - - return results - } - - -} - export { CNVpytorVCF } diff --git a/js/vendor/lm-esm.js b/js/vendor/lm-esm.js new file mode 100644 index 000000000..466d921bc --- /dev/null +++ b/js/vendor/lm-esm.js @@ -0,0 +1,6238 @@ +// source: https://github.com/mljs/levenberg-marquardt/tree/rollup-web + +const toString$4 = Object.prototype.toString; + +function isAnyArray$5(object) { + return toString$4.call(object).endsWith('Array]'); +} + +function checkOptions(data, parameterizedFunction, options) { + let { + timeout, + minValues, + maxValues, + initialValues, + weights = 1, + damping = 1e-2, + dampingStepUp = 11, + dampingStepDown = 9, + maxIterations = 100, + errorTolerance = 1e-7, + centralDifference = false, + gradientDifference = 10e-2, + improvementThreshold = 1e-3, + } = options; + + if (damping <= 0) { + throw new Error('The damping option must be a positive number'); + } else if (!data.x || !data.y) { + throw new Error('The data parameter must have x and y elements'); + } else if ( + !isAnyArray$5(data.x) || + data.x.length < 2 || + !isAnyArray$5(data.y) || + data.y.length < 2 + ) { + throw new Error( + 'The data parameter elements must be an array with more than 2 points', + ); + } else if (data.x.length !== data.y.length) { + throw new Error('The data parameter elements must have the same size'); + } + + let parameters = + initialValues || new Array(parameterizedFunction.length).fill(1); + + let nbPoints = data.y.length; + let parLen = parameters.length; + maxValues = maxValues || new Array(parLen).fill(Number.MAX_SAFE_INTEGER); + minValues = minValues || new Array(parLen).fill(Number.MIN_SAFE_INTEGER); + + if (maxValues.length !== minValues.length) { + throw new Error('minValues and maxValues must be the same size'); + } + + if (!isAnyArray$5(parameters)) { + throw new Error('initialValues must be an array'); + } + + if (typeof gradientDifference === 'number') { + gradientDifference = new Array(parameters.length).fill(gradientDifference); + } else if (isAnyArray$5(gradientDifference)) { + if (gradientDifference.length !== parLen) { + gradientDifference = new Array(parLen).fill(gradientDifference[0]); + } + } else { + throw new Error( + 'gradientDifference should be a number or array with length equal to the number of parameters', + ); + } + + let filler; + if (typeof weights === 'number') { + let value = 1 / weights ** 2; + filler = () => value; + } else if (isAnyArray$5(weights)) { + if (weights.length < data.x.length) { + let value = 1 / weights[0] ** 2; + filler = () => value; + } else { + filler = (i) => 1 / weights[i] ** 2; + } + } else { + throw new Error( + 'weights should be a number or array with length equal to the number of data points', + ); + } + + let checkTimeout; + if (timeout !== undefined) { + if (typeof timeout !== 'number') { + throw new Error('timeout should be a number'); + } + let endTime = Date.now() + timeout * 1000; + checkTimeout = () => Date.now() > endTime; + } else { + checkTimeout = () => false; + } + + let weightSquare = new Array(data.x.length); + for (let i = 0; i < nbPoints; i++) { + weightSquare[i] = filler(i); + } + + return { + checkTimeout, + minValues, + maxValues, + parameters, + weightSquare, + damping, + dampingStepUp, + dampingStepDown, + maxIterations, + errorTolerance, + centralDifference, + gradientDifference, + improvementThreshold, + }; +} + +/** + * the sum of the weighted squares of the errors (or weighted residuals) between the data.y + * and the curve-fit function. + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} parameters - Array of current parameter values + * @param {function} parameterizedFunction - The parameters and returns a function with the independent variable as a parameter + * @param {Array} weightSquare - Square of weights + * @return {number} + */ +function errorCalculation( + data, + parameters, + parameterizedFunction, + weightSquare, +) { + let error = 0; + const func = parameterizedFunction(parameters); + for (let i = 0; i < data.x.length; i++) { + error += Math.pow(data.y[i] - func(data.x[i]), 2) / weightSquare[i]; + } + + return error; +} + +function getAugmentedNamespace(n) { + if (n.__esModule) return n; + var a = Object.defineProperty({}, '__esModule', {value: true}); + Object.keys(n).forEach(function (k) { + var d = Object.getOwnPropertyDescriptor(n, k); + Object.defineProperty(a, k, d.get ? d : { + enumerable: true, + get: function () { + return n[k]; + } + }); + }); + return a; +} + +var matrix = {}; + +// eslint-disable-next-line @typescript-eslint/unbound-method +const toString$3 = Object.prototype.toString; +/** + * Checks if an object is an instance of an Array (array or typed array, except those that contain bigint values). + * + * @param value - Object to check. + * @returns True if the object is an array or a typed array. + */ +function isAnyArray$4(value) { + const tag = toString$3.call(value); + return tag.endsWith('Array]') && !tag.includes('Big'); +} + +var libEsm = /*#__PURE__*/Object.freeze({ + __proto__: null, + isAnyArray: isAnyArray$4 +}); + +var require$$0 = /*@__PURE__*/getAugmentedNamespace(libEsm); + +// eslint-disable-next-line @typescript-eslint/unbound-method +const toString$2 = Object.prototype.toString; +/** + * Checks if an object is an instance of an Array (array or typed array, except those that contain bigint values). + * + * @param value - Object to check. + * @returns True if the object is an array or a typed array. + */ +function isAnyArray$3(value) { + const tag = toString$2.call(value); + return tag.endsWith('Array]') && !tag.includes('Big'); +} + +// eslint-disable-next-line @typescript-eslint/unbound-method +const toString$1 = Object.prototype.toString; +/** + * Checks if an object is an instance of an Array (array or typed array, except those that contain bigint values). + * + * @param value - Object to check. + * @returns True if the object is an array or a typed array. + */ +function isAnyArray$2(value) { + const tag = toString$1.call(value); + return tag.endsWith('Array]') && !tag.includes('Big'); +} + +function max(input) { + var options = arguments.length > 1 && arguments[1] !== undefined ? arguments[1] : {}; + + if (!isAnyArray$2(input)) { + throw new TypeError('input must be an array'); + } + + if (input.length === 0) { + throw new TypeError('input must not be empty'); + } + + var _options$fromIndex = options.fromIndex, + fromIndex = _options$fromIndex === void 0 ? 0 : _options$fromIndex, + _options$toIndex = options.toIndex, + toIndex = _options$toIndex === void 0 ? input.length : _options$toIndex; + + if (fromIndex < 0 || fromIndex >= input.length || !Number.isInteger(fromIndex)) { + throw new Error('fromIndex must be a positive integer smaller than length'); + } + + if (toIndex <= fromIndex || toIndex > input.length || !Number.isInteger(toIndex)) { + throw new Error('toIndex must be an integer greater than fromIndex and at most equal to length'); + } + + var maxValue = input[fromIndex]; + + for (var i = fromIndex + 1; i < toIndex; i++) { + if (input[i] > maxValue) maxValue = input[i]; + } + + return maxValue; +} + +// eslint-disable-next-line @typescript-eslint/unbound-method +const toString = Object.prototype.toString; +/** + * Checks if an object is an instance of an Array (array or typed array, except those that contain bigint values). + * + * @param value - Object to check. + * @returns True if the object is an array or a typed array. + */ +function isAnyArray$1(value) { + const tag = toString.call(value); + return tag.endsWith('Array]') && !tag.includes('Big'); +} + +function min(input) { + var options = arguments.length > 1 && arguments[1] !== undefined ? arguments[1] : {}; + + if (!isAnyArray$1(input)) { + throw new TypeError('input must be an array'); + } + + if (input.length === 0) { + throw new TypeError('input must not be empty'); + } + + var _options$fromIndex = options.fromIndex, + fromIndex = _options$fromIndex === void 0 ? 0 : _options$fromIndex, + _options$toIndex = options.toIndex, + toIndex = _options$toIndex === void 0 ? input.length : _options$toIndex; + + if (fromIndex < 0 || fromIndex >= input.length || !Number.isInteger(fromIndex)) { + throw new Error('fromIndex must be a positive integer smaller than length'); + } + + if (toIndex <= fromIndex || toIndex > input.length || !Number.isInteger(toIndex)) { + throw new Error('toIndex must be an integer greater than fromIndex and at most equal to length'); + } + + var minValue = input[fromIndex]; + + for (var i = fromIndex + 1; i < toIndex; i++) { + if (input[i] < minValue) minValue = input[i]; + } + + return minValue; +} + +function rescale$1(input) { + var options = arguments.length > 1 && arguments[1] !== undefined ? arguments[1] : {}; + + if (!isAnyArray$3(input)) { + throw new TypeError('input must be an array'); + } else if (input.length === 0) { + throw new TypeError('input must not be empty'); + } + + var output; + + if (options.output !== undefined) { + if (!isAnyArray$3(options.output)) { + throw new TypeError('output option must be an array if specified'); + } + + output = options.output; + } else { + output = new Array(input.length); + } + + var currentMin = min(input); + var currentMax = max(input); + + if (currentMin === currentMax) { + throw new RangeError('minimum and maximum input values are equal. Cannot rescale a constant array'); + } + + var _options$min = options.min, + minValue = _options$min === void 0 ? options.autoMinMax ? currentMin : 0 : _options$min, + _options$max = options.max, + maxValue = _options$max === void 0 ? options.autoMinMax ? currentMax : 1 : _options$max; + + if (minValue >= maxValue) { + throw new RangeError('min option must be smaller than max option'); + } + + var factor = (maxValue - minValue) / (currentMax - currentMin); + + for (var i = 0; i < input.length; i++) { + output[i] = (input[i] - currentMin) * factor + minValue; + } + + return output; +} + +var libEs6 = /*#__PURE__*/Object.freeze({ + __proto__: null, + 'default': rescale$1 +}); + +var require$$1 = /*@__PURE__*/getAugmentedNamespace(libEs6); + +Object.defineProperty(matrix, '__esModule', { value: true }); + +var isAnyArray = require$$0; +var rescale = require$$1; + +const indent = ' '.repeat(2); +const indentData = ' '.repeat(4); + +/** + * @this {Matrix} + * @returns {string} + */ +function inspectMatrix() { + return inspectMatrixWithOptions(this); +} + +function inspectMatrixWithOptions(matrix, options = {}) { + const { + maxRows = 15, + maxColumns = 10, + maxNumSize = 8, + padMinus = 'auto', + } = options; + return `${matrix.constructor.name} { +${indent}[ +${indentData}${inspectData(matrix, maxRows, maxColumns, maxNumSize, padMinus)} +${indent}] +${indent}rows: ${matrix.rows} +${indent}columns: ${matrix.columns} +}`; +} + +function inspectData(matrix, maxRows, maxColumns, maxNumSize, padMinus) { + const { rows, columns } = matrix; + const maxI = Math.min(rows, maxRows); + const maxJ = Math.min(columns, maxColumns); + const result = []; + + if (padMinus === 'auto') { + padMinus = false; + loop: for (let i = 0; i < maxI; i++) { + for (let j = 0; j < maxJ; j++) { + if (matrix.get(i, j) < 0) { + padMinus = true; + break loop; + } + } + } + } + + for (let i = 0; i < maxI; i++) { + let line = []; + for (let j = 0; j < maxJ; j++) { + line.push(formatNumber(matrix.get(i, j), maxNumSize, padMinus)); + } + result.push(`${line.join(' ')}`); + } + if (maxJ !== columns) { + result[result.length - 1] += ` ... ${columns - maxColumns} more columns`; + } + if (maxI !== rows) { + result.push(`... ${rows - maxRows} more rows`); + } + return result.join(`\n${indentData}`); +} + +function formatNumber(num, maxNumSize, padMinus) { + return ( + num >= 0 && padMinus + ? ` ${formatNumber2(num, maxNumSize - 1)}` + : formatNumber2(num, maxNumSize) + ).padEnd(maxNumSize); +} + +function formatNumber2(num, len) { + // small.length numbers should be as is + let str = num.toString(); + if (str.length <= len) return str; + + // (7)'0.00123' is better then (7)'1.23e-2' + // (8)'0.000123' is worse then (7)'1.23e-3', + let fix = num.toFixed(len); + if (fix.length > len) { + fix = num.toFixed(Math.max(0, len - (fix.length - len))); + } + if ( + fix.length <= len && + !fix.startsWith('0.000') && + !fix.startsWith('-0.000') + ) { + return fix; + } + + // well, if it's still too long the user should've used longer numbers + let exp = num.toExponential(len); + if (exp.length > len) { + exp = num.toExponential(Math.max(0, len - (exp.length - len))); + } + return exp.slice(0); +} + +function installMathOperations(AbstractMatrix, Matrix) { + AbstractMatrix.prototype.add = function add(value) { + if (typeof value === 'number') return this.addS(value); + return this.addM(value); + }; + + AbstractMatrix.prototype.addS = function addS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) + value); + } + } + return this; + }; + + AbstractMatrix.prototype.addM = function addM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) + matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.add = function add(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.add(value); + }; + + AbstractMatrix.prototype.sub = function sub(value) { + if (typeof value === 'number') return this.subS(value); + return this.subM(value); + }; + + AbstractMatrix.prototype.subS = function subS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) - value); + } + } + return this; + }; + + AbstractMatrix.prototype.subM = function subM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) - matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.sub = function sub(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.sub(value); + }; + AbstractMatrix.prototype.subtract = AbstractMatrix.prototype.sub; + AbstractMatrix.prototype.subtractS = AbstractMatrix.prototype.subS; + AbstractMatrix.prototype.subtractM = AbstractMatrix.prototype.subM; + AbstractMatrix.subtract = AbstractMatrix.sub; + + AbstractMatrix.prototype.mul = function mul(value) { + if (typeof value === 'number') return this.mulS(value); + return this.mulM(value); + }; + + AbstractMatrix.prototype.mulS = function mulS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) * value); + } + } + return this; + }; + + AbstractMatrix.prototype.mulM = function mulM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) * matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.mul = function mul(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.mul(value); + }; + AbstractMatrix.prototype.multiply = AbstractMatrix.prototype.mul; + AbstractMatrix.prototype.multiplyS = AbstractMatrix.prototype.mulS; + AbstractMatrix.prototype.multiplyM = AbstractMatrix.prototype.mulM; + AbstractMatrix.multiply = AbstractMatrix.mul; + + AbstractMatrix.prototype.div = function div(value) { + if (typeof value === 'number') return this.divS(value); + return this.divM(value); + }; + + AbstractMatrix.prototype.divS = function divS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) / value); + } + } + return this; + }; + + AbstractMatrix.prototype.divM = function divM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) / matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.div = function div(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.div(value); + }; + AbstractMatrix.prototype.divide = AbstractMatrix.prototype.div; + AbstractMatrix.prototype.divideS = AbstractMatrix.prototype.divS; + AbstractMatrix.prototype.divideM = AbstractMatrix.prototype.divM; + AbstractMatrix.divide = AbstractMatrix.div; + + AbstractMatrix.prototype.mod = function mod(value) { + if (typeof value === 'number') return this.modS(value); + return this.modM(value); + }; + + AbstractMatrix.prototype.modS = function modS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) % value); + } + } + return this; + }; + + AbstractMatrix.prototype.modM = function modM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) % matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.mod = function mod(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.mod(value); + }; + AbstractMatrix.prototype.modulus = AbstractMatrix.prototype.mod; + AbstractMatrix.prototype.modulusS = AbstractMatrix.prototype.modS; + AbstractMatrix.prototype.modulusM = AbstractMatrix.prototype.modM; + AbstractMatrix.modulus = AbstractMatrix.mod; + + AbstractMatrix.prototype.and = function and(value) { + if (typeof value === 'number') return this.andS(value); + return this.andM(value); + }; + + AbstractMatrix.prototype.andS = function andS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) & value); + } + } + return this; + }; + + AbstractMatrix.prototype.andM = function andM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) & matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.and = function and(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.and(value); + }; + + AbstractMatrix.prototype.or = function or(value) { + if (typeof value === 'number') return this.orS(value); + return this.orM(value); + }; + + AbstractMatrix.prototype.orS = function orS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) | value); + } + } + return this; + }; + + AbstractMatrix.prototype.orM = function orM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) | matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.or = function or(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.or(value); + }; + + AbstractMatrix.prototype.xor = function xor(value) { + if (typeof value === 'number') return this.xorS(value); + return this.xorM(value); + }; + + AbstractMatrix.prototype.xorS = function xorS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) ^ value); + } + } + return this; + }; + + AbstractMatrix.prototype.xorM = function xorM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) ^ matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.xor = function xor(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.xor(value); + }; + + AbstractMatrix.prototype.leftShift = function leftShift(value) { + if (typeof value === 'number') return this.leftShiftS(value); + return this.leftShiftM(value); + }; + + AbstractMatrix.prototype.leftShiftS = function leftShiftS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) << value); + } + } + return this; + }; + + AbstractMatrix.prototype.leftShiftM = function leftShiftM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) << matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.leftShift = function leftShift(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.leftShift(value); + }; + + AbstractMatrix.prototype.signPropagatingRightShift = function signPropagatingRightShift(value) { + if (typeof value === 'number') return this.signPropagatingRightShiftS(value); + return this.signPropagatingRightShiftM(value); + }; + + AbstractMatrix.prototype.signPropagatingRightShiftS = function signPropagatingRightShiftS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) >> value); + } + } + return this; + }; + + AbstractMatrix.prototype.signPropagatingRightShiftM = function signPropagatingRightShiftM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) >> matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.signPropagatingRightShift = function signPropagatingRightShift(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.signPropagatingRightShift(value); + }; + + AbstractMatrix.prototype.rightShift = function rightShift(value) { + if (typeof value === 'number') return this.rightShiftS(value); + return this.rightShiftM(value); + }; + + AbstractMatrix.prototype.rightShiftS = function rightShiftS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) >>> value); + } + } + return this; + }; + + AbstractMatrix.prototype.rightShiftM = function rightShiftM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) >>> matrix.get(i, j)); + } + } + return this; + }; + + AbstractMatrix.rightShift = function rightShift(matrix, value) { + const newMatrix = new Matrix(matrix); + return newMatrix.rightShift(value); + }; + AbstractMatrix.prototype.zeroFillRightShift = AbstractMatrix.prototype.rightShift; + AbstractMatrix.prototype.zeroFillRightShiftS = AbstractMatrix.prototype.rightShiftS; + AbstractMatrix.prototype.zeroFillRightShiftM = AbstractMatrix.prototype.rightShiftM; + AbstractMatrix.zeroFillRightShift = AbstractMatrix.rightShift; + + AbstractMatrix.prototype.not = function not() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, ~(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.not = function not(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.not(); + }; + + AbstractMatrix.prototype.abs = function abs() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.abs(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.abs = function abs(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.abs(); + }; + + AbstractMatrix.prototype.acos = function acos() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.acos(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.acos = function acos(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.acos(); + }; + + AbstractMatrix.prototype.acosh = function acosh() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.acosh(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.acosh = function acosh(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.acosh(); + }; + + AbstractMatrix.prototype.asin = function asin() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.asin(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.asin = function asin(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.asin(); + }; + + AbstractMatrix.prototype.asinh = function asinh() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.asinh(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.asinh = function asinh(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.asinh(); + }; + + AbstractMatrix.prototype.atan = function atan() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.atan(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.atan = function atan(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.atan(); + }; + + AbstractMatrix.prototype.atanh = function atanh() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.atanh(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.atanh = function atanh(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.atanh(); + }; + + AbstractMatrix.prototype.cbrt = function cbrt() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.cbrt(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.cbrt = function cbrt(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.cbrt(); + }; + + AbstractMatrix.prototype.ceil = function ceil() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.ceil(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.ceil = function ceil(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.ceil(); + }; + + AbstractMatrix.prototype.clz32 = function clz32() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.clz32(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.clz32 = function clz32(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.clz32(); + }; + + AbstractMatrix.prototype.cos = function cos() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.cos(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.cos = function cos(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.cos(); + }; + + AbstractMatrix.prototype.cosh = function cosh() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.cosh(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.cosh = function cosh(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.cosh(); + }; + + AbstractMatrix.prototype.exp = function exp() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.exp(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.exp = function exp(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.exp(); + }; + + AbstractMatrix.prototype.expm1 = function expm1() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.expm1(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.expm1 = function expm1(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.expm1(); + }; + + AbstractMatrix.prototype.floor = function floor() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.floor(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.floor = function floor(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.floor(); + }; + + AbstractMatrix.prototype.fround = function fround() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.fround(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.fround = function fround(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.fround(); + }; + + AbstractMatrix.prototype.log = function log() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.log(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.log = function log(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.log(); + }; + + AbstractMatrix.prototype.log1p = function log1p() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.log1p(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.log1p = function log1p(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.log1p(); + }; + + AbstractMatrix.prototype.log10 = function log10() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.log10(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.log10 = function log10(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.log10(); + }; + + AbstractMatrix.prototype.log2 = function log2() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.log2(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.log2 = function log2(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.log2(); + }; + + AbstractMatrix.prototype.round = function round() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.round(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.round = function round(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.round(); + }; + + AbstractMatrix.prototype.sign = function sign() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.sign(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.sign = function sign(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.sign(); + }; + + AbstractMatrix.prototype.sin = function sin() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.sin(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.sin = function sin(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.sin(); + }; + + AbstractMatrix.prototype.sinh = function sinh() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.sinh(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.sinh = function sinh(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.sinh(); + }; + + AbstractMatrix.prototype.sqrt = function sqrt() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.sqrt(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.sqrt = function sqrt(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.sqrt(); + }; + + AbstractMatrix.prototype.tan = function tan() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.tan(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.tan = function tan(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.tan(); + }; + + AbstractMatrix.prototype.tanh = function tanh() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.tanh(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.tanh = function tanh(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.tanh(); + }; + + AbstractMatrix.prototype.trunc = function trunc() { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, Math.trunc(this.get(i, j))); + } + } + return this; + }; + + AbstractMatrix.trunc = function trunc(matrix) { + const newMatrix = new Matrix(matrix); + return newMatrix.trunc(); + }; + + AbstractMatrix.pow = function pow(matrix, arg0) { + const newMatrix = new Matrix(matrix); + return newMatrix.pow(arg0); + }; + + AbstractMatrix.prototype.pow = function pow(value) { + if (typeof value === 'number') return this.powS(value); + return this.powM(value); + }; + + AbstractMatrix.prototype.powS = function powS(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) ** value); + } + } + return this; + }; + + AbstractMatrix.prototype.powM = function powM(matrix) { + matrix = Matrix.checkMatrix(matrix); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) ** matrix.get(i, j)); + } + } + return this; + }; +} + +/** + * @private + * Check that a row index is not out of bounds + * @param {Matrix} matrix + * @param {number} index + * @param {boolean} [outer] + */ +function checkRowIndex(matrix, index, outer) { + let max = outer ? matrix.rows : matrix.rows - 1; + if (index < 0 || index > max) { + throw new RangeError('Row index out of range'); + } +} + +/** + * @private + * Check that a column index is not out of bounds + * @param {Matrix} matrix + * @param {number} index + * @param {boolean} [outer] + */ +function checkColumnIndex(matrix, index, outer) { + let max = outer ? matrix.columns : matrix.columns - 1; + if (index < 0 || index > max) { + throw new RangeError('Column index out of range'); + } +} + +/** + * @private + * Check that the provided vector is an array with the right length + * @param {Matrix} matrix + * @param {Array|Matrix} vector + * @return {Array} + * @throws {RangeError} + */ +function checkRowVector(matrix, vector) { + if (vector.to1DArray) { + vector = vector.to1DArray(); + } + if (vector.length !== matrix.columns) { + throw new RangeError( + 'vector size must be the same as the number of columns', + ); + } + return vector; +} + +/** + * @private + * Check that the provided vector is an array with the right length + * @param {Matrix} matrix + * @param {Array|Matrix} vector + * @return {Array} + * @throws {RangeError} + */ +function checkColumnVector(matrix, vector) { + if (vector.to1DArray) { + vector = vector.to1DArray(); + } + if (vector.length !== matrix.rows) { + throw new RangeError('vector size must be the same as the number of rows'); + } + return vector; +} + +function checkRowIndices(matrix, rowIndices) { + if (!isAnyArray.isAnyArray(rowIndices)) { + throw new TypeError('row indices must be an array'); + } + + for (let i = 0; i < rowIndices.length; i++) { + if (rowIndices[i] < 0 || rowIndices[i] >= matrix.rows) { + throw new RangeError('row indices are out of range'); + } + } +} + +function checkColumnIndices(matrix, columnIndices) { + if (!isAnyArray.isAnyArray(columnIndices)) { + throw new TypeError('column indices must be an array'); + } + + for (let i = 0; i < columnIndices.length; i++) { + if (columnIndices[i] < 0 || columnIndices[i] >= matrix.columns) { + throw new RangeError('column indices are out of range'); + } + } +} + +function checkRange(matrix, startRow, endRow, startColumn, endColumn) { + if (arguments.length !== 5) { + throw new RangeError('expected 4 arguments'); + } + checkNumber('startRow', startRow); + checkNumber('endRow', endRow); + checkNumber('startColumn', startColumn); + checkNumber('endColumn', endColumn); + if ( + startRow > endRow || + startColumn > endColumn || + startRow < 0 || + startRow >= matrix.rows || + endRow < 0 || + endRow >= matrix.rows || + startColumn < 0 || + startColumn >= matrix.columns || + endColumn < 0 || + endColumn >= matrix.columns + ) { + throw new RangeError('Submatrix indices are out of range'); + } +} + +function newArray(length, value = 0) { + let array = []; + for (let i = 0; i < length; i++) { + array.push(value); + } + return array; +} + +function checkNumber(name, value) { + if (typeof value !== 'number') { + throw new TypeError(`${name} must be a number`); + } +} + +function checkNonEmpty(matrix) { + if (matrix.isEmpty()) { + throw new Error('Empty matrix has no elements to index'); + } +} + +function sumByRow(matrix) { + let sum = newArray(matrix.rows); + for (let i = 0; i < matrix.rows; ++i) { + for (let j = 0; j < matrix.columns; ++j) { + sum[i] += matrix.get(i, j); + } + } + return sum; +} + +function sumByColumn(matrix) { + let sum = newArray(matrix.columns); + for (let i = 0; i < matrix.rows; ++i) { + for (let j = 0; j < matrix.columns; ++j) { + sum[j] += matrix.get(i, j); + } + } + return sum; +} + +function sumAll(matrix) { + let v = 0; + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + v += matrix.get(i, j); + } + } + return v; +} + +function productByRow(matrix) { + let sum = newArray(matrix.rows, 1); + for (let i = 0; i < matrix.rows; ++i) { + for (let j = 0; j < matrix.columns; ++j) { + sum[i] *= matrix.get(i, j); + } + } + return sum; +} + +function productByColumn(matrix) { + let sum = newArray(matrix.columns, 1); + for (let i = 0; i < matrix.rows; ++i) { + for (let j = 0; j < matrix.columns; ++j) { + sum[j] *= matrix.get(i, j); + } + } + return sum; +} + +function productAll(matrix) { + let v = 1; + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + v *= matrix.get(i, j); + } + } + return v; +} + +function varianceByRow(matrix, unbiased, mean) { + const rows = matrix.rows; + const cols = matrix.columns; + const variance = []; + + for (let i = 0; i < rows; i++) { + let sum1 = 0; + let sum2 = 0; + let x = 0; + for (let j = 0; j < cols; j++) { + x = matrix.get(i, j) - mean[i]; + sum1 += x; + sum2 += x * x; + } + if (unbiased) { + variance.push((sum2 - (sum1 * sum1) / cols) / (cols - 1)); + } else { + variance.push((sum2 - (sum1 * sum1) / cols) / cols); + } + } + return variance; +} + +function varianceByColumn(matrix, unbiased, mean) { + const rows = matrix.rows; + const cols = matrix.columns; + const variance = []; + + for (let j = 0; j < cols; j++) { + let sum1 = 0; + let sum2 = 0; + let x = 0; + for (let i = 0; i < rows; i++) { + x = matrix.get(i, j) - mean[j]; + sum1 += x; + sum2 += x * x; + } + if (unbiased) { + variance.push((sum2 - (sum1 * sum1) / rows) / (rows - 1)); + } else { + variance.push((sum2 - (sum1 * sum1) / rows) / rows); + } + } + return variance; +} + +function varianceAll(matrix, unbiased, mean) { + const rows = matrix.rows; + const cols = matrix.columns; + const size = rows * cols; + + let sum1 = 0; + let sum2 = 0; + let x = 0; + for (let i = 0; i < rows; i++) { + for (let j = 0; j < cols; j++) { + x = matrix.get(i, j) - mean; + sum1 += x; + sum2 += x * x; + } + } + if (unbiased) { + return (sum2 - (sum1 * sum1) / size) / (size - 1); + } else { + return (sum2 - (sum1 * sum1) / size) / size; + } +} + +function centerByRow(matrix, mean) { + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + matrix.set(i, j, matrix.get(i, j) - mean[i]); + } + } +} + +function centerByColumn(matrix, mean) { + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + matrix.set(i, j, matrix.get(i, j) - mean[j]); + } + } +} + +function centerAll(matrix, mean) { + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + matrix.set(i, j, matrix.get(i, j) - mean); + } + } +} + +function getScaleByRow(matrix) { + const scale = []; + for (let i = 0; i < matrix.rows; i++) { + let sum = 0; + for (let j = 0; j < matrix.columns; j++) { + sum += matrix.get(i, j) ** 2 / (matrix.columns - 1); + } + scale.push(Math.sqrt(sum)); + } + return scale; +} + +function scaleByRow(matrix, scale) { + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + matrix.set(i, j, matrix.get(i, j) / scale[i]); + } + } +} + +function getScaleByColumn(matrix) { + const scale = []; + for (let j = 0; j < matrix.columns; j++) { + let sum = 0; + for (let i = 0; i < matrix.rows; i++) { + sum += matrix.get(i, j) ** 2 / (matrix.rows - 1); + } + scale.push(Math.sqrt(sum)); + } + return scale; +} + +function scaleByColumn(matrix, scale) { + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + matrix.set(i, j, matrix.get(i, j) / scale[j]); + } + } +} + +function getScaleAll(matrix) { + const divider = matrix.size - 1; + let sum = 0; + for (let j = 0; j < matrix.columns; j++) { + for (let i = 0; i < matrix.rows; i++) { + sum += matrix.get(i, j) ** 2 / divider; + } + } + return Math.sqrt(sum); +} + +function scaleAll(matrix, scale) { + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + matrix.set(i, j, matrix.get(i, j) / scale); + } + } +} + +class AbstractMatrix { + static from1DArray(newRows, newColumns, newData) { + let length = newRows * newColumns; + if (length !== newData.length) { + throw new RangeError('data length does not match given dimensions'); + } + let newMatrix = new Matrix$1(newRows, newColumns); + for (let row = 0; row < newRows; row++) { + for (let column = 0; column < newColumns; column++) { + newMatrix.set(row, column, newData[row * newColumns + column]); + } + } + return newMatrix; + } + + static rowVector(newData) { + let vector = new Matrix$1(1, newData.length); + for (let i = 0; i < newData.length; i++) { + vector.set(0, i, newData[i]); + } + return vector; + } + + static columnVector(newData) { + let vector = new Matrix$1(newData.length, 1); + for (let i = 0; i < newData.length; i++) { + vector.set(i, 0, newData[i]); + } + return vector; + } + + static zeros(rows, columns) { + return new Matrix$1(rows, columns); + } + + static ones(rows, columns) { + return new Matrix$1(rows, columns).fill(1); + } + + static rand(rows, columns, options = {}) { + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + const { random = Math.random } = options; + let matrix = new Matrix$1(rows, columns); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + matrix.set(i, j, random()); + } + } + return matrix; + } + + static randInt(rows, columns, options = {}) { + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + const { min = 0, max = 1000, random = Math.random } = options; + if (!Number.isInteger(min)) throw new TypeError('min must be an integer'); + if (!Number.isInteger(max)) throw new TypeError('max must be an integer'); + if (min >= max) throw new RangeError('min must be smaller than max'); + let interval = max - min; + let matrix = new Matrix$1(rows, columns); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + let value = min + Math.round(random() * interval); + matrix.set(i, j, value); + } + } + return matrix; + } + + static eye(rows, columns, value) { + if (columns === undefined) columns = rows; + if (value === undefined) value = 1; + let min = Math.min(rows, columns); + let matrix = this.zeros(rows, columns); + for (let i = 0; i < min; i++) { + matrix.set(i, i, value); + } + return matrix; + } + + static diag(data, rows, columns) { + let l = data.length; + if (rows === undefined) rows = l; + if (columns === undefined) columns = rows; + let min = Math.min(l, rows, columns); + let matrix = this.zeros(rows, columns); + for (let i = 0; i < min; i++) { + matrix.set(i, i, data[i]); + } + return matrix; + } + + static min(matrix1, matrix2) { + matrix1 = this.checkMatrix(matrix1); + matrix2 = this.checkMatrix(matrix2); + let rows = matrix1.rows; + let columns = matrix1.columns; + let result = new Matrix$1(rows, columns); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + result.set(i, j, Math.min(matrix1.get(i, j), matrix2.get(i, j))); + } + } + return result; + } + + static max(matrix1, matrix2) { + matrix1 = this.checkMatrix(matrix1); + matrix2 = this.checkMatrix(matrix2); + let rows = matrix1.rows; + let columns = matrix1.columns; + let result = new this(rows, columns); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + result.set(i, j, Math.max(matrix1.get(i, j), matrix2.get(i, j))); + } + } + return result; + } + + static checkMatrix(value) { + return AbstractMatrix.isMatrix(value) ? value : new Matrix$1(value); + } + + static isMatrix(value) { + return value != null && value.klass === 'Matrix'; + } + + get size() { + return this.rows * this.columns; + } + + apply(callback) { + if (typeof callback !== 'function') { + throw new TypeError('callback must be a function'); + } + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + callback.call(this, i, j); + } + } + return this; + } + + to1DArray() { + let array = []; + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + array.push(this.get(i, j)); + } + } + return array; + } + + to2DArray() { + let copy = []; + for (let i = 0; i < this.rows; i++) { + copy.push([]); + for (let j = 0; j < this.columns; j++) { + copy[i].push(this.get(i, j)); + } + } + return copy; + } + + toJSON() { + return this.to2DArray(); + } + + isRowVector() { + return this.rows === 1; + } + + isColumnVector() { + return this.columns === 1; + } + + isVector() { + return this.rows === 1 || this.columns === 1; + } + + isSquare() { + return this.rows === this.columns; + } + + isEmpty() { + return this.rows === 0 || this.columns === 0; + } + + isSymmetric() { + if (this.isSquare()) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j <= i; j++) { + if (this.get(i, j) !== this.get(j, i)) { + return false; + } + } + } + return true; + } + return false; + } + + isDistance() { + if (!this.isSymmetric()) return false; + + for (let i = 0; i < this.rows; i++) { + if (this.get(i, i) !== 0) return false; + } + + return true; + } + + isEchelonForm() { + let i = 0; + let j = 0; + let previousColumn = -1; + let isEchelonForm = true; + let checked = false; + while (i < this.rows && isEchelonForm) { + j = 0; + checked = false; + while (j < this.columns && checked === false) { + if (this.get(i, j) === 0) { + j++; + } else if (this.get(i, j) === 1 && j > previousColumn) { + checked = true; + previousColumn = j; + } else { + isEchelonForm = false; + checked = true; + } + } + i++; + } + return isEchelonForm; + } + + isReducedEchelonForm() { + let i = 0; + let j = 0; + let previousColumn = -1; + let isReducedEchelonForm = true; + let checked = false; + while (i < this.rows && isReducedEchelonForm) { + j = 0; + checked = false; + while (j < this.columns && checked === false) { + if (this.get(i, j) === 0) { + j++; + } else if (this.get(i, j) === 1 && j > previousColumn) { + checked = true; + previousColumn = j; + } else { + isReducedEchelonForm = false; + checked = true; + } + } + for (let k = j + 1; k < this.rows; k++) { + if (this.get(i, k) !== 0) { + isReducedEchelonForm = false; + } + } + i++; + } + return isReducedEchelonForm; + } + + echelonForm() { + let result = this.clone(); + let h = 0; + let k = 0; + while (h < result.rows && k < result.columns) { + let iMax = h; + for (let i = h; i < result.rows; i++) { + if (result.get(i, k) > result.get(iMax, k)) { + iMax = i; + } + } + if (result.get(iMax, k) === 0) { + k++; + } else { + result.swapRows(h, iMax); + let tmp = result.get(h, k); + for (let j = k; j < result.columns; j++) { + result.set(h, j, result.get(h, j) / tmp); + } + for (let i = h + 1; i < result.rows; i++) { + let factor = result.get(i, k) / result.get(h, k); + result.set(i, k, 0); + for (let j = k + 1; j < result.columns; j++) { + result.set(i, j, result.get(i, j) - result.get(h, j) * factor); + } + } + h++; + k++; + } + } + return result; + } + + reducedEchelonForm() { + let result = this.echelonForm(); + let m = result.columns; + let n = result.rows; + let h = n - 1; + while (h >= 0) { + if (result.maxRow(h) === 0) { + h--; + } else { + let p = 0; + let pivot = false; + while (p < n && pivot === false) { + if (result.get(h, p) === 1) { + pivot = true; + } else { + p++; + } + } + for (let i = 0; i < h; i++) { + let factor = result.get(i, p); + for (let j = p; j < m; j++) { + let tmp = result.get(i, j) - factor * result.get(h, j); + result.set(i, j, tmp); + } + } + h--; + } + } + return result; + } + + set() { + throw new Error('set method is unimplemented'); + } + + get() { + throw new Error('get method is unimplemented'); + } + + repeat(options = {}) { + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + const { rows = 1, columns = 1 } = options; + if (!Number.isInteger(rows) || rows <= 0) { + throw new TypeError('rows must be a positive integer'); + } + if (!Number.isInteger(columns) || columns <= 0) { + throw new TypeError('columns must be a positive integer'); + } + let matrix = new Matrix$1(this.rows * rows, this.columns * columns); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + matrix.setSubMatrix(this, this.rows * i, this.columns * j); + } + } + return matrix; + } + + fill(value) { + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, value); + } + } + return this; + } + + neg() { + return this.mulS(-1); + } + + getRow(index) { + checkRowIndex(this, index); + let row = []; + for (let i = 0; i < this.columns; i++) { + row.push(this.get(index, i)); + } + return row; + } + + getRowVector(index) { + return Matrix$1.rowVector(this.getRow(index)); + } + + setRow(index, array) { + checkRowIndex(this, index); + array = checkRowVector(this, array); + for (let i = 0; i < this.columns; i++) { + this.set(index, i, array[i]); + } + return this; + } + + swapRows(row1, row2) { + checkRowIndex(this, row1); + checkRowIndex(this, row2); + for (let i = 0; i < this.columns; i++) { + let temp = this.get(row1, i); + this.set(row1, i, this.get(row2, i)); + this.set(row2, i, temp); + } + return this; + } + + getColumn(index) { + checkColumnIndex(this, index); + let column = []; + for (let i = 0; i < this.rows; i++) { + column.push(this.get(i, index)); + } + return column; + } + + getColumnVector(index) { + return Matrix$1.columnVector(this.getColumn(index)); + } + + setColumn(index, array) { + checkColumnIndex(this, index); + array = checkColumnVector(this, array); + for (let i = 0; i < this.rows; i++) { + this.set(i, index, array[i]); + } + return this; + } + + swapColumns(column1, column2) { + checkColumnIndex(this, column1); + checkColumnIndex(this, column2); + for (let i = 0; i < this.rows; i++) { + let temp = this.get(i, column1); + this.set(i, column1, this.get(i, column2)); + this.set(i, column2, temp); + } + return this; + } + + addRowVector(vector) { + vector = checkRowVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) + vector[j]); + } + } + return this; + } + + subRowVector(vector) { + vector = checkRowVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) - vector[j]); + } + } + return this; + } + + mulRowVector(vector) { + vector = checkRowVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) * vector[j]); + } + } + return this; + } + + divRowVector(vector) { + vector = checkRowVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) / vector[j]); + } + } + return this; + } + + addColumnVector(vector) { + vector = checkColumnVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) + vector[i]); + } + } + return this; + } + + subColumnVector(vector) { + vector = checkColumnVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) - vector[i]); + } + } + return this; + } + + mulColumnVector(vector) { + vector = checkColumnVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) * vector[i]); + } + } + return this; + } + + divColumnVector(vector) { + vector = checkColumnVector(this, vector); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + this.set(i, j, this.get(i, j) / vector[i]); + } + } + return this; + } + + mulRow(index, value) { + checkRowIndex(this, index); + for (let i = 0; i < this.columns; i++) { + this.set(index, i, this.get(index, i) * value); + } + return this; + } + + mulColumn(index, value) { + checkColumnIndex(this, index); + for (let i = 0; i < this.rows; i++) { + this.set(i, index, this.get(i, index) * value); + } + return this; + } + + max(by) { + if (this.isEmpty()) { + return NaN; + } + switch (by) { + case 'row': { + const max = new Array(this.rows).fill(Number.NEGATIVE_INFINITY); + for (let row = 0; row < this.rows; row++) { + for (let column = 0; column < this.columns; column++) { + if (this.get(row, column) > max[row]) { + max[row] = this.get(row, column); + } + } + } + return max; + } + case 'column': { + const max = new Array(this.columns).fill(Number.NEGATIVE_INFINITY); + for (let row = 0; row < this.rows; row++) { + for (let column = 0; column < this.columns; column++) { + if (this.get(row, column) > max[column]) { + max[column] = this.get(row, column); + } + } + } + return max; + } + case undefined: { + let max = this.get(0, 0); + for (let row = 0; row < this.rows; row++) { + for (let column = 0; column < this.columns; column++) { + if (this.get(row, column) > max) { + max = this.get(row, column); + } + } + } + return max; + } + default: + throw new Error(`invalid option: ${by}`); + } + } + + maxIndex() { + checkNonEmpty(this); + let v = this.get(0, 0); + let idx = [0, 0]; + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + if (this.get(i, j) > v) { + v = this.get(i, j); + idx[0] = i; + idx[1] = j; + } + } + } + return idx; + } + + min(by) { + if (this.isEmpty()) { + return NaN; + } + + switch (by) { + case 'row': { + const min = new Array(this.rows).fill(Number.POSITIVE_INFINITY); + for (let row = 0; row < this.rows; row++) { + for (let column = 0; column < this.columns; column++) { + if (this.get(row, column) < min[row]) { + min[row] = this.get(row, column); + } + } + } + return min; + } + case 'column': { + const min = new Array(this.columns).fill(Number.POSITIVE_INFINITY); + for (let row = 0; row < this.rows; row++) { + for (let column = 0; column < this.columns; column++) { + if (this.get(row, column) < min[column]) { + min[column] = this.get(row, column); + } + } + } + return min; + } + case undefined: { + let min = this.get(0, 0); + for (let row = 0; row < this.rows; row++) { + for (let column = 0; column < this.columns; column++) { + if (this.get(row, column) < min) { + min = this.get(row, column); + } + } + } + return min; + } + default: + throw new Error(`invalid option: ${by}`); + } + } + + minIndex() { + checkNonEmpty(this); + let v = this.get(0, 0); + let idx = [0, 0]; + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + if (this.get(i, j) < v) { + v = this.get(i, j); + idx[0] = i; + idx[1] = j; + } + } + } + return idx; + } + + maxRow(row) { + checkRowIndex(this, row); + if (this.isEmpty()) { + return NaN; + } + let v = this.get(row, 0); + for (let i = 1; i < this.columns; i++) { + if (this.get(row, i) > v) { + v = this.get(row, i); + } + } + return v; + } + + maxRowIndex(row) { + checkRowIndex(this, row); + checkNonEmpty(this); + let v = this.get(row, 0); + let idx = [row, 0]; + for (let i = 1; i < this.columns; i++) { + if (this.get(row, i) > v) { + v = this.get(row, i); + idx[1] = i; + } + } + return idx; + } + + minRow(row) { + checkRowIndex(this, row); + if (this.isEmpty()) { + return NaN; + } + let v = this.get(row, 0); + for (let i = 1; i < this.columns; i++) { + if (this.get(row, i) < v) { + v = this.get(row, i); + } + } + return v; + } + + minRowIndex(row) { + checkRowIndex(this, row); + checkNonEmpty(this); + let v = this.get(row, 0); + let idx = [row, 0]; + for (let i = 1; i < this.columns; i++) { + if (this.get(row, i) < v) { + v = this.get(row, i); + idx[1] = i; + } + } + return idx; + } + + maxColumn(column) { + checkColumnIndex(this, column); + if (this.isEmpty()) { + return NaN; + } + let v = this.get(0, column); + for (let i = 1; i < this.rows; i++) { + if (this.get(i, column) > v) { + v = this.get(i, column); + } + } + return v; + } + + maxColumnIndex(column) { + checkColumnIndex(this, column); + checkNonEmpty(this); + let v = this.get(0, column); + let idx = [0, column]; + for (let i = 1; i < this.rows; i++) { + if (this.get(i, column) > v) { + v = this.get(i, column); + idx[0] = i; + } + } + return idx; + } + + minColumn(column) { + checkColumnIndex(this, column); + if (this.isEmpty()) { + return NaN; + } + let v = this.get(0, column); + for (let i = 1; i < this.rows; i++) { + if (this.get(i, column) < v) { + v = this.get(i, column); + } + } + return v; + } + + minColumnIndex(column) { + checkColumnIndex(this, column); + checkNonEmpty(this); + let v = this.get(0, column); + let idx = [0, column]; + for (let i = 1; i < this.rows; i++) { + if (this.get(i, column) < v) { + v = this.get(i, column); + idx[0] = i; + } + } + return idx; + } + + diag() { + let min = Math.min(this.rows, this.columns); + let diag = []; + for (let i = 0; i < min; i++) { + diag.push(this.get(i, i)); + } + return diag; + } + + norm(type = 'frobenius') { + switch (type) { + case 'max': + return this.max(); + case 'frobenius': + return Math.sqrt(this.dot(this)); + default: + throw new RangeError(`unknown norm type: ${type}`); + } + } + + cumulativeSum() { + let sum = 0; + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + sum += this.get(i, j); + this.set(i, j, sum); + } + } + return this; + } + + dot(vector2) { + if (AbstractMatrix.isMatrix(vector2)) vector2 = vector2.to1DArray(); + let vector1 = this.to1DArray(); + if (vector1.length !== vector2.length) { + throw new RangeError('vectors do not have the same size'); + } + let dot = 0; + for (let i = 0; i < vector1.length; i++) { + dot += vector1[i] * vector2[i]; + } + return dot; + } + + mmul(other) { + other = Matrix$1.checkMatrix(other); + + let m = this.rows; + let n = this.columns; + let p = other.columns; + + let result = new Matrix$1(m, p); + + let Bcolj = new Float64Array(n); + for (let j = 0; j < p; j++) { + for (let k = 0; k < n; k++) { + Bcolj[k] = other.get(k, j); + } + + for (let i = 0; i < m; i++) { + let s = 0; + for (let k = 0; k < n; k++) { + s += this.get(i, k) * Bcolj[k]; + } + + result.set(i, j, s); + } + } + return result; + } + + mpow(scalar) { + if (!this.isSquare()) { + throw new RangeError('Matrix must be square'); + } + if (!Number.isInteger(scalar) || scalar < 0) { + throw new RangeError('Exponent must be a non-negative integer'); + } + // Russian Peasant exponentiation, i.e. exponentiation by squaring + let result = Matrix$1.eye(this.rows); + let bb = this; + // Note: Don't bit shift. In JS, that would truncate at 32 bits + for (let e = scalar; e > 1; e /= 2) { + if ((e & 1) !== 0) { + result = result.mmul(bb); + } + bb = bb.mmul(bb); + } + return result; + } + + strassen2x2(other) { + other = Matrix$1.checkMatrix(other); + let result = new Matrix$1(2, 2); + const a11 = this.get(0, 0); + const b11 = other.get(0, 0); + const a12 = this.get(0, 1); + const b12 = other.get(0, 1); + const a21 = this.get(1, 0); + const b21 = other.get(1, 0); + const a22 = this.get(1, 1); + const b22 = other.get(1, 1); + + // Compute intermediate values. + const m1 = (a11 + a22) * (b11 + b22); + const m2 = (a21 + a22) * b11; + const m3 = a11 * (b12 - b22); + const m4 = a22 * (b21 - b11); + const m5 = (a11 + a12) * b22; + const m6 = (a21 - a11) * (b11 + b12); + const m7 = (a12 - a22) * (b21 + b22); + + // Combine intermediate values into the output. + const c00 = m1 + m4 - m5 + m7; + const c01 = m3 + m5; + const c10 = m2 + m4; + const c11 = m1 - m2 + m3 + m6; + + result.set(0, 0, c00); + result.set(0, 1, c01); + result.set(1, 0, c10); + result.set(1, 1, c11); + return result; + } + + strassen3x3(other) { + other = Matrix$1.checkMatrix(other); + let result = new Matrix$1(3, 3); + + const a00 = this.get(0, 0); + const a01 = this.get(0, 1); + const a02 = this.get(0, 2); + const a10 = this.get(1, 0); + const a11 = this.get(1, 1); + const a12 = this.get(1, 2); + const a20 = this.get(2, 0); + const a21 = this.get(2, 1); + const a22 = this.get(2, 2); + + const b00 = other.get(0, 0); + const b01 = other.get(0, 1); + const b02 = other.get(0, 2); + const b10 = other.get(1, 0); + const b11 = other.get(1, 1); + const b12 = other.get(1, 2); + const b20 = other.get(2, 0); + const b21 = other.get(2, 1); + const b22 = other.get(2, 2); + + const m1 = (a00 + a01 + a02 - a10 - a11 - a21 - a22) * b11; + const m2 = (a00 - a10) * (-b01 + b11); + const m3 = a11 * (-b00 + b01 + b10 - b11 - b12 - b20 + b22); + const m4 = (-a00 + a10 + a11) * (b00 - b01 + b11); + const m5 = (a10 + a11) * (-b00 + b01); + const m6 = a00 * b00; + const m7 = (-a00 + a20 + a21) * (b00 - b02 + b12); + const m8 = (-a00 + a20) * (b02 - b12); + const m9 = (a20 + a21) * (-b00 + b02); + const m10 = (a00 + a01 + a02 - a11 - a12 - a20 - a21) * b12; + const m11 = a21 * (-b00 + b02 + b10 - b11 - b12 - b20 + b21); + const m12 = (-a02 + a21 + a22) * (b11 + b20 - b21); + const m13 = (a02 - a22) * (b11 - b21); + const m14 = a02 * b20; + const m15 = (a21 + a22) * (-b20 + b21); + const m16 = (-a02 + a11 + a12) * (b12 + b20 - b22); + const m17 = (a02 - a12) * (b12 - b22); + const m18 = (a11 + a12) * (-b20 + b22); + const m19 = a01 * b10; + const m20 = a12 * b21; + const m21 = a10 * b02; + const m22 = a20 * b01; + const m23 = a22 * b22; + + const c00 = m6 + m14 + m19; + const c01 = m1 + m4 + m5 + m6 + m12 + m14 + m15; + const c02 = m6 + m7 + m9 + m10 + m14 + m16 + m18; + const c10 = m2 + m3 + m4 + m6 + m14 + m16 + m17; + const c11 = m2 + m4 + m5 + m6 + m20; + const c12 = m14 + m16 + m17 + m18 + m21; + const c20 = m6 + m7 + m8 + m11 + m12 + m13 + m14; + const c21 = m12 + m13 + m14 + m15 + m22; + const c22 = m6 + m7 + m8 + m9 + m23; + + result.set(0, 0, c00); + result.set(0, 1, c01); + result.set(0, 2, c02); + result.set(1, 0, c10); + result.set(1, 1, c11); + result.set(1, 2, c12); + result.set(2, 0, c20); + result.set(2, 1, c21); + result.set(2, 2, c22); + return result; + } + + mmulStrassen(y) { + y = Matrix$1.checkMatrix(y); + let x = this.clone(); + let r1 = x.rows; + let c1 = x.columns; + let r2 = y.rows; + let c2 = y.columns; + if (c1 !== r2) { + // eslint-disable-next-line no-console + console.warn( + `Multiplying ${r1} x ${c1} and ${r2} x ${c2} matrix: dimensions do not match.`, + ); + } + + // Put a matrix into the top left of a matrix of zeros. + // `rows` and `cols` are the dimensions of the output matrix. + function embed(mat, rows, cols) { + let r = mat.rows; + let c = mat.columns; + if (r === rows && c === cols) { + return mat; + } else { + let resultat = AbstractMatrix.zeros(rows, cols); + resultat = resultat.setSubMatrix(mat, 0, 0); + return resultat; + } + } + + // Make sure both matrices are the same size. + // This is exclusively for simplicity: + // this algorithm can be implemented with matrices of different sizes. + + let r = Math.max(r1, r2); + let c = Math.max(c1, c2); + x = embed(x, r, c); + y = embed(y, r, c); + + // Our recursive multiplication function. + function blockMult(a, b, rows, cols) { + // For small matrices, resort to naive multiplication. + if (rows <= 512 || cols <= 512) { + return a.mmul(b); // a is equivalent to this + } + + // Apply dynamic padding. + if (rows % 2 === 1 && cols % 2 === 1) { + a = embed(a, rows + 1, cols + 1); + b = embed(b, rows + 1, cols + 1); + } else if (rows % 2 === 1) { + a = embed(a, rows + 1, cols); + b = embed(b, rows + 1, cols); + } else if (cols % 2 === 1) { + a = embed(a, rows, cols + 1); + b = embed(b, rows, cols + 1); + } + + let halfRows = parseInt(a.rows / 2, 10); + let halfCols = parseInt(a.columns / 2, 10); + // Subdivide input matrices. + let a11 = a.subMatrix(0, halfRows - 1, 0, halfCols - 1); + let b11 = b.subMatrix(0, halfRows - 1, 0, halfCols - 1); + + let a12 = a.subMatrix(0, halfRows - 1, halfCols, a.columns - 1); + let b12 = b.subMatrix(0, halfRows - 1, halfCols, b.columns - 1); + + let a21 = a.subMatrix(halfRows, a.rows - 1, 0, halfCols - 1); + let b21 = b.subMatrix(halfRows, b.rows - 1, 0, halfCols - 1); + + let a22 = a.subMatrix(halfRows, a.rows - 1, halfCols, a.columns - 1); + let b22 = b.subMatrix(halfRows, b.rows - 1, halfCols, b.columns - 1); + + // Compute intermediate values. + let m1 = blockMult( + AbstractMatrix.add(a11, a22), + AbstractMatrix.add(b11, b22), + halfRows, + halfCols, + ); + let m2 = blockMult(AbstractMatrix.add(a21, a22), b11, halfRows, halfCols); + let m3 = blockMult(a11, AbstractMatrix.sub(b12, b22), halfRows, halfCols); + let m4 = blockMult(a22, AbstractMatrix.sub(b21, b11), halfRows, halfCols); + let m5 = blockMult(AbstractMatrix.add(a11, a12), b22, halfRows, halfCols); + let m6 = blockMult( + AbstractMatrix.sub(a21, a11), + AbstractMatrix.add(b11, b12), + halfRows, + halfCols, + ); + let m7 = blockMult( + AbstractMatrix.sub(a12, a22), + AbstractMatrix.add(b21, b22), + halfRows, + halfCols, + ); + + // Combine intermediate values into the output. + let c11 = AbstractMatrix.add(m1, m4); + c11.sub(m5); + c11.add(m7); + let c12 = AbstractMatrix.add(m3, m5); + let c21 = AbstractMatrix.add(m2, m4); + let c22 = AbstractMatrix.sub(m1, m2); + c22.add(m3); + c22.add(m6); + + // Crop output to the desired size (undo dynamic padding). + let result = AbstractMatrix.zeros(2 * c11.rows, 2 * c11.columns); + result = result.setSubMatrix(c11, 0, 0); + result = result.setSubMatrix(c12, c11.rows, 0); + result = result.setSubMatrix(c21, 0, c11.columns); + result = result.setSubMatrix(c22, c11.rows, c11.columns); + return result.subMatrix(0, rows - 1, 0, cols - 1); + } + + return blockMult(x, y, r, c); + } + + scaleRows(options = {}) { + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + const { min = 0, max = 1 } = options; + if (!Number.isFinite(min)) throw new TypeError('min must be a number'); + if (!Number.isFinite(max)) throw new TypeError('max must be a number'); + if (min >= max) throw new RangeError('min must be smaller than max'); + let newMatrix = new Matrix$1(this.rows, this.columns); + for (let i = 0; i < this.rows; i++) { + const row = this.getRow(i); + if (row.length > 0) { + rescale(row, { min, max, output: row }); + } + newMatrix.setRow(i, row); + } + return newMatrix; + } + + scaleColumns(options = {}) { + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + const { min = 0, max = 1 } = options; + if (!Number.isFinite(min)) throw new TypeError('min must be a number'); + if (!Number.isFinite(max)) throw new TypeError('max must be a number'); + if (min >= max) throw new RangeError('min must be smaller than max'); + let newMatrix = new Matrix$1(this.rows, this.columns); + for (let i = 0; i < this.columns; i++) { + const column = this.getColumn(i); + if (column.length) { + rescale(column, { + min, + max, + output: column, + }); + } + newMatrix.setColumn(i, column); + } + return newMatrix; + } + + flipRows() { + const middle = Math.ceil(this.columns / 2); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < middle; j++) { + let first = this.get(i, j); + let last = this.get(i, this.columns - 1 - j); + this.set(i, j, last); + this.set(i, this.columns - 1 - j, first); + } + } + return this; + } + + flipColumns() { + const middle = Math.ceil(this.rows / 2); + for (let j = 0; j < this.columns; j++) { + for (let i = 0; i < middle; i++) { + let first = this.get(i, j); + let last = this.get(this.rows - 1 - i, j); + this.set(i, j, last); + this.set(this.rows - 1 - i, j, first); + } + } + return this; + } + + kroneckerProduct(other) { + other = Matrix$1.checkMatrix(other); + + let m = this.rows; + let n = this.columns; + let p = other.rows; + let q = other.columns; + + let result = new Matrix$1(m * p, n * q); + for (let i = 0; i < m; i++) { + for (let j = 0; j < n; j++) { + for (let k = 0; k < p; k++) { + for (let l = 0; l < q; l++) { + result.set(p * i + k, q * j + l, this.get(i, j) * other.get(k, l)); + } + } + } + } + return result; + } + + kroneckerSum(other) { + other = Matrix$1.checkMatrix(other); + if (!this.isSquare() || !other.isSquare()) { + throw new Error('Kronecker Sum needs two Square Matrices'); + } + let m = this.rows; + let n = other.rows; + let AxI = this.kroneckerProduct(Matrix$1.eye(n, n)); + let IxB = Matrix$1.eye(m, m).kroneckerProduct(other); + return AxI.add(IxB); + } + + transpose() { + let result = new Matrix$1(this.columns, this.rows); + for (let i = 0; i < this.rows; i++) { + for (let j = 0; j < this.columns; j++) { + result.set(j, i, this.get(i, j)); + } + } + return result; + } + + sortRows(compareFunction = compareNumbers) { + for (let i = 0; i < this.rows; i++) { + this.setRow(i, this.getRow(i).sort(compareFunction)); + } + return this; + } + + sortColumns(compareFunction = compareNumbers) { + for (let i = 0; i < this.columns; i++) { + this.setColumn(i, this.getColumn(i).sort(compareFunction)); + } + return this; + } + + subMatrix(startRow, endRow, startColumn, endColumn) { + checkRange(this, startRow, endRow, startColumn, endColumn); + let newMatrix = new Matrix$1( + endRow - startRow + 1, + endColumn - startColumn + 1, + ); + for (let i = startRow; i <= endRow; i++) { + for (let j = startColumn; j <= endColumn; j++) { + newMatrix.set(i - startRow, j - startColumn, this.get(i, j)); + } + } + return newMatrix; + } + + subMatrixRow(indices, startColumn, endColumn) { + if (startColumn === undefined) startColumn = 0; + if (endColumn === undefined) endColumn = this.columns - 1; + if ( + startColumn > endColumn || + startColumn < 0 || + startColumn >= this.columns || + endColumn < 0 || + endColumn >= this.columns + ) { + throw new RangeError('Argument out of range'); + } + + let newMatrix = new Matrix$1(indices.length, endColumn - startColumn + 1); + for (let i = 0; i < indices.length; i++) { + for (let j = startColumn; j <= endColumn; j++) { + if (indices[i] < 0 || indices[i] >= this.rows) { + throw new RangeError(`Row index out of range: ${indices[i]}`); + } + newMatrix.set(i, j - startColumn, this.get(indices[i], j)); + } + } + return newMatrix; + } + + subMatrixColumn(indices, startRow, endRow) { + if (startRow === undefined) startRow = 0; + if (endRow === undefined) endRow = this.rows - 1; + if ( + startRow > endRow || + startRow < 0 || + startRow >= this.rows || + endRow < 0 || + endRow >= this.rows + ) { + throw new RangeError('Argument out of range'); + } + + let newMatrix = new Matrix$1(endRow - startRow + 1, indices.length); + for (let i = 0; i < indices.length; i++) { + for (let j = startRow; j <= endRow; j++) { + if (indices[i] < 0 || indices[i] >= this.columns) { + throw new RangeError(`Column index out of range: ${indices[i]}`); + } + newMatrix.set(j - startRow, i, this.get(j, indices[i])); + } + } + return newMatrix; + } + + setSubMatrix(matrix, startRow, startColumn) { + matrix = Matrix$1.checkMatrix(matrix); + if (matrix.isEmpty()) { + return this; + } + let endRow = startRow + matrix.rows - 1; + let endColumn = startColumn + matrix.columns - 1; + checkRange(this, startRow, endRow, startColumn, endColumn); + for (let i = 0; i < matrix.rows; i++) { + for (let j = 0; j < matrix.columns; j++) { + this.set(startRow + i, startColumn + j, matrix.get(i, j)); + } + } + return this; + } + + selection(rowIndices, columnIndices) { + checkRowIndices(this, rowIndices); + checkColumnIndices(this, columnIndices); + let newMatrix = new Matrix$1(rowIndices.length, columnIndices.length); + for (let i = 0; i < rowIndices.length; i++) { + let rowIndex = rowIndices[i]; + for (let j = 0; j < columnIndices.length; j++) { + let columnIndex = columnIndices[j]; + newMatrix.set(i, j, this.get(rowIndex, columnIndex)); + } + } + return newMatrix; + } + + trace() { + let min = Math.min(this.rows, this.columns); + let trace = 0; + for (let i = 0; i < min; i++) { + trace += this.get(i, i); + } + return trace; + } + + clone() { + return this.constructor.copy(this, new Matrix$1(this.rows, this.columns)); + } + + /** + * @template {AbstractMatrix} M + * @param {AbstractMatrix} from + * @param {M} to + * @return {M} + */ + static copy(from, to) { + for (const [row, column, value] of from.entries()) { + to.set(row, column, value); + } + + return to; + } + + sum(by) { + switch (by) { + case 'row': + return sumByRow(this); + case 'column': + return sumByColumn(this); + case undefined: + return sumAll(this); + default: + throw new Error(`invalid option: ${by}`); + } + } + + product(by) { + switch (by) { + case 'row': + return productByRow(this); + case 'column': + return productByColumn(this); + case undefined: + return productAll(this); + default: + throw new Error(`invalid option: ${by}`); + } + } + + mean(by) { + const sum = this.sum(by); + switch (by) { + case 'row': { + for (let i = 0; i < this.rows; i++) { + sum[i] /= this.columns; + } + return sum; + } + case 'column': { + for (let i = 0; i < this.columns; i++) { + sum[i] /= this.rows; + } + return sum; + } + case undefined: + return sum / this.size; + default: + throw new Error(`invalid option: ${by}`); + } + } + + variance(by, options = {}) { + if (typeof by === 'object') { + options = by; + by = undefined; + } + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + const { unbiased = true, mean = this.mean(by) } = options; + if (typeof unbiased !== 'boolean') { + throw new TypeError('unbiased must be a boolean'); + } + switch (by) { + case 'row': { + if (!isAnyArray.isAnyArray(mean)) { + throw new TypeError('mean must be an array'); + } + return varianceByRow(this, unbiased, mean); + } + case 'column': { + if (!isAnyArray.isAnyArray(mean)) { + throw new TypeError('mean must be an array'); + } + return varianceByColumn(this, unbiased, mean); + } + case undefined: { + if (typeof mean !== 'number') { + throw new TypeError('mean must be a number'); + } + return varianceAll(this, unbiased, mean); + } + default: + throw new Error(`invalid option: ${by}`); + } + } + + standardDeviation(by, options) { + if (typeof by === 'object') { + options = by; + by = undefined; + } + const variance = this.variance(by, options); + if (by === undefined) { + return Math.sqrt(variance); + } else { + for (let i = 0; i < variance.length; i++) { + variance[i] = Math.sqrt(variance[i]); + } + return variance; + } + } + + center(by, options = {}) { + if (typeof by === 'object') { + options = by; + by = undefined; + } + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + const { center = this.mean(by) } = options; + switch (by) { + case 'row': { + if (!isAnyArray.isAnyArray(center)) { + throw new TypeError('center must be an array'); + } + centerByRow(this, center); + return this; + } + case 'column': { + if (!isAnyArray.isAnyArray(center)) { + throw new TypeError('center must be an array'); + } + centerByColumn(this, center); + return this; + } + case undefined: { + if (typeof center !== 'number') { + throw new TypeError('center must be a number'); + } + centerAll(this, center); + return this; + } + default: + throw new Error(`invalid option: ${by}`); + } + } + + scale(by, options = {}) { + if (typeof by === 'object') { + options = by; + by = undefined; + } + if (typeof options !== 'object') { + throw new TypeError('options must be an object'); + } + let scale = options.scale; + switch (by) { + case 'row': { + if (scale === undefined) { + scale = getScaleByRow(this); + } else if (!isAnyArray.isAnyArray(scale)) { + throw new TypeError('scale must be an array'); + } + scaleByRow(this, scale); + return this; + } + case 'column': { + if (scale === undefined) { + scale = getScaleByColumn(this); + } else if (!isAnyArray.isAnyArray(scale)) { + throw new TypeError('scale must be an array'); + } + scaleByColumn(this, scale); + return this; + } + case undefined: { + if (scale === undefined) { + scale = getScaleAll(this); + } else if (typeof scale !== 'number') { + throw new TypeError('scale must be a number'); + } + scaleAll(this, scale); + return this; + } + default: + throw new Error(`invalid option: ${by}`); + } + } + + toString(options) { + return inspectMatrixWithOptions(this, options); + } + + [Symbol.iterator]() { + return this.entries(); + } + + /** + * iterator from left to right, from top to bottom + * yield [row, column, value] + * @returns {Generator<[number, number, number], void, void>} + */ + *entries() { + for (let row = 0; row < this.rows; row++) { + for (let col = 0; col < this.columns; col++) { + yield [row, col, this.get(row, col)]; + } + } + } + + /** + * iterator from left to right, from top to bottom + * yield value + * @returns {Generator} + */ + *values() { + for (let row = 0; row < this.rows; row++) { + for (let col = 0; col < this.columns; col++) { + yield this.get(row, col); + } + } + } +} + +AbstractMatrix.prototype.klass = 'Matrix'; +if (typeof Symbol !== 'undefined') { + AbstractMatrix.prototype[Symbol.for('nodejs.util.inspect.custom')] = + inspectMatrix; +} + +function compareNumbers(a, b) { + return a - b; +} + +function isArrayOfNumbers(array) { + return array.every((element) => { + return typeof element === 'number'; + }); +} + +// Synonyms +AbstractMatrix.random = AbstractMatrix.rand; +AbstractMatrix.randomInt = AbstractMatrix.randInt; +AbstractMatrix.diagonal = AbstractMatrix.diag; +AbstractMatrix.prototype.diagonal = AbstractMatrix.prototype.diag; +AbstractMatrix.identity = AbstractMatrix.eye; +AbstractMatrix.prototype.negate = AbstractMatrix.prototype.neg; +AbstractMatrix.prototype.tensorProduct = + AbstractMatrix.prototype.kroneckerProduct; + +class Matrix$1 extends AbstractMatrix { + /** + * @type {Float64Array[]} + */ + data; + + /** + * Init an empty matrix + * @param {number} nRows + * @param {number} nColumns + */ + #initData(nRows, nColumns) { + this.data = []; + + if (Number.isInteger(nColumns) && nColumns >= 0) { + for (let i = 0; i < nRows; i++) { + this.data.push(new Float64Array(nColumns)); + } + } else { + throw new TypeError('nColumns must be a positive integer'); + } + + this.rows = nRows; + this.columns = nColumns; + } + + constructor(nRows, nColumns) { + super(); + if (Matrix$1.isMatrix(nRows)) { + this.#initData(nRows.rows, nRows.columns); + Matrix$1.copy(nRows, this); + } else if (Number.isInteger(nRows) && nRows >= 0) { + this.#initData(nRows, nColumns); + } else if (isAnyArray.isAnyArray(nRows)) { + // Copy the values from the 2D array + const arrayData = nRows; + nRows = arrayData.length; + nColumns = nRows ? arrayData[0].length : 0; + if (typeof nColumns !== 'number') { + throw new TypeError( + 'Data must be a 2D array with at least one element', + ); + } + this.data = []; + + for (let i = 0; i < nRows; i++) { + if (arrayData[i].length !== nColumns) { + throw new RangeError('Inconsistent array dimensions'); + } + if (!isArrayOfNumbers(arrayData[i])) { + throw new TypeError('Input data contains non-numeric values'); + } + this.data.push(Float64Array.from(arrayData[i])); + } + + this.rows = nRows; + this.columns = nColumns; + } else { + throw new TypeError( + 'First argument must be a positive number or an array', + ); + } + } + + set(rowIndex, columnIndex, value) { + this.data[rowIndex][columnIndex] = value; + return this; + } + + get(rowIndex, columnIndex) { + return this.data[rowIndex][columnIndex]; + } + + removeRow(index) { + checkRowIndex(this, index); + this.data.splice(index, 1); + this.rows -= 1; + return this; + } + + addRow(index, array) { + if (array === undefined) { + array = index; + index = this.rows; + } + checkRowIndex(this, index, true); + array = Float64Array.from(checkRowVector(this, array)); + this.data.splice(index, 0, array); + this.rows += 1; + return this; + } + + removeColumn(index) { + checkColumnIndex(this, index); + for (let i = 0; i < this.rows; i++) { + const newRow = new Float64Array(this.columns - 1); + for (let j = 0; j < index; j++) { + newRow[j] = this.data[i][j]; + } + for (let j = index + 1; j < this.columns; j++) { + newRow[j - 1] = this.data[i][j]; + } + this.data[i] = newRow; + } + this.columns -= 1; + return this; + } + + addColumn(index, array) { + if (typeof array === 'undefined') { + array = index; + index = this.columns; + } + checkColumnIndex(this, index, true); + array = checkColumnVector(this, array); + for (let i = 0; i < this.rows; i++) { + const newRow = new Float64Array(this.columns + 1); + let j = 0; + for (; j < index; j++) { + newRow[j] = this.data[i][j]; + } + newRow[j++] = array[i]; + for (; j < this.columns + 1; j++) { + newRow[j] = this.data[i][j - 1]; + } + this.data[i] = newRow; + } + this.columns += 1; + return this; + } +} + +installMathOperations(AbstractMatrix, Matrix$1); + +/** + * @typedef {0 | 1 | number | boolean} Mask + */ + +class SymmetricMatrix extends AbstractMatrix { + /** @type {Matrix} */ + #matrix; + + get size() { + return this.#matrix.size; + } + + get rows() { + return this.#matrix.rows; + } + + get columns() { + return this.#matrix.columns; + } + + get diagonalSize() { + return this.rows; + } + + /** + * not the same as matrix.isSymmetric() + * Here is to check if it's instanceof SymmetricMatrix without bundling issues + * + * @param value + * @returns {boolean} + */ + static isSymmetricMatrix(value) { + return Matrix$1.isMatrix(value) && value.klassType === 'SymmetricMatrix'; + } + + /** + * @param diagonalSize + * @return {SymmetricMatrix} + */ + static zeros(diagonalSize) { + return new this(diagonalSize); + } + + /** + * @param diagonalSize + * @return {SymmetricMatrix} + */ + static ones(diagonalSize) { + return new this(diagonalSize).fill(1); + } + + /** + * @param {number | AbstractMatrix | ArrayLike>} diagonalSize + * @return {this} + */ + constructor(diagonalSize) { + super(); + + if (Matrix$1.isMatrix(diagonalSize)) { + if (!diagonalSize.isSymmetric()) { + throw new TypeError('not symmetric data'); + } + + this.#matrix = Matrix$1.copy( + diagonalSize, + new Matrix$1(diagonalSize.rows, diagonalSize.rows), + ); + } else if (Number.isInteger(diagonalSize) && diagonalSize >= 0) { + this.#matrix = new Matrix$1(diagonalSize, diagonalSize); + } else { + this.#matrix = new Matrix$1(diagonalSize); + + if (!this.isSymmetric()) { + throw new TypeError('not symmetric data'); + } + } + } + + clone() { + const matrix = new SymmetricMatrix(this.diagonalSize); + + for (const [row, col, value] of this.upperRightEntries()) { + matrix.set(row, col, value); + } + + return matrix; + } + + toMatrix() { + return new Matrix$1(this); + } + + get(rowIndex, columnIndex) { + return this.#matrix.get(rowIndex, columnIndex); + } + set(rowIndex, columnIndex, value) { + // symmetric set + this.#matrix.set(rowIndex, columnIndex, value); + this.#matrix.set(columnIndex, rowIndex, value); + + return this; + } + + removeCross(index) { + // symmetric remove side + this.#matrix.removeRow(index); + this.#matrix.removeColumn(index); + + return this; + } + + addCross(index, array) { + if (array === undefined) { + array = index; + index = this.diagonalSize; + } + + const row = array.slice(); + row.splice(index, 1); + + this.#matrix.addRow(index, row); + this.#matrix.addColumn(index, array); + + return this; + } + + /** + * @param {Mask[]} mask + */ + applyMask(mask) { + if (mask.length !== this.diagonalSize) { + throw new RangeError('Mask size do not match with matrix size'); + } + + // prepare sides to remove from matrix from mask + /** @type {number[]} */ + const sidesToRemove = []; + for (const [index, passthroughs] of mask.entries()) { + if (passthroughs) continue; + sidesToRemove.push(index); + } + // to remove from highest to lowest for no mutation shifting + sidesToRemove.reverse(); + + // remove sides + for (const sideIndex of sidesToRemove) { + this.removeCross(sideIndex); + } + + return this; + } + + /** + * Compact format upper-right corner of matrix + * iterate from left to right, from top to bottom. + * + * ``` + * A B C D + * A 1 2 3 4 + * B 2 5 6 7 + * C 3 6 8 9 + * D 4 7 9 10 + * ``` + * + * will return compact 1D array `[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]` + * + * length is S(i=0, n=sideSize) => 10 for a 4 sideSized matrix + * + * @returns {number[]} + */ + toCompact() { + const { diagonalSize } = this; + + /** @type {number[]} */ + const compact = new Array((diagonalSize * (diagonalSize + 1)) / 2); + for (let col = 0, row = 0, index = 0; index < compact.length; index++) { + compact[index] = this.get(row, col); + + if (++col >= diagonalSize) col = ++row; + } + + return compact; + } + + /** + * @param {number[]} compact + * @return {SymmetricMatrix} + */ + static fromCompact(compact) { + const compactSize = compact.length; + // compactSize = (sideSize * (sideSize + 1)) / 2 + // https://mathsolver.microsoft.com/fr/solve-problem/y%20%3D%20%20x%20%60cdot%20%20%20%60frac%7B%20%20%60left(%20x%2B1%20%20%60right)%20%20%20%20%7D%7B%202%20%20%7D + // sideSize = (Sqrt(8 × compactSize + 1) - 1) / 2 + const diagonalSize = (Math.sqrt(8 * compactSize + 1) - 1) / 2; + + if (!Number.isInteger(diagonalSize)) { + throw new TypeError( + `This array is not a compact representation of a Symmetric Matrix, ${JSON.stringify( + compact, + )}`, + ); + } + + const matrix = new SymmetricMatrix(diagonalSize); + for (let col = 0, row = 0, index = 0; index < compactSize; index++) { + matrix.set(col, row, compact[index]); + if (++col >= diagonalSize) col = ++row; + } + + return matrix; + } + + /** + * half iterator upper-right-corner from left to right, from top to bottom + * yield [row, column, value] + * + * @returns {Generator<[number, number, number], void, void>} + */ + *upperRightEntries() { + for (let row = 0, col = 0; row < this.diagonalSize; void 0) { + const value = this.get(row, col); + + yield [row, col, value]; + + // at the end of row, move cursor to next row at diagonal position + if (++col >= this.diagonalSize) col = ++row; + } + } + + /** + * half iterator upper-right-corner from left to right, from top to bottom + * yield value + * + * @returns {Generator<[number, number, number], void, void>} + */ + *upperRightValues() { + for (let row = 0, col = 0; row < this.diagonalSize; void 0) { + const value = this.get(row, col); + + yield value; + + // at the end of row, move cursor to next row at diagonal position + if (++col >= this.diagonalSize) col = ++row; + } + } +} +SymmetricMatrix.prototype.klassType = 'SymmetricMatrix'; + +class DistanceMatrix extends SymmetricMatrix { + /** + * not the same as matrix.isSymmetric() + * Here is to check if it's instanceof SymmetricMatrix without bundling issues + * + * @param value + * @returns {boolean} + */ + static isDistanceMatrix(value) { + return ( + SymmetricMatrix.isSymmetricMatrix(value) && + value.klassSubType === 'DistanceMatrix' + ); + } + + constructor(sideSize) { + super(sideSize); + + if (!this.isDistance()) { + throw new TypeError('Provided arguments do no produce a distance matrix'); + } + } + + set(rowIndex, columnIndex, value) { + // distance matrix diagonal is 0 + if (rowIndex === columnIndex) value = 0; + + return super.set(rowIndex, columnIndex, value); + } + + addCross(index, array) { + if (array === undefined) { + array = index; + index = this.diagonalSize; + } + + // ensure distance + array = array.slice(); + array[index] = 0; + + return super.addCross(index, array); + } + + toSymmetricMatrix() { + return new SymmetricMatrix(this); + } + + clone() { + const matrix = new DistanceMatrix(this.diagonalSize); + + for (const [row, col, value] of this.upperRightEntries()) { + if (row === col) continue; + matrix.set(row, col, value); + } + + return matrix; + } + + /** + * Compact format upper-right corner of matrix + * no diagonal (only zeros) + * iterable from left to right, from top to bottom. + * + * ``` + * A B C D + * A 0 1 2 3 + * B 1 0 4 5 + * C 2 4 0 6 + * D 3 5 6 0 + * ``` + * + * will return compact 1D array `[1, 2, 3, 4, 5, 6]` + * + * length is S(i=0, n=sideSize-1) => 6 for a 4 side sized matrix + * + * @returns {number[]} + */ + toCompact() { + const { diagonalSize } = this; + const compactLength = ((diagonalSize - 1) * diagonalSize) / 2; + + /** @type {number[]} */ + const compact = new Array(compactLength); + for (let col = 1, row = 0, index = 0; index < compact.length; index++) { + compact[index] = this.get(row, col); + + if (++col >= diagonalSize) col = ++row + 1; + } + + return compact; + } + + /** + * @param {number[]} compact + */ + static fromCompact(compact) { + const compactSize = compact.length; + + if (compactSize === 0) { + return new this(0); + } + + // compactSize in Natural integer range ]0;∞] + // compactSize = (sideSize * (sideSize - 1)) / 2 + // sideSize = (Sqrt(8 × compactSize + 1) + 1) / 2 + const diagonalSize = (Math.sqrt(8 * compactSize + 1) + 1) / 2; + + if (!Number.isInteger(diagonalSize)) { + throw new TypeError( + `This array is not a compact representation of a DistanceMatrix, ${JSON.stringify( + compact, + )}`, + ); + } + + const matrix = new this(diagonalSize); + for (let col = 1, row = 0, index = 0; index < compactSize; index++) { + matrix.set(col, row, compact[index]); + if (++col >= diagonalSize) col = ++row + 1; + } + + return matrix; + } +} +DistanceMatrix.prototype.klassSubType = 'DistanceMatrix'; + +class BaseView extends AbstractMatrix { + constructor(matrix, rows, columns) { + super(); + this.matrix = matrix; + this.rows = rows; + this.columns = columns; + } +} + +class MatrixColumnView extends BaseView { + constructor(matrix, column) { + checkColumnIndex(matrix, column); + super(matrix, matrix.rows, 1); + this.column = column; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(rowIndex, this.column, value); + return this; + } + + get(rowIndex) { + return this.matrix.get(rowIndex, this.column); + } +} + +class MatrixColumnSelectionView extends BaseView { + constructor(matrix, columnIndices) { + checkColumnIndices(matrix, columnIndices); + super(matrix, matrix.rows, columnIndices.length); + this.columnIndices = columnIndices; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(rowIndex, this.columnIndices[columnIndex], value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(rowIndex, this.columnIndices[columnIndex]); + } +} + +class MatrixFlipColumnView extends BaseView { + constructor(matrix) { + super(matrix, matrix.rows, matrix.columns); + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(rowIndex, this.columns - columnIndex - 1, value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(rowIndex, this.columns - columnIndex - 1); + } +} + +class MatrixFlipRowView extends BaseView { + constructor(matrix) { + super(matrix, matrix.rows, matrix.columns); + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(this.rows - rowIndex - 1, columnIndex, value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(this.rows - rowIndex - 1, columnIndex); + } +} + +class MatrixRowView extends BaseView { + constructor(matrix, row) { + checkRowIndex(matrix, row); + super(matrix, 1, matrix.columns); + this.row = row; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(this.row, columnIndex, value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(this.row, columnIndex); + } +} + +class MatrixRowSelectionView extends BaseView { + constructor(matrix, rowIndices) { + checkRowIndices(matrix, rowIndices); + super(matrix, rowIndices.length, matrix.columns); + this.rowIndices = rowIndices; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(this.rowIndices[rowIndex], columnIndex, value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(this.rowIndices[rowIndex], columnIndex); + } +} + +class MatrixSelectionView extends BaseView { + constructor(matrix, rowIndices, columnIndices) { + checkRowIndices(matrix, rowIndices); + checkColumnIndices(matrix, columnIndices); + super(matrix, rowIndices.length, columnIndices.length); + this.rowIndices = rowIndices; + this.columnIndices = columnIndices; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set( + this.rowIndices[rowIndex], + this.columnIndices[columnIndex], + value, + ); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get( + this.rowIndices[rowIndex], + this.columnIndices[columnIndex], + ); + } +} + +class MatrixSubView extends BaseView { + constructor(matrix, startRow, endRow, startColumn, endColumn) { + checkRange(matrix, startRow, endRow, startColumn, endColumn); + super(matrix, endRow - startRow + 1, endColumn - startColumn + 1); + this.startRow = startRow; + this.startColumn = startColumn; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set( + this.startRow + rowIndex, + this.startColumn + columnIndex, + value, + ); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get( + this.startRow + rowIndex, + this.startColumn + columnIndex, + ); + } +} + +class MatrixTransposeView extends BaseView { + constructor(matrix) { + super(matrix, matrix.columns, matrix.rows); + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(columnIndex, rowIndex, value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(columnIndex, rowIndex); + } +} + +class WrapperMatrix1D extends AbstractMatrix { + constructor(data, options = {}) { + const { rows = 1 } = options; + + if (data.length % rows !== 0) { + throw new Error('the data length is not divisible by the number of rows'); + } + super(); + this.rows = rows; + this.columns = data.length / rows; + this.data = data; + } + + set(rowIndex, columnIndex, value) { + let index = this._calculateIndex(rowIndex, columnIndex); + this.data[index] = value; + return this; + } + + get(rowIndex, columnIndex) { + let index = this._calculateIndex(rowIndex, columnIndex); + return this.data[index]; + } + + _calculateIndex(row, column) { + return row * this.columns + column; + } +} + +class WrapperMatrix2D extends AbstractMatrix { + constructor(data) { + super(); + this.data = data; + this.rows = data.length; + this.columns = data[0].length; + } + + set(rowIndex, columnIndex, value) { + this.data[rowIndex][columnIndex] = value; + return this; + } + + get(rowIndex, columnIndex) { + return this.data[rowIndex][columnIndex]; + } +} + +function wrap(array, options) { + if (isAnyArray.isAnyArray(array)) { + if (array[0] && isAnyArray.isAnyArray(array[0])) { + return new WrapperMatrix2D(array); + } else { + return new WrapperMatrix1D(array, options); + } + } else { + throw new Error('the argument is not an array'); + } +} + +class LuDecomposition { + constructor(matrix) { + matrix = WrapperMatrix2D.checkMatrix(matrix); + + let lu = matrix.clone(); + let rows = lu.rows; + let columns = lu.columns; + let pivotVector = new Float64Array(rows); + let pivotSign = 1; + let i, j, k, p, s, t, v; + let LUcolj, kmax; + + for (i = 0; i < rows; i++) { + pivotVector[i] = i; + } + + LUcolj = new Float64Array(rows); + + for (j = 0; j < columns; j++) { + for (i = 0; i < rows; i++) { + LUcolj[i] = lu.get(i, j); + } + + for (i = 0; i < rows; i++) { + kmax = Math.min(i, j); + s = 0; + for (k = 0; k < kmax; k++) { + s += lu.get(i, k) * LUcolj[k]; + } + LUcolj[i] -= s; + lu.set(i, j, LUcolj[i]); + } + + p = j; + for (i = j + 1; i < rows; i++) { + if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) { + p = i; + } + } + + if (p !== j) { + for (k = 0; k < columns; k++) { + t = lu.get(p, k); + lu.set(p, k, lu.get(j, k)); + lu.set(j, k, t); + } + + v = pivotVector[p]; + pivotVector[p] = pivotVector[j]; + pivotVector[j] = v; + + pivotSign = -pivotSign; + } + + if (j < rows && lu.get(j, j) !== 0) { + for (i = j + 1; i < rows; i++) { + lu.set(i, j, lu.get(i, j) / lu.get(j, j)); + } + } + } + + this.LU = lu; + this.pivotVector = pivotVector; + this.pivotSign = pivotSign; + } + + isSingular() { + let data = this.LU; + let col = data.columns; + for (let j = 0; j < col; j++) { + if (data.get(j, j) === 0) { + return true; + } + } + return false; + } + + solve(value) { + value = Matrix$1.checkMatrix(value); + + let lu = this.LU; + let rows = lu.rows; + + if (rows !== value.rows) { + throw new Error('Invalid matrix dimensions'); + } + if (this.isSingular()) { + throw new Error('LU matrix is singular'); + } + + let count = value.columns; + let X = value.subMatrixRow(this.pivotVector, 0, count - 1); + let columns = lu.columns; + let i, j, k; + + for (k = 0; k < columns; k++) { + for (i = k + 1; i < columns; i++) { + for (j = 0; j < count; j++) { + X.set(i, j, X.get(i, j) - X.get(k, j) * lu.get(i, k)); + } + } + } + for (k = columns - 1; k >= 0; k--) { + for (j = 0; j < count; j++) { + X.set(k, j, X.get(k, j) / lu.get(k, k)); + } + for (i = 0; i < k; i++) { + for (j = 0; j < count; j++) { + X.set(i, j, X.get(i, j) - X.get(k, j) * lu.get(i, k)); + } + } + } + return X; + } + + get determinant() { + let data = this.LU; + if (!data.isSquare()) { + throw new Error('Matrix must be square'); + } + let determinant = this.pivotSign; + let col = data.columns; + for (let j = 0; j < col; j++) { + determinant *= data.get(j, j); + } + return determinant; + } + + get lowerTriangularMatrix() { + let data = this.LU; + let rows = data.rows; + let columns = data.columns; + let X = new Matrix$1(rows, columns); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + if (i > j) { + X.set(i, j, data.get(i, j)); + } else if (i === j) { + X.set(i, j, 1); + } else { + X.set(i, j, 0); + } + } + } + return X; + } + + get upperTriangularMatrix() { + let data = this.LU; + let rows = data.rows; + let columns = data.columns; + let X = new Matrix$1(rows, columns); + for (let i = 0; i < rows; i++) { + for (let j = 0; j < columns; j++) { + if (i <= j) { + X.set(i, j, data.get(i, j)); + } else { + X.set(i, j, 0); + } + } + } + return X; + } + + get pivotPermutationVector() { + return Array.from(this.pivotVector); + } +} + +function hypotenuse(a, b) { + let r = 0; + if (Math.abs(a) > Math.abs(b)) { + r = b / a; + return Math.abs(a) * Math.sqrt(1 + r * r); + } + if (b !== 0) { + r = a / b; + return Math.abs(b) * Math.sqrt(1 + r * r); + } + return 0; +} + +class QrDecomposition { + constructor(value) { + value = WrapperMatrix2D.checkMatrix(value); + + let qr = value.clone(); + let m = value.rows; + let n = value.columns; + let rdiag = new Float64Array(n); + let i, j, k, s; + + for (k = 0; k < n; k++) { + let nrm = 0; + for (i = k; i < m; i++) { + nrm = hypotenuse(nrm, qr.get(i, k)); + } + if (nrm !== 0) { + if (qr.get(k, k) < 0) { + nrm = -nrm; + } + for (i = k; i < m; i++) { + qr.set(i, k, qr.get(i, k) / nrm); + } + qr.set(k, k, qr.get(k, k) + 1); + for (j = k + 1; j < n; j++) { + s = 0; + for (i = k; i < m; i++) { + s += qr.get(i, k) * qr.get(i, j); + } + s = -s / qr.get(k, k); + for (i = k; i < m; i++) { + qr.set(i, j, qr.get(i, j) + s * qr.get(i, k)); + } + } + } + rdiag[k] = -nrm; + } + + this.QR = qr; + this.Rdiag = rdiag; + } + + solve(value) { + value = Matrix$1.checkMatrix(value); + + let qr = this.QR; + let m = qr.rows; + + if (value.rows !== m) { + throw new Error('Matrix row dimensions must agree'); + } + if (!this.isFullRank()) { + throw new Error('Matrix is rank deficient'); + } + + let count = value.columns; + let X = value.clone(); + let n = qr.columns; + let i, j, k, s; + + for (k = 0; k < n; k++) { + for (j = 0; j < count; j++) { + s = 0; + for (i = k; i < m; i++) { + s += qr.get(i, k) * X.get(i, j); + } + s = -s / qr.get(k, k); + for (i = k; i < m; i++) { + X.set(i, j, X.get(i, j) + s * qr.get(i, k)); + } + } + } + for (k = n - 1; k >= 0; k--) { + for (j = 0; j < count; j++) { + X.set(k, j, X.get(k, j) / this.Rdiag[k]); + } + for (i = 0; i < k; i++) { + for (j = 0; j < count; j++) { + X.set(i, j, X.get(i, j) - X.get(k, j) * qr.get(i, k)); + } + } + } + + return X.subMatrix(0, n - 1, 0, count - 1); + } + + isFullRank() { + let columns = this.QR.columns; + for (let i = 0; i < columns; i++) { + if (this.Rdiag[i] === 0) { + return false; + } + } + return true; + } + + get upperTriangularMatrix() { + let qr = this.QR; + let n = qr.columns; + let X = new Matrix$1(n, n); + let i, j; + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + if (i < j) { + X.set(i, j, qr.get(i, j)); + } else if (i === j) { + X.set(i, j, this.Rdiag[i]); + } else { + X.set(i, j, 0); + } + } + } + return X; + } + + get orthogonalMatrix() { + let qr = this.QR; + let rows = qr.rows; + let columns = qr.columns; + let X = new Matrix$1(rows, columns); + let i, j, k, s; + + for (k = columns - 1; k >= 0; k--) { + for (i = 0; i < rows; i++) { + X.set(i, k, 0); + } + X.set(k, k, 1); + for (j = k; j < columns; j++) { + if (qr.get(k, k) !== 0) { + s = 0; + for (i = k; i < rows; i++) { + s += qr.get(i, k) * X.get(i, j); + } + + s = -s / qr.get(k, k); + + for (i = k; i < rows; i++) { + X.set(i, j, X.get(i, j) + s * qr.get(i, k)); + } + } + } + } + return X; + } +} + +class SingularValueDecomposition { + constructor(value, options = {}) { + value = WrapperMatrix2D.checkMatrix(value); + + if (value.isEmpty()) { + throw new Error('Matrix must be non-empty'); + } + + let m = value.rows; + let n = value.columns; + + const { + computeLeftSingularVectors = true, + computeRightSingularVectors = true, + autoTranspose = false, + } = options; + + let wantu = Boolean(computeLeftSingularVectors); + let wantv = Boolean(computeRightSingularVectors); + + let swapped = false; + let a; + if (m < n) { + if (!autoTranspose) { + a = value.clone(); + // eslint-disable-next-line no-console + console.warn( + 'Computing SVD on a matrix with more columns than rows. Consider enabling autoTranspose', + ); + } else { + a = value.transpose(); + m = a.rows; + n = a.columns; + swapped = true; + let aux = wantu; + wantu = wantv; + wantv = aux; + } + } else { + a = value.clone(); + } + + let nu = Math.min(m, n); + let ni = Math.min(m + 1, n); + let s = new Float64Array(ni); + let U = new Matrix$1(m, nu); + let V = new Matrix$1(n, n); + + let e = new Float64Array(n); + let work = new Float64Array(m); + + let si = new Float64Array(ni); + for (let i = 0; i < ni; i++) si[i] = i; + + let nct = Math.min(m - 1, n); + let nrt = Math.max(0, Math.min(n - 2, m)); + let mrc = Math.max(nct, nrt); + + for (let k = 0; k < mrc; k++) { + if (k < nct) { + s[k] = 0; + for (let i = k; i < m; i++) { + s[k] = hypotenuse(s[k], a.get(i, k)); + } + if (s[k] !== 0) { + if (a.get(k, k) < 0) { + s[k] = -s[k]; + } + for (let i = k; i < m; i++) { + a.set(i, k, a.get(i, k) / s[k]); + } + a.set(k, k, a.get(k, k) + 1); + } + s[k] = -s[k]; + } + + for (let j = k + 1; j < n; j++) { + if (k < nct && s[k] !== 0) { + let t = 0; + for (let i = k; i < m; i++) { + t += a.get(i, k) * a.get(i, j); + } + t = -t / a.get(k, k); + for (let i = k; i < m; i++) { + a.set(i, j, a.get(i, j) + t * a.get(i, k)); + } + } + e[j] = a.get(k, j); + } + + if (wantu && k < nct) { + for (let i = k; i < m; i++) { + U.set(i, k, a.get(i, k)); + } + } + + if (k < nrt) { + e[k] = 0; + for (let i = k + 1; i < n; i++) { + e[k] = hypotenuse(e[k], e[i]); + } + if (e[k] !== 0) { + if (e[k + 1] < 0) { + e[k] = 0 - e[k]; + } + for (let i = k + 1; i < n; i++) { + e[i] /= e[k]; + } + e[k + 1] += 1; + } + e[k] = -e[k]; + if (k + 1 < m && e[k] !== 0) { + for (let i = k + 1; i < m; i++) { + work[i] = 0; + } + for (let i = k + 1; i < m; i++) { + for (let j = k + 1; j < n; j++) { + work[i] += e[j] * a.get(i, j); + } + } + for (let j = k + 1; j < n; j++) { + let t = -e[j] / e[k + 1]; + for (let i = k + 1; i < m; i++) { + a.set(i, j, a.get(i, j) + t * work[i]); + } + } + } + if (wantv) { + for (let i = k + 1; i < n; i++) { + V.set(i, k, e[i]); + } + } + } + } + + let p = Math.min(n, m + 1); + if (nct < n) { + s[nct] = a.get(nct, nct); + } + if (m < p) { + s[p - 1] = 0; + } + if (nrt + 1 < p) { + e[nrt] = a.get(nrt, p - 1); + } + e[p - 1] = 0; + + if (wantu) { + for (let j = nct; j < nu; j++) { + for (let i = 0; i < m; i++) { + U.set(i, j, 0); + } + U.set(j, j, 1); + } + for (let k = nct - 1; k >= 0; k--) { + if (s[k] !== 0) { + for (let j = k + 1; j < nu; j++) { + let t = 0; + for (let i = k; i < m; i++) { + t += U.get(i, k) * U.get(i, j); + } + t = -t / U.get(k, k); + for (let i = k; i < m; i++) { + U.set(i, j, U.get(i, j) + t * U.get(i, k)); + } + } + for (let i = k; i < m; i++) { + U.set(i, k, -U.get(i, k)); + } + U.set(k, k, 1 + U.get(k, k)); + for (let i = 0; i < k - 1; i++) { + U.set(i, k, 0); + } + } else { + for (let i = 0; i < m; i++) { + U.set(i, k, 0); + } + U.set(k, k, 1); + } + } + } + + if (wantv) { + for (let k = n - 1; k >= 0; k--) { + if (k < nrt && e[k] !== 0) { + for (let j = k + 1; j < n; j++) { + let t = 0; + for (let i = k + 1; i < n; i++) { + t += V.get(i, k) * V.get(i, j); + } + t = -t / V.get(k + 1, k); + for (let i = k + 1; i < n; i++) { + V.set(i, j, V.get(i, j) + t * V.get(i, k)); + } + } + } + for (let i = 0; i < n; i++) { + V.set(i, k, 0); + } + V.set(k, k, 1); + } + } + + let pp = p - 1; + let eps = Number.EPSILON; + while (p > 0) { + let k, kase; + for (k = p - 2; k >= -1; k--) { + if (k === -1) { + break; + } + const alpha = + Number.MIN_VALUE + eps * Math.abs(s[k] + Math.abs(s[k + 1])); + if (Math.abs(e[k]) <= alpha || Number.isNaN(e[k])) { + e[k] = 0; + break; + } + } + if (k === p - 2) { + kase = 4; + } else { + let ks; + for (ks = p - 1; ks >= k; ks--) { + if (ks === k) { + break; + } + let t = + (ks !== p ? Math.abs(e[ks]) : 0) + + (ks !== k + 1 ? Math.abs(e[ks - 1]) : 0); + if (Math.abs(s[ks]) <= eps * t) { + s[ks] = 0; + break; + } + } + if (ks === k) { + kase = 3; + } else if (ks === p - 1) { + kase = 1; + } else { + kase = 2; + k = ks; + } + } + + k++; + + switch (kase) { + case 1: { + let f = e[p - 2]; + e[p - 2] = 0; + for (let j = p - 2; j >= k; j--) { + let t = hypotenuse(s[j], f); + let cs = s[j] / t; + let sn = f / t; + s[j] = t; + if (j !== k) { + f = -sn * e[j - 1]; + e[j - 1] = cs * e[j - 1]; + } + if (wantv) { + for (let i = 0; i < n; i++) { + t = cs * V.get(i, j) + sn * V.get(i, p - 1); + V.set(i, p - 1, -sn * V.get(i, j) + cs * V.get(i, p - 1)); + V.set(i, j, t); + } + } + } + break; + } + case 2: { + let f = e[k - 1]; + e[k - 1] = 0; + for (let j = k; j < p; j++) { + let t = hypotenuse(s[j], f); + let cs = s[j] / t; + let sn = f / t; + s[j] = t; + f = -sn * e[j]; + e[j] = cs * e[j]; + if (wantu) { + for (let i = 0; i < m; i++) { + t = cs * U.get(i, j) + sn * U.get(i, k - 1); + U.set(i, k - 1, -sn * U.get(i, j) + cs * U.get(i, k - 1)); + U.set(i, j, t); + } + } + } + break; + } + case 3: { + const scale = Math.max( + Math.abs(s[p - 1]), + Math.abs(s[p - 2]), + Math.abs(e[p - 2]), + Math.abs(s[k]), + Math.abs(e[k]), + ); + const sp = s[p - 1] / scale; + const spm1 = s[p - 2] / scale; + const epm1 = e[p - 2] / scale; + const sk = s[k] / scale; + const ek = e[k] / scale; + const b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2; + const c = sp * epm1 * (sp * epm1); + let shift = 0; + if (b !== 0 || c !== 0) { + if (b < 0) { + shift = 0 - Math.sqrt(b * b + c); + } else { + shift = Math.sqrt(b * b + c); + } + shift = c / (b + shift); + } + let f = (sk + sp) * (sk - sp) + shift; + let g = sk * ek; + for (let j = k; j < p - 1; j++) { + let t = hypotenuse(f, g); + if (t === 0) t = Number.MIN_VALUE; + let cs = f / t; + let sn = g / t; + if (j !== k) { + e[j - 1] = t; + } + f = cs * s[j] + sn * e[j]; + e[j] = cs * e[j] - sn * s[j]; + g = sn * s[j + 1]; + s[j + 1] = cs * s[j + 1]; + if (wantv) { + for (let i = 0; i < n; i++) { + t = cs * V.get(i, j) + sn * V.get(i, j + 1); + V.set(i, j + 1, -sn * V.get(i, j) + cs * V.get(i, j + 1)); + V.set(i, j, t); + } + } + t = hypotenuse(f, g); + if (t === 0) t = Number.MIN_VALUE; + cs = f / t; + sn = g / t; + s[j] = t; + f = cs * e[j] + sn * s[j + 1]; + s[j + 1] = -sn * e[j] + cs * s[j + 1]; + g = sn * e[j + 1]; + e[j + 1] = cs * e[j + 1]; + if (wantu && j < m - 1) { + for (let i = 0; i < m; i++) { + t = cs * U.get(i, j) + sn * U.get(i, j + 1); + U.set(i, j + 1, -sn * U.get(i, j) + cs * U.get(i, j + 1)); + U.set(i, j, t); + } + } + } + e[p - 2] = f; + break; + } + case 4: { + if (s[k] <= 0) { + s[k] = s[k] < 0 ? -s[k] : 0; + if (wantv) { + for (let i = 0; i <= pp; i++) { + V.set(i, k, -V.get(i, k)); + } + } + } + while (k < pp) { + if (s[k] >= s[k + 1]) { + break; + } + let t = s[k]; + s[k] = s[k + 1]; + s[k + 1] = t; + if (wantv && k < n - 1) { + for (let i = 0; i < n; i++) { + t = V.get(i, k + 1); + V.set(i, k + 1, V.get(i, k)); + V.set(i, k, t); + } + } + if (wantu && k < m - 1) { + for (let i = 0; i < m; i++) { + t = U.get(i, k + 1); + U.set(i, k + 1, U.get(i, k)); + U.set(i, k, t); + } + } + k++; + } + p--; + break; + } + // no default + } + } + + if (swapped) { + let tmp = V; + V = U; + U = tmp; + } + + this.m = m; + this.n = n; + this.s = s; + this.U = U; + this.V = V; + } + + solve(value) { + let Y = value; + let e = this.threshold; + let scols = this.s.length; + let Ls = Matrix$1.zeros(scols, scols); + + for (let i = 0; i < scols; i++) { + if (Math.abs(this.s[i]) <= e) { + Ls.set(i, i, 0); + } else { + Ls.set(i, i, 1 / this.s[i]); + } + } + + let U = this.U; + let V = this.rightSingularVectors; + + let VL = V.mmul(Ls); + let vrows = V.rows; + let urows = U.rows; + let VLU = Matrix$1.zeros(vrows, urows); + + for (let i = 0; i < vrows; i++) { + for (let j = 0; j < urows; j++) { + let sum = 0; + for (let k = 0; k < scols; k++) { + sum += VL.get(i, k) * U.get(j, k); + } + VLU.set(i, j, sum); + } + } + + return VLU.mmul(Y); + } + + solveForDiagonal(value) { + return this.solve(Matrix$1.diag(value)); + } + + inverse() { + let V = this.V; + let e = this.threshold; + let vrows = V.rows; + let vcols = V.columns; + let X = new Matrix$1(vrows, this.s.length); + + for (let i = 0; i < vrows; i++) { + for (let j = 0; j < vcols; j++) { + if (Math.abs(this.s[j]) > e) { + X.set(i, j, V.get(i, j) / this.s[j]); + } + } + } + + let U = this.U; + + let urows = U.rows; + let ucols = U.columns; + let Y = new Matrix$1(vrows, urows); + + for (let i = 0; i < vrows; i++) { + for (let j = 0; j < urows; j++) { + let sum = 0; + for (let k = 0; k < ucols; k++) { + sum += X.get(i, k) * U.get(j, k); + } + Y.set(i, j, sum); + } + } + + return Y; + } + + get condition() { + return this.s[0] / this.s[Math.min(this.m, this.n) - 1]; + } + + get norm2() { + return this.s[0]; + } + + get rank() { + let tol = Math.max(this.m, this.n) * this.s[0] * Number.EPSILON; + let r = 0; + let s = this.s; + for (let i = 0, ii = s.length; i < ii; i++) { + if (s[i] > tol) { + r++; + } + } + return r; + } + + get diagonal() { + return Array.from(this.s); + } + + get threshold() { + return (Number.EPSILON / 2) * Math.max(this.m, this.n) * this.s[0]; + } + + get leftSingularVectors() { + return this.U; + } + + get rightSingularVectors() { + return this.V; + } + + get diagonalMatrix() { + return Matrix$1.diag(this.s); + } +} + +function inverse$1(matrix, useSVD = false) { + matrix = WrapperMatrix2D.checkMatrix(matrix); + if (useSVD) { + return new SingularValueDecomposition(matrix).inverse(); + } else { + return solve(matrix, Matrix$1.eye(matrix.rows)); + } +} + +function solve(leftHandSide, rightHandSide, useSVD = false) { + leftHandSide = WrapperMatrix2D.checkMatrix(leftHandSide); + rightHandSide = WrapperMatrix2D.checkMatrix(rightHandSide); + if (useSVD) { + return new SingularValueDecomposition(leftHandSide).solve(rightHandSide); + } else { + return leftHandSide.isSquare() + ? new LuDecomposition(leftHandSide).solve(rightHandSide) + : new QrDecomposition(leftHandSide).solve(rightHandSide); + } +} + +function determinant(matrix) { + matrix = Matrix$1.checkMatrix(matrix); + if (matrix.isSquare()) { + if (matrix.columns === 0) { + return 1; + } + + let a, b, c, d; + if (matrix.columns === 2) { + // 2 x 2 matrix + a = matrix.get(0, 0); + b = matrix.get(0, 1); + c = matrix.get(1, 0); + d = matrix.get(1, 1); + + return a * d - b * c; + } else if (matrix.columns === 3) { + // 3 x 3 matrix + let subMatrix0, subMatrix1, subMatrix2; + subMatrix0 = new MatrixSelectionView(matrix, [1, 2], [1, 2]); + subMatrix1 = new MatrixSelectionView(matrix, [1, 2], [0, 2]); + subMatrix2 = new MatrixSelectionView(matrix, [1, 2], [0, 1]); + a = matrix.get(0, 0); + b = matrix.get(0, 1); + c = matrix.get(0, 2); + + return ( + a * determinant(subMatrix0) - + b * determinant(subMatrix1) + + c * determinant(subMatrix2) + ); + } else { + // general purpose determinant using the LU decomposition + return new LuDecomposition(matrix).determinant; + } + } else { + throw Error('determinant can only be calculated for a square matrix'); + } +} + +function xrange(n, exception) { + let range = []; + for (let i = 0; i < n; i++) { + if (i !== exception) { + range.push(i); + } + } + return range; +} + +function dependenciesOneRow( + error, + matrix, + index, + thresholdValue = 10e-10, + thresholdError = 10e-10, +) { + if (error > thresholdError) { + return new Array(matrix.rows + 1).fill(0); + } else { + let returnArray = matrix.addRow(index, [0]); + for (let i = 0; i < returnArray.rows; i++) { + if (Math.abs(returnArray.get(i, 0)) < thresholdValue) { + returnArray.set(i, 0, 0); + } + } + return returnArray.to1DArray(); + } +} + +function linearDependencies(matrix, options = {}) { + const { thresholdValue = 10e-10, thresholdError = 10e-10 } = options; + matrix = Matrix$1.checkMatrix(matrix); + + let n = matrix.rows; + let results = new Matrix$1(n, n); + + for (let i = 0; i < n; i++) { + let b = Matrix$1.columnVector(matrix.getRow(i)); + let Abis = matrix.subMatrixRow(xrange(n, i)).transpose(); + let svd = new SingularValueDecomposition(Abis); + let x = svd.solve(b); + let error = Matrix$1.sub(b, Abis.mmul(x)).abs().max(); + results.setRow( + i, + dependenciesOneRow(error, x, i, thresholdValue, thresholdError), + ); + } + return results; +} + +function pseudoInverse(matrix, threshold = Number.EPSILON) { + matrix = Matrix$1.checkMatrix(matrix); + if (matrix.isEmpty()) { + // with a zero dimension, the pseudo-inverse is the transpose, since all 0xn and nx0 matrices are singular + // (0xn)*(nx0)*(0xn) = 0xn + // (nx0)*(0xn)*(nx0) = nx0 + return matrix.transpose(); + } + let svdSolution = new SingularValueDecomposition(matrix, { autoTranspose: true }); + + let U = svdSolution.leftSingularVectors; + let V = svdSolution.rightSingularVectors; + let s = svdSolution.diagonal; + + for (let i = 0; i < s.length; i++) { + if (Math.abs(s[i]) > threshold) { + s[i] = 1.0 / s[i]; + } else { + s[i] = 0.0; + } + } + + return V.mmul(Matrix$1.diag(s).mmul(U.transpose())); +} + +function covariance(xMatrix, yMatrix = xMatrix, options = {}) { + xMatrix = new Matrix$1(xMatrix); + let yIsSame = false; + if ( + typeof yMatrix === 'object' && + !Matrix$1.isMatrix(yMatrix) && + !isAnyArray.isAnyArray(yMatrix) + ) { + options = yMatrix; + yMatrix = xMatrix; + yIsSame = true; + } else { + yMatrix = new Matrix$1(yMatrix); + } + if (xMatrix.rows !== yMatrix.rows) { + throw new TypeError('Both matrices must have the same number of rows'); + } + const { center = true } = options; + if (center) { + xMatrix = xMatrix.center('column'); + if (!yIsSame) { + yMatrix = yMatrix.center('column'); + } + } + const cov = xMatrix.transpose().mmul(yMatrix); + for (let i = 0; i < cov.rows; i++) { + for (let j = 0; j < cov.columns; j++) { + cov.set(i, j, cov.get(i, j) * (1 / (xMatrix.rows - 1))); + } + } + return cov; +} + +function correlation(xMatrix, yMatrix = xMatrix, options = {}) { + xMatrix = new Matrix$1(xMatrix); + let yIsSame = false; + if ( + typeof yMatrix === 'object' && + !Matrix$1.isMatrix(yMatrix) && + !isAnyArray.isAnyArray(yMatrix) + ) { + options = yMatrix; + yMatrix = xMatrix; + yIsSame = true; + } else { + yMatrix = new Matrix$1(yMatrix); + } + if (xMatrix.rows !== yMatrix.rows) { + throw new TypeError('Both matrices must have the same number of rows'); + } + + const { center = true, scale = true } = options; + if (center) { + xMatrix.center('column'); + if (!yIsSame) { + yMatrix.center('column'); + } + } + if (scale) { + xMatrix.scale('column'); + if (!yIsSame) { + yMatrix.scale('column'); + } + } + + const sdx = xMatrix.standardDeviation('column', { unbiased: true }); + const sdy = yIsSame + ? sdx + : yMatrix.standardDeviation('column', { unbiased: true }); + + const corr = xMatrix.transpose().mmul(yMatrix); + for (let i = 0; i < corr.rows; i++) { + for (let j = 0; j < corr.columns; j++) { + corr.set( + i, + j, + corr.get(i, j) * (1 / (sdx[i] * sdy[j])) * (1 / (xMatrix.rows - 1)), + ); + } + } + return corr; +} + +class EigenvalueDecomposition { + constructor(matrix, options = {}) { + const { assumeSymmetric = false } = options; + + matrix = WrapperMatrix2D.checkMatrix(matrix); + if (!matrix.isSquare()) { + throw new Error('Matrix is not a square matrix'); + } + + if (matrix.isEmpty()) { + throw new Error('Matrix must be non-empty'); + } + + let n = matrix.columns; + let V = new Matrix$1(n, n); + let d = new Float64Array(n); + let e = new Float64Array(n); + let value = matrix; + let i, j; + + let isSymmetric = false; + if (assumeSymmetric) { + isSymmetric = true; + } else { + isSymmetric = matrix.isSymmetric(); + } + + if (isSymmetric) { + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + V.set(i, j, value.get(i, j)); + } + } + tred2(n, e, d, V); + tql2(n, e, d, V); + } else { + let H = new Matrix$1(n, n); + let ort = new Float64Array(n); + for (j = 0; j < n; j++) { + for (i = 0; i < n; i++) { + H.set(i, j, value.get(i, j)); + } + } + orthes(n, H, ort, V); + hqr2(n, e, d, V, H); + } + + this.n = n; + this.e = e; + this.d = d; + this.V = V; + } + + get realEigenvalues() { + return Array.from(this.d); + } + + get imaginaryEigenvalues() { + return Array.from(this.e); + } + + get eigenvectorMatrix() { + return this.V; + } + + get diagonalMatrix() { + let n = this.n; + let e = this.e; + let d = this.d; + let X = new Matrix$1(n, n); + let i, j; + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + X.set(i, j, 0); + } + X.set(i, i, d[i]); + if (e[i] > 0) { + X.set(i, i + 1, e[i]); + } else if (e[i] < 0) { + X.set(i, i - 1, e[i]); + } + } + return X; + } +} + +function tred2(n, e, d, V) { + let f, g, h, i, j, k, hh, scale; + + for (j = 0; j < n; j++) { + d[j] = V.get(n - 1, j); + } + + for (i = n - 1; i > 0; i--) { + scale = 0; + h = 0; + for (k = 0; k < i; k++) { + scale = scale + Math.abs(d[k]); + } + + if (scale === 0) { + e[i] = d[i - 1]; + for (j = 0; j < i; j++) { + d[j] = V.get(i - 1, j); + V.set(i, j, 0); + V.set(j, i, 0); + } + } else { + for (k = 0; k < i; k++) { + d[k] /= scale; + h += d[k] * d[k]; + } + + f = d[i - 1]; + g = Math.sqrt(h); + if (f > 0) { + g = -g; + } + + e[i] = scale * g; + h = h - f * g; + d[i - 1] = f - g; + for (j = 0; j < i; j++) { + e[j] = 0; + } + + for (j = 0; j < i; j++) { + f = d[j]; + V.set(j, i, f); + g = e[j] + V.get(j, j) * f; + for (k = j + 1; k <= i - 1; k++) { + g += V.get(k, j) * d[k]; + e[k] += V.get(k, j) * f; + } + e[j] = g; + } + + f = 0; + for (j = 0; j < i; j++) { + e[j] /= h; + f += e[j] * d[j]; + } + + hh = f / (h + h); + for (j = 0; j < i; j++) { + e[j] -= hh * d[j]; + } + + for (j = 0; j < i; j++) { + f = d[j]; + g = e[j]; + for (k = j; k <= i - 1; k++) { + V.set(k, j, V.get(k, j) - (f * e[k] + g * d[k])); + } + d[j] = V.get(i - 1, j); + V.set(i, j, 0); + } + } + d[i] = h; + } + + for (i = 0; i < n - 1; i++) { + V.set(n - 1, i, V.get(i, i)); + V.set(i, i, 1); + h = d[i + 1]; + if (h !== 0) { + for (k = 0; k <= i; k++) { + d[k] = V.get(k, i + 1) / h; + } + + for (j = 0; j <= i; j++) { + g = 0; + for (k = 0; k <= i; k++) { + g += V.get(k, i + 1) * V.get(k, j); + } + for (k = 0; k <= i; k++) { + V.set(k, j, V.get(k, j) - g * d[k]); + } + } + } + + for (k = 0; k <= i; k++) { + V.set(k, i + 1, 0); + } + } + + for (j = 0; j < n; j++) { + d[j] = V.get(n - 1, j); + V.set(n - 1, j, 0); + } + + V.set(n - 1, n - 1, 1); + e[0] = 0; +} + +function tql2(n, e, d, V) { + let g, h, i, j, k, l, m, p, r, dl1, c, c2, c3, el1, s, s2; + + for (i = 1; i < n; i++) { + e[i - 1] = e[i]; + } + + e[n - 1] = 0; + + let f = 0; + let tst1 = 0; + let eps = Number.EPSILON; + + for (l = 0; l < n; l++) { + tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l])); + m = l; + while (m < n) { + if (Math.abs(e[m]) <= eps * tst1) { + break; + } + m++; + } + + if (m > l) { + do { + + g = d[l]; + p = (d[l + 1] - g) / (2 * e[l]); + r = hypotenuse(p, 1); + if (p < 0) { + r = -r; + } + + d[l] = e[l] / (p + r); + d[l + 1] = e[l] * (p + r); + dl1 = d[l + 1]; + h = g - d[l]; + for (i = l + 2; i < n; i++) { + d[i] -= h; + } + + f = f + h; + + p = d[m]; + c = 1; + c2 = c; + c3 = c; + el1 = e[l + 1]; + s = 0; + s2 = 0; + for (i = m - 1; i >= l; i--) { + c3 = c2; + c2 = c; + s2 = s; + g = c * e[i]; + h = c * p; + r = hypotenuse(p, e[i]); + e[i + 1] = s * r; + s = e[i] / r; + c = p / r; + p = c * d[i] - s * g; + d[i + 1] = h + s * (c * g + s * d[i]); + + for (k = 0; k < n; k++) { + h = V.get(k, i + 1); + V.set(k, i + 1, s * V.get(k, i) + c * h); + V.set(k, i, c * V.get(k, i) - s * h); + } + } + + p = (-s * s2 * c3 * el1 * e[l]) / dl1; + e[l] = s * p; + d[l] = c * p; + } while (Math.abs(e[l]) > eps * tst1); + } + d[l] = d[l] + f; + e[l] = 0; + } + + for (i = 0; i < n - 1; i++) { + k = i; + p = d[i]; + for (j = i + 1; j < n; j++) { + if (d[j] < p) { + k = j; + p = d[j]; + } + } + + if (k !== i) { + d[k] = d[i]; + d[i] = p; + for (j = 0; j < n; j++) { + p = V.get(j, i); + V.set(j, i, V.get(j, k)); + V.set(j, k, p); + } + } + } +} + +function orthes(n, H, ort, V) { + let low = 0; + let high = n - 1; + let f, g, h, i, j, m; + let scale; + + for (m = low + 1; m <= high - 1; m++) { + scale = 0; + for (i = m; i <= high; i++) { + scale = scale + Math.abs(H.get(i, m - 1)); + } + + if (scale !== 0) { + h = 0; + for (i = high; i >= m; i--) { + ort[i] = H.get(i, m - 1) / scale; + h += ort[i] * ort[i]; + } + + g = Math.sqrt(h); + if (ort[m] > 0) { + g = -g; + } + + h = h - ort[m] * g; + ort[m] = ort[m] - g; + + for (j = m; j < n; j++) { + f = 0; + for (i = high; i >= m; i--) { + f += ort[i] * H.get(i, j); + } + + f = f / h; + for (i = m; i <= high; i++) { + H.set(i, j, H.get(i, j) - f * ort[i]); + } + } + + for (i = 0; i <= high; i++) { + f = 0; + for (j = high; j >= m; j--) { + f += ort[j] * H.get(i, j); + } + + f = f / h; + for (j = m; j <= high; j++) { + H.set(i, j, H.get(i, j) - f * ort[j]); + } + } + + ort[m] = scale * ort[m]; + H.set(m, m - 1, scale * g); + } + } + + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + V.set(i, j, i === j ? 1 : 0); + } + } + + for (m = high - 1; m >= low + 1; m--) { + if (H.get(m, m - 1) !== 0) { + for (i = m + 1; i <= high; i++) { + ort[i] = H.get(i, m - 1); + } + + for (j = m; j <= high; j++) { + g = 0; + for (i = m; i <= high; i++) { + g += ort[i] * V.get(i, j); + } + + g = g / ort[m] / H.get(m, m - 1); + for (i = m; i <= high; i++) { + V.set(i, j, V.get(i, j) + g * ort[i]); + } + } + } + } +} + +function hqr2(nn, e, d, V, H) { + let n = nn - 1; + let low = 0; + let high = nn - 1; + let eps = Number.EPSILON; + let exshift = 0; + let norm = 0; + let p = 0; + let q = 0; + let r = 0; + let s = 0; + let z = 0; + let iter = 0; + let i, j, k, l, m, t, w, x, y; + let ra, sa, vr, vi; + let notlast, cdivres; + + for (i = 0; i < nn; i++) { + if (i < low || i > high) { + d[i] = H.get(i, i); + e[i] = 0; + } + + for (j = Math.max(i - 1, 0); j < nn; j++) { + norm = norm + Math.abs(H.get(i, j)); + } + } + + while (n >= low) { + l = n; + while (l > low) { + s = Math.abs(H.get(l - 1, l - 1)) + Math.abs(H.get(l, l)); + if (s === 0) { + s = norm; + } + if (Math.abs(H.get(l, l - 1)) < eps * s) { + break; + } + l--; + } + + if (l === n) { + H.set(n, n, H.get(n, n) + exshift); + d[n] = H.get(n, n); + e[n] = 0; + n--; + iter = 0; + } else if (l === n - 1) { + w = H.get(n, n - 1) * H.get(n - 1, n); + p = (H.get(n - 1, n - 1) - H.get(n, n)) / 2; + q = p * p + w; + z = Math.sqrt(Math.abs(q)); + H.set(n, n, H.get(n, n) + exshift); + H.set(n - 1, n - 1, H.get(n - 1, n - 1) + exshift); + x = H.get(n, n); + + if (q >= 0) { + z = p >= 0 ? p + z : p - z; + d[n - 1] = x + z; + d[n] = d[n - 1]; + if (z !== 0) { + d[n] = x - w / z; + } + e[n - 1] = 0; + e[n] = 0; + x = H.get(n, n - 1); + s = Math.abs(x) + Math.abs(z); + p = x / s; + q = z / s; + r = Math.sqrt(p * p + q * q); + p = p / r; + q = q / r; + + for (j = n - 1; j < nn; j++) { + z = H.get(n - 1, j); + H.set(n - 1, j, q * z + p * H.get(n, j)); + H.set(n, j, q * H.get(n, j) - p * z); + } + + for (i = 0; i <= n; i++) { + z = H.get(i, n - 1); + H.set(i, n - 1, q * z + p * H.get(i, n)); + H.set(i, n, q * H.get(i, n) - p * z); + } + + for (i = low; i <= high; i++) { + z = V.get(i, n - 1); + V.set(i, n - 1, q * z + p * V.get(i, n)); + V.set(i, n, q * V.get(i, n) - p * z); + } + } else { + d[n - 1] = x + p; + d[n] = x + p; + e[n - 1] = z; + e[n] = -z; + } + + n = n - 2; + iter = 0; + } else { + x = H.get(n, n); + y = 0; + w = 0; + if (l < n) { + y = H.get(n - 1, n - 1); + w = H.get(n, n - 1) * H.get(n - 1, n); + } + + if (iter === 10) { + exshift += x; + for (i = low; i <= n; i++) { + H.set(i, i, H.get(i, i) - x); + } + s = Math.abs(H.get(n, n - 1)) + Math.abs(H.get(n - 1, n - 2)); + // eslint-disable-next-line no-multi-assign + x = y = 0.75 * s; + w = -0.4375 * s * s; + } + + if (iter === 30) { + s = (y - x) / 2; + s = s * s + w; + if (s > 0) { + s = Math.sqrt(s); + if (y < x) { + s = -s; + } + s = x - w / ((y - x) / 2 + s); + for (i = low; i <= n; i++) { + H.set(i, i, H.get(i, i) - s); + } + exshift += s; + // eslint-disable-next-line no-multi-assign + x = y = w = 0.964; + } + } + + iter = iter + 1; + + m = n - 2; + while (m >= l) { + z = H.get(m, m); + r = x - z; + s = y - z; + p = (r * s - w) / H.get(m + 1, m) + H.get(m, m + 1); + q = H.get(m + 1, m + 1) - z - r - s; + r = H.get(m + 2, m + 1); + s = Math.abs(p) + Math.abs(q) + Math.abs(r); + p = p / s; + q = q / s; + r = r / s; + if (m === l) { + break; + } + if ( + Math.abs(H.get(m, m - 1)) * (Math.abs(q) + Math.abs(r)) < + eps * + (Math.abs(p) * + (Math.abs(H.get(m - 1, m - 1)) + + Math.abs(z) + + Math.abs(H.get(m + 1, m + 1)))) + ) { + break; + } + m--; + } + + for (i = m + 2; i <= n; i++) { + H.set(i, i - 2, 0); + if (i > m + 2) { + H.set(i, i - 3, 0); + } + } + + for (k = m; k <= n - 1; k++) { + notlast = k !== n - 1; + if (k !== m) { + p = H.get(k, k - 1); + q = H.get(k + 1, k - 1); + r = notlast ? H.get(k + 2, k - 1) : 0; + x = Math.abs(p) + Math.abs(q) + Math.abs(r); + if (x !== 0) { + p = p / x; + q = q / x; + r = r / x; + } + } + + if (x === 0) { + break; + } + + s = Math.sqrt(p * p + q * q + r * r); + if (p < 0) { + s = -s; + } + + if (s !== 0) { + if (k !== m) { + H.set(k, k - 1, -s * x); + } else if (l !== m) { + H.set(k, k - 1, -H.get(k, k - 1)); + } + + p = p + s; + x = p / s; + y = q / s; + z = r / s; + q = q / p; + r = r / p; + + for (j = k; j < nn; j++) { + p = H.get(k, j) + q * H.get(k + 1, j); + if (notlast) { + p = p + r * H.get(k + 2, j); + H.set(k + 2, j, H.get(k + 2, j) - p * z); + } + + H.set(k, j, H.get(k, j) - p * x); + H.set(k + 1, j, H.get(k + 1, j) - p * y); + } + + for (i = 0; i <= Math.min(n, k + 3); i++) { + p = x * H.get(i, k) + y * H.get(i, k + 1); + if (notlast) { + p = p + z * H.get(i, k + 2); + H.set(i, k + 2, H.get(i, k + 2) - p * r); + } + + H.set(i, k, H.get(i, k) - p); + H.set(i, k + 1, H.get(i, k + 1) - p * q); + } + + for (i = low; i <= high; i++) { + p = x * V.get(i, k) + y * V.get(i, k + 1); + if (notlast) { + p = p + z * V.get(i, k + 2); + V.set(i, k + 2, V.get(i, k + 2) - p * r); + } + + V.set(i, k, V.get(i, k) - p); + V.set(i, k + 1, V.get(i, k + 1) - p * q); + } + } + } + } + } + + if (norm === 0) { + return; + } + + for (n = nn - 1; n >= 0; n--) { + p = d[n]; + q = e[n]; + + if (q === 0) { + l = n; + H.set(n, n, 1); + for (i = n - 1; i >= 0; i--) { + w = H.get(i, i) - p; + r = 0; + for (j = l; j <= n; j++) { + r = r + H.get(i, j) * H.get(j, n); + } + + if (e[i] < 0) { + z = w; + s = r; + } else { + l = i; + if (e[i] === 0) { + H.set(i, n, w !== 0 ? -r / w : -r / (eps * norm)); + } else { + x = H.get(i, i + 1); + y = H.get(i + 1, i); + q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; + t = (x * s - z * r) / q; + H.set(i, n, t); + H.set( + i + 1, + n, + Math.abs(x) > Math.abs(z) ? (-r - w * t) / x : (-s - y * t) / z, + ); + } + + t = Math.abs(H.get(i, n)); + if (eps * t * t > 1) { + for (j = i; j <= n; j++) { + H.set(j, n, H.get(j, n) / t); + } + } + } + } + } else if (q < 0) { + l = n - 1; + + if (Math.abs(H.get(n, n - 1)) > Math.abs(H.get(n - 1, n))) { + H.set(n - 1, n - 1, q / H.get(n, n - 1)); + H.set(n - 1, n, -(H.get(n, n) - p) / H.get(n, n - 1)); + } else { + cdivres = cdiv(0, -H.get(n - 1, n), H.get(n - 1, n - 1) - p, q); + H.set(n - 1, n - 1, cdivres[0]); + H.set(n - 1, n, cdivres[1]); + } + + H.set(n, n - 1, 0); + H.set(n, n, 1); + for (i = n - 2; i >= 0; i--) { + ra = 0; + sa = 0; + for (j = l; j <= n; j++) { + ra = ra + H.get(i, j) * H.get(j, n - 1); + sa = sa + H.get(i, j) * H.get(j, n); + } + + w = H.get(i, i) - p; + + if (e[i] < 0) { + z = w; + r = ra; + s = sa; + } else { + l = i; + if (e[i] === 0) { + cdivres = cdiv(-ra, -sa, w, q); + H.set(i, n - 1, cdivres[0]); + H.set(i, n, cdivres[1]); + } else { + x = H.get(i, i + 1); + y = H.get(i + 1, i); + vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; + vi = (d[i] - p) * 2 * q; + if (vr === 0 && vi === 0) { + vr = + eps * + norm * + (Math.abs(w) + + Math.abs(q) + + Math.abs(x) + + Math.abs(y) + + Math.abs(z)); + } + cdivres = cdiv( + x * r - z * ra + q * sa, + x * s - z * sa - q * ra, + vr, + vi, + ); + H.set(i, n - 1, cdivres[0]); + H.set(i, n, cdivres[1]); + if (Math.abs(x) > Math.abs(z) + Math.abs(q)) { + H.set( + i + 1, + n - 1, + (-ra - w * H.get(i, n - 1) + q * H.get(i, n)) / x, + ); + H.set( + i + 1, + n, + (-sa - w * H.get(i, n) - q * H.get(i, n - 1)) / x, + ); + } else { + cdivres = cdiv( + -r - y * H.get(i, n - 1), + -s - y * H.get(i, n), + z, + q, + ); + H.set(i + 1, n - 1, cdivres[0]); + H.set(i + 1, n, cdivres[1]); + } + } + + t = Math.max(Math.abs(H.get(i, n - 1)), Math.abs(H.get(i, n))); + if (eps * t * t > 1) { + for (j = i; j <= n; j++) { + H.set(j, n - 1, H.get(j, n - 1) / t); + H.set(j, n, H.get(j, n) / t); + } + } + } + } + } + } + + for (i = 0; i < nn; i++) { + if (i < low || i > high) { + for (j = i; j < nn; j++) { + V.set(i, j, H.get(i, j)); + } + } + } + + for (j = nn - 1; j >= low; j--) { + for (i = low; i <= high; i++) { + z = 0; + for (k = low; k <= Math.min(j, high); k++) { + z = z + V.get(i, k) * H.get(k, j); + } + V.set(i, j, z); + } + } +} + +function cdiv(xr, xi, yr, yi) { + let r, d; + if (Math.abs(yr) > Math.abs(yi)) { + r = yi / yr; + d = yr + r * yi; + return [(xr + r * xi) / d, (xi - r * xr) / d]; + } else { + r = yr / yi; + d = yi + r * yr; + return [(r * xr + xi) / d, (r * xi - xr) / d]; + } +} + +class CholeskyDecomposition { + constructor(value) { + value = WrapperMatrix2D.checkMatrix(value); + if (!value.isSymmetric()) { + throw new Error('Matrix is not symmetric'); + } + + let a = value; + let dimension = a.rows; + let l = new Matrix$1(dimension, dimension); + let positiveDefinite = true; + let i, j, k; + + for (j = 0; j < dimension; j++) { + let d = 0; + for (k = 0; k < j; k++) { + let s = 0; + for (i = 0; i < k; i++) { + s += l.get(k, i) * l.get(j, i); + } + s = (a.get(j, k) - s) / l.get(k, k); + l.set(j, k, s); + d = d + s * s; + } + + d = a.get(j, j) - d; + + positiveDefinite &&= d > 0; + l.set(j, j, Math.sqrt(Math.max(d, 0))); + for (k = j + 1; k < dimension; k++) { + l.set(j, k, 0); + } + } + + this.L = l; + this.positiveDefinite = positiveDefinite; + } + + isPositiveDefinite() { + return this.positiveDefinite; + } + + solve(value) { + value = WrapperMatrix2D.checkMatrix(value); + + let l = this.L; + let dimension = l.rows; + + if (value.rows !== dimension) { + throw new Error('Matrix dimensions do not match'); + } + if (this.isPositiveDefinite() === false) { + throw new Error('Matrix is not positive definite'); + } + + let count = value.columns; + let B = value.clone(); + let i, j, k; + + for (k = 0; k < dimension; k++) { + for (j = 0; j < count; j++) { + for (i = 0; i < k; i++) { + B.set(k, j, B.get(k, j) - B.get(i, j) * l.get(k, i)); + } + B.set(k, j, B.get(k, j) / l.get(k, k)); + } + } + + for (k = dimension - 1; k >= 0; k--) { + for (j = 0; j < count; j++) { + for (i = k + 1; i < dimension; i++) { + B.set(k, j, B.get(k, j) - B.get(i, j) * l.get(i, k)); + } + B.set(k, j, B.get(k, j) / l.get(k, k)); + } + } + + return B; + } + + get lowerTriangularMatrix() { + return this.L; + } +} + +class nipals { + constructor(X, options = {}) { + X = WrapperMatrix2D.checkMatrix(X); + let { Y } = options; + const { + scaleScores = false, + maxIterations = 1000, + terminationCriteria = 1e-10, + } = options; + + let u; + if (Y) { + if (isAnyArray.isAnyArray(Y) && typeof Y[0] === 'number') { + Y = Matrix$1.columnVector(Y); + } else { + Y = WrapperMatrix2D.checkMatrix(Y); + } + if (Y.rows !== X.rows) { + throw new Error('Y should have the same number of rows as X'); + } + u = Y.getColumnVector(0); + } else { + u = X.getColumnVector(0); + } + + let diff = 1; + let t, q, w, tOld; + + for ( + let counter = 0; + counter < maxIterations && diff > terminationCriteria; + counter++ + ) { + w = X.transpose().mmul(u).div(u.transpose().mmul(u).get(0, 0)); + w = w.div(w.norm()); + + t = X.mmul(w).div(w.transpose().mmul(w).get(0, 0)); + + if (counter > 0) { + diff = t.clone().sub(tOld).pow(2).sum(); + } + tOld = t.clone(); + + if (Y) { + q = Y.transpose().mmul(t).div(t.transpose().mmul(t).get(0, 0)); + q = q.div(q.norm()); + + u = Y.mmul(q).div(q.transpose().mmul(q).get(0, 0)); + } else { + u = t; + } + } + + if (Y) { + let p = X.transpose().mmul(t).div(t.transpose().mmul(t).get(0, 0)); + p = p.div(p.norm()); + let xResidual = X.clone().sub(t.clone().mmul(p.transpose())); + let residual = u.transpose().mmul(t).div(t.transpose().mmul(t).get(0, 0)); + let yResidual = Y.clone().sub( + t.clone().mulS(residual.get(0, 0)).mmul(q.transpose()), + ); + + this.t = t; + this.p = p.transpose(); + this.w = w.transpose(); + this.q = q; + this.u = u; + this.s = t.transpose().mmul(t); + this.xResidual = xResidual; + this.yResidual = yResidual; + this.betas = residual; + } else { + this.w = w.transpose(); + this.s = t.transpose().mmul(t).sqrt(); + if (scaleScores) { + this.t = t.clone().div(this.s.get(0, 0)); + } else { + this.t = t; + } + this.xResidual = X.sub(t.mmul(w.transpose())); + } + } +} + +matrix.AbstractMatrix = AbstractMatrix; +matrix.CHO = CholeskyDecomposition; +matrix.CholeskyDecomposition = CholeskyDecomposition; +matrix.DistanceMatrix = DistanceMatrix; +matrix.EVD = EigenvalueDecomposition; +matrix.EigenvalueDecomposition = EigenvalueDecomposition; +matrix.LU = LuDecomposition; +matrix.LuDecomposition = LuDecomposition; +var Matrix_1 = matrix.Matrix = Matrix$1; +matrix.MatrixColumnSelectionView = MatrixColumnSelectionView; +matrix.MatrixColumnView = MatrixColumnView; +matrix.MatrixFlipColumnView = MatrixFlipColumnView; +matrix.MatrixFlipRowView = MatrixFlipRowView; +matrix.MatrixRowSelectionView = MatrixRowSelectionView; +matrix.MatrixRowView = MatrixRowView; +matrix.MatrixSelectionView = MatrixSelectionView; +matrix.MatrixSubView = MatrixSubView; +matrix.MatrixTransposeView = MatrixTransposeView; +matrix.NIPALS = nipals; +matrix.Nipals = nipals; +matrix.QR = QrDecomposition; +matrix.QrDecomposition = QrDecomposition; +matrix.SVD = SingularValueDecomposition; +matrix.SingularValueDecomposition = SingularValueDecomposition; +matrix.SymmetricMatrix = SymmetricMatrix; +matrix.WrapperMatrix1D = WrapperMatrix1D; +matrix.WrapperMatrix2D = WrapperMatrix2D; +matrix.correlation = correlation; +matrix.covariance = covariance; +var _default = matrix.default = Matrix$1; +matrix.determinant = determinant; +var inverse_1 = matrix.inverse = inverse$1; +matrix.linearDependencies = linearDependencies; +matrix.pseudoInverse = pseudoInverse; +matrix.solve = solve; +matrix.wrap = wrap; + +const Matrix = Matrix_1; +_default.Matrix ? _default.Matrix : Matrix_1; +const inverse = inverse_1; + +/** + * Difference of the matrix function over the parameters + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} evaluatedData - Array of previous evaluated function values + * @param {Array} params - Array of previous parameter values + * @param {number|array} gradientDifference - The step size to approximate the jacobian matrix + * @param {boolean} centralDifference - If true the jacobian matrix is approximated by central differences otherwise by forward differences + * @param {function} paramFunction - The parameters and returns a function with the independent variable as a parameter + * @return {Matrix} + */ + +function gradientFunction( + data, + evaluatedData, + params, + gradientDifference, + paramFunction, + centralDifference, +) { + const nbParams = params.length; + const nbPoints = data.x.length; + let ans = Matrix.zeros(nbParams, nbPoints); + + let rowIndex = 0; + for (let param = 0; param < nbParams; param++) { + if (gradientDifference[param] === 0) continue; + let delta = gradientDifference[param]; + let auxParams = params.slice(); + auxParams[param] += delta; + let funcParam = paramFunction(auxParams); + if (!centralDifference) { + for (let point = 0; point < nbPoints; point++) { + ans.set( + rowIndex, + point, + (evaluatedData[point] - funcParam(data.x[point])) / delta, + ); + } + } else { + auxParams = params.slice(); + auxParams[param] -= delta; + delta *= 2; + let funcParam2 = paramFunction(auxParams); + for (let point = 0; point < nbPoints; point++) { + ans.set( + rowIndex, + point, + (funcParam2(data.x[point]) - funcParam(data.x[point])) / delta, + ); + } + } + rowIndex++; + } + + return ans; +} + +/** + * Matrix function over the samples + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} evaluatedData - Array of previous evaluated function values + * @return {Matrix} + */ +function matrixFunction(data, evaluatedData) { + const m = data.x.length; + + let ans = new Matrix(m, 1); + + for (let point = 0; point < m; point++) { + ans.set(point, 0, data.y[point] - evaluatedData[point]); + } + return ans; +} + +/** + * Iteration for Levenberg-Marquardt + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} params - Array of previous parameter values + * @param {number} damping - Levenberg-Marquardt parameter + * @param {number|array} gradientDifference - The step size to approximate the jacobian matrix + * @param {boolean} centralDifference - If true the jacobian matrix is approximated by central differences otherwise by forward differences + * @param {function} parameterizedFunction - The parameters and returns a function with the independent variable as a parameter + * @return {Array} + */ +function step( + data, + params, + damping, + gradientDifference, + parameterizedFunction, + centralDifference, + weights, +) { + let value = damping; + let identity = Matrix.eye(params.length, params.length, value); + + const func = parameterizedFunction(params); + + let evaluatedData = new Float64Array(data.x.length); + for (let i = 0; i < data.x.length; i++) { + evaluatedData[i] = func(data.x[i]); + } + + let gradientFunc = gradientFunction( + data, + evaluatedData, + params, + gradientDifference, + parameterizedFunction, + centralDifference, + ); + let residualError = matrixFunction(data, evaluatedData); + + let inverseMatrix = inverse( + identity.add( + gradientFunc.mmul( + gradientFunc.transpose().scale('row', { scale: weights }), + ), + ), + ); + + let jacobianWeigthResidualError = gradientFunc.mmul( + residualError.scale('row', { scale: weights }), + ); + + let perturbations = inverseMatrix.mmul(jacobianWeigthResidualError); + + return { + perturbations, + jacobianWeigthResidualError, + }; +} + +/** + * Curve fitting algorithm + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {function} parameterizedFunction - The parameters and returns a function with the independent variable as a parameter + * @param {object} [options] - Options object + * @param {number|array} [options.weights = 1] - weighting vector, if the length does not match with the number of data points, the vector is reconstructed with first value. + * @param {number} [options.damping = 1e-2] - Levenberg-Marquardt parameter, small values of the damping parameter λ result in a Gauss-Newton update and large +values of λ result in a gradient descent update + * @param {number} [options.dampingStepDown = 9] - factor to reduce the damping (Levenberg-Marquardt parameter) when there is not an improvement when updating parameters. + * @param {number} [options.dampingStepUp = 11] - factor to increase the damping (Levenberg-Marquardt parameter) when there is an improvement when updating parameters. + * @param {number} [options.improvementThreshold = 1e-3] - the threshold to define an improvement through an update of parameters + * @param {number|array} [options.gradientDifference = 10e-2] - The step size to approximate the jacobian matrix + * @param {boolean} [options.centralDifference = false] - If true the jacobian matrix is approximated by central differences otherwise by forward differences + * @param {Array} [options.minValues] - Minimum allowed values for parameters + * @param {Array} [options.maxValues] - Maximum allowed values for parameters + * @param {Array} [options.initialValues] - Array of initial parameter values + * @param {number} [options.maxIterations = 100] - Maximum of allowed iterations + * @param {number} [options.errorTolerance = 10e-3] - Minimum uncertainty allowed for each point. + * @param {number} [options.timeout] - maximum time running before throw in seconds. + * @return {{parameterValues: Array, parameterError: number, iterations: number}} + */ +function levenbergMarquardt( + data, + parameterizedFunction, + options = {}, +) { + let { + checkTimeout, + minValues, + maxValues, + parameters, + weightSquare, + damping, + dampingStepUp, + dampingStepDown, + maxIterations, + errorTolerance, + centralDifference, + gradientDifference, + improvementThreshold, + } = checkOptions(data, parameterizedFunction, options); + + let error = errorCalculation( + data, + parameters, + parameterizedFunction, + weightSquare, + ); + + let converged = error <= errorTolerance; + + let iteration = 0; + for (; iteration < maxIterations && !converged; iteration++) { + let previousError = error; + + let { perturbations, jacobianWeigthResidualError } = step( + data, + parameters, + damping, + gradientDifference, + parameterizedFunction, + centralDifference, + weightSquare, + ); + + for (let k = 0; k < parameters.length; k++) { + parameters[k] = Math.min( + Math.max(minValues[k], parameters[k] - perturbations.get(k, 0)), + maxValues[k], + ); + } + + error = errorCalculation( + data, + parameters, + parameterizedFunction, + weightSquare, + ); + + if (isNaN(error)) break; + + let improvementMetric = + (previousError - error) / + perturbations + .transpose() + .mmul(perturbations.mulS(damping).add(jacobianWeigthResidualError)) + .get(0, 0); + + if (improvementMetric > improvementThreshold) { + damping = Math.max(damping / dampingStepDown, 1e-7); + } else { + error = previousError; + damping = Math.min(damping * dampingStepUp, 1e7); + } + + if (checkTimeout()) { + throw new Error( + `The execution time is over to ${options.timeout} seconds`, + ); + } + + converged = error <= errorTolerance; + } + + return { + parameterValues: parameters, + parameterError: error, + iterations: iteration, + }; +} + +export { levenbergMarquardt as default };