commonMain.jetbrains.datalore.plot.base.stat.DensityStatUtil.kt Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of lets-plot-common Show documentation
Lets-Plot JVM package without rendering part
There is a newer version: 4.5.3-alpha1
/*
 * Copyright (c) 2019. JetBrains s.r.o.
 * Use of this source code is governed by the MIT license that can be found in the LICENSE file.
 */

package jetbrains.datalore.plot.base.stat

import jetbrains.datalore.base.gcommon.collect.ClosedRange
import jetbrains.datalore.plot.common.data.SeriesUtil
import kotlin.math.*

object DensityStatUtil {

    private const val DEF_STEP_SIZE = 0.5

    private fun stdDev(data: List): Double {
        var sum = 0.0
        var counter = 0.0

        for (i in data) {
            sum += i
        }
        val mean = sum / data.size
        for (i in data) {
            counter += (i - mean).pow(2.0)
        }
        return sqrt(counter / data.size)
    }

    fun bandWidth(bw: DensityStat.BandWidthMethod, valuesX: List): Double {
        val mySize = valuesX.size

        @Suppress("UNCHECKED_CAST")
        val valuesXFinite = valuesX.filter { SeriesUtil.isFinite(it) } as List
        val dataSummary = FiveNumberSummary(valuesXFinite)
        val myIQR = dataSummary.thirdQuartile - dataSummary.firstQuartile
        val myStdD = stdDev(valuesXFinite)

        when (bw) {
            DensityStat.BandWidthMethod.NRD0 -> {
                if (myIQR > 0) {
                    return 0.9 * min(myStdD, myIQR / 1.34) * mySize.toDouble().pow(-0.2)
                }
                if (myStdD > 0) {
                    return 0.9 * myStdD * mySize.toDouble().pow(-0.2)
                }
            }
            DensityStat.BandWidthMethod.NRD -> {
                if (myIQR > 0) {
                    return 1.06 * min(myStdD, myIQR / 1.34) * mySize.toDouble().pow(-0.2)
                }
                if (myStdD > 0) {
                    return 1.06 * myStdD * mySize.toDouble().pow(-0.2)
                }
            }
        }
        return 1.0
    }

    fun kernel(ker: DensityStat.Kernel): (Double) -> Double {
        return when (ker) {
            DensityStat.Kernel.GAUSSIAN -> { value -> 1 / sqrt(2 * PI) * exp(-0.5 * value.pow(2.0)) }
            DensityStat.Kernel.RECTANGULAR -> { value -> if (abs(value) <= 1) 0.5 else 0.0 }
            DensityStat.Kernel.TRIANGULAR -> { value -> if (abs(value) <= 1) 1 - abs(value) else 0.0 }
            DensityStat.Kernel.BIWEIGHT -> { value -> if (abs(value) <= 1) .9375 * (1 - value * value).pow(2.0) else 0.0 }
            DensityStat.Kernel.EPANECHNIKOV -> { value -> if (abs(value) <= 1) .75 * (1 - value * value) else 0.0 }
            DensityStat.Kernel.OPTCOSINE -> { value -> if (abs(value) <= 1) PI / 4 * cos(PI / 2 * value) else 0.0 }
            else //case COSINE
            -> { value -> if (abs(value) <= 1) (cos(PI * value) + 1) / 2 else 0.0 }
        }
    }

    internal fun densityFunctionFullScan(
        xs: List,
        weights: List,
        ker: (Double) -> Double,
        bw: Double,
        ad: Double
    ): (Double) -> Double {
        val h = bw * ad
        return { x ->
            var sum = 0.0
            for (i in xs.indices) {
                sum += ker((x - xs[i]) / h) * weights[i]
            }
            sum / h
        }
    }

    internal fun densityFunctionFast(
        xs: List,  // must be ordered!
        weights: List,
        ker: (Double) -> Double,
        bw: Double,
        ad: Double
    ): (Double) -> Double {
        val h = bw * ad
        val cutoff = h * 5

        return { x ->
            var sum = 0.0
            var from = xs.binarySearch(x - cutoff)
            if (from < 0) {
                from = -from - 1
            }
            var to = xs.binarySearch(x + cutoff)
            if (to < 0) {
                to = -to - 1
            }

            for (i in (from until to)) {
                sum += ker((x - xs[i]) / h) * weights[i]
            }
            sum / h
        }
    }

    fun createStepValues(range: ClosedRange, n: Int): List {
        val x = ArrayList()
        var min = range.lowerEnd
        var max = range.upperEnd
        val step: Double

        if (max == min) {
            max += DEF_STEP_SIZE
            min -= DEF_STEP_SIZE
        }
        step = (max - min) / (n - 1)
        for (i in 0 until n) {
            x.add(min + step * i)
        }
        return x
    }

    fun toKernel(method: String): DensityStat.Kernel {
        return when (method) {
            "gaussian" -> DensityStat.Kernel.GAUSSIAN
            "rectangular", "uniform" -> DensityStat.Kernel.RECTANGULAR
            "triangular" -> DensityStat.Kernel.TRIANGULAR
            "biweight", "quartic" -> DensityStat.Kernel.BIWEIGHT
            "epanechikov", "parabolic" -> DensityStat.Kernel.EPANECHNIKOV
            "optcosine" -> DensityStat.Kernel.OPTCOSINE
            "cosine" -> DensityStat.Kernel.COSINE
            else -> throw IllegalArgumentException(
                "Unsupported kernel method: '$method'.\n" +
                        "Use one of: gaussian, rectangular, triangular, biweight, epanechikov, optcosine, cos."
            )
        }
    }

    fun toBandWidthMethod(bw: String): DensityStat.BandWidthMethod {
        return when (bw) {
            "nrd0" -> DensityStat.BandWidthMethod.NRD0
            "nrd" -> DensityStat.BandWidthMethod.NRD
            else -> throw IllegalArgumentException(
                "Unsupported bandwidth method: '$bw'.\n" +
                        "Use one of: nrd0, nrd."
            )
        }
    }

    fun createRawMatrix(
        values: List,
        list: List,
        ker: (Double) -> Double,
        bw: Double,
        ad: Double,
        weight: List
    ): Array {
        val a = bw * ad
        val n = values.size
        val x = list.size
        val result = Array(x) { DoubleArray(n) }

        for (row in 0 until x) {
            for (col in 0 until n) {
                result[row][col] = ker((list[row] - values[col]!!) / a) * sqrt(weight[col]!!) / a
            }
        }
        return result
    }
}