Merge branch 're' of https://github.com/Quafadas/vecxt into re

Author: Quafadas

experiments/src/index.scala

@@ -8,4 +8,5 @@ import io.github.quafadas.plots.SetupVegaBrowser.{*, given}
   val calYrIdx = vecxt_re.CalendarYearIndex(2025, idx.year, idx.idx)
   println(calYrIdx)
   calYrIdx.plotIndex(1.0)
-  println("finished")
+  println("finished")
+end plotIndex

vecxt_re/src-jvm/dist/Dist.scala

@@ -1,55 +1,51 @@
 package vecxt_re
 
-trait Density[T] {
+trait Density[T]:
 
-  /** Returns the unnormalized value of the measure*/
+  /** Returns the unnormalized value of the measure */
   def apply(x: T): Double
 
-  /** Returns the log unnormalized value of the measure*/
+  /** Returns the log unnormalized value of the measure */
   def logApply(x: T): Double = math.log(apply(x))
-}
+end Density
 
-/**
- * Represents a continuous Distribution.
- */
-trait ContinuousDistr[T] extends Density[T] with Rand[T] {
+/** Represents a continuous Distribution.
+  */
+trait ContinuousDistr[T] extends Density[T] with Rand[T]:
 
-  /** Returns the probability density function at that point.*/
+  /** Returns the probability density function at that point. */
   def pdf(x: T): Double = math.exp(logPdf(x))
   def logPdf(x: T): Double = unnormalizedLogPdf(x) - logNormalizer
 
-  /** Returns the probability density function up to a constant at that point.*/
+  /** Returns the probability density function up to a constant at that point. */
   def unnormalizedPdf(x: T): Double = math.exp(unnormalizedLogPdf(x))
 
   def unnormalizedLogPdf(x: T): Double
   def logNormalizer: Double
 
   // 1/Z where Z = exp(logNormalizer)
-  lazy val normalizer
-    : Double = math.exp(-logNormalizer)
+  lazy val normalizer: Double = math.exp(-logNormalizer)
 
   def apply(x: T) = unnormalizedPdf(x)
   override def logApply(x: T) = unnormalizedLogPdf(x)
-}
+end ContinuousDistr
 
-trait HasCdf {
+trait HasCdf:
   def probability(x: Double, y: Double): Double // Probability that P(x < X <= y)
   def cdf(x: Double): Double
 
   // experimental plotting support
   def plot(using viz.LowPriorityPlotTarget): viz.VizReturn
   def plotCdf(using viz.LowPriorityPlotTarget): viz.VizReturn
+end HasCdf
 
-}
+trait HasInverseCdf:
+  def inverseCdf(p: Double): Double // Compute the quantile of p
+end HasInverseCdf
 
-trait HasInverseCdf {
-  def inverseCdf(p: Double): Double //Compute the quantile of p
-}
-
-/**
- * Represents a discrete Distribution
- */
-trait DiscreteDistr[T] extends Density[T] with Rand[T] {
+/** Represents a discrete Distribution
+  */
+trait DiscreteDistr[T] extends Density[T] with Rand[T]:
 
   /** Returns the probability of that draw. */
   def probabilityOf(x: T): Double
@@ -61,12 +57,12 @@ trait DiscreteDistr[T] extends Density[T] with Rand[T] {
 
   def apply(x: T) = unnormalizedProbabilityOf(x)
   override def logApply(x: T) = unnormalizedLogProbabilityOf(x)
-}
+end DiscreteDistr
 
-trait HasMean[T] {
+trait HasMean[T]:
   def mean: T
-}
+end HasMean
 
-trait HasVariance[T] {
+trait HasVariance[T]:
   def variance: T
-}
+end HasVariance

vecxt_re/src-jvm/dist/Empirical.scala

@@ -6,14 +6,13 @@ import io.circe.syntax.*
 
 /** Empirical distribution (JVM only).
   *
-  * This is a nonparametric distribution built directly from observed samples.
-  * It supports positive weights $w_i$.
+  * This is a nonparametric distribution built directly from observed samples. It supports positive weights $w_i$.
   *
-  * The distribution is represented as a discrete measure on the (possibly repeated) sample values:
-  * $$\mathbb{P}(X = x) = \sum_{i: x_i = x} \frac{w_i}{\sum_k w_k}.$$
+  * The distribution is represented as a discrete measure on the (possibly repeated) sample values: $$\mathbb{P}(X = x) =
+  * \sum_{i: x_i = x} \frac{w_i}{\sum_k w_k}.$$
   *
-  * Consequently, the CDF is a right-continuous step function
-  * $$F(t) = \mathbb{P}(X \le t) = \sum_{x \le t} \mathbb{P}(X=x).$$
+  * Consequently, the CDF is a right-continuous step function $$F(t) = \mathbb{P}(X \le t) = \sum_{x \le t}
+  * \mathbb{P}(X=x).$$
   *
   * Sampling is performed by inverse-transform sampling on the cumulative weights.
   *
@@ -44,6 +43,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
       j += 1
     end while
     out
+  end pairs
 
   scala.util.Sorting.stableSort(pairs, (a: (Double, Double), b: (Double, Double)) => a._1 < b._1)
 
@@ -79,6 +79,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
       j += 1
     end while
     out
+  end cdfVals
 
   private val meanVal: Double =
     var s = 0.0
@@ -88,6 +89,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
       j += 1
     end while
     s
+  end meanVal
 
   private val varVal: Double =
     var s2 = 0.0
@@ -98,6 +100,7 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
       j += 1
     end while
     s2
+  end varVal
 
   def mean: Double = meanVal
 
@@ -107,11 +110,14 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
   def probabilityOf(x: Double): Double =
     val idx = java.util.Arrays.binarySearch(xs, x)
     if idx >= 0 then probs(idx) else 0.0
+    end if
+  end probabilityOf
 
   /** Draw a sample using inverse CDF sampling over the atomic masses. */
   def draw: Double =
     val u = rng.nextDouble()
     inverseCdf(u)
+  end draw
 
   /** CDF $F(t)=P(X\le t)$ (right-continuous). */
   def cdf(x: Double): Double =
@@ -136,15 +142,18 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
       val ip = java.util.Arrays.binarySearch(cdfVals, p)
       val idx = if ip >= 0 then ip else -ip - 1
       xs(math.min(idx, xs.length - 1))
+    end if
+  end inverseCdf
 
   /** Plot a (weighted) histogram density estimate. */
   def plot(using viz.LowPriorityPlotTarget) =
     val plot = VegaPlot.fromResource("empiricalPdf.vl.json")
-    val data = (0 until n).map { i => (x = values(i), w = weights(i)) }
+    val data = (0 until n).map(i => (x = values(i), w = weights(i)))
     plot.plot(
       _.data.values := data.asJson,
       _ += (title = s"Empirical Distribution (n=$n)").asJson
     )
+  end plot
 
   /** Plot the empirical CDF (step function). */
   def plotCdf(using viz.LowPriorityPlotTarget) =
@@ -160,19 +169,21 @@ case class Empirical(values: IArray[Double], weights: IArray[Double])
         j += 1
       end while
       pts.toVector
+    end points
 
     plot.plot(
       _.data.values := points.asJson,
       _ += (title = s"Empirical CDF (n=$n)").asJson
     )
+  end plotCdf
 
 end Empirical
 
 object Empirical:
   /** Construct an unweighted empirical distribution (all weights equal to $1$).
     *
-    * Note: We intentionally avoid `apply` overloads here because `IArray[Double]` erases to `Array[Double]`
-    * on the JVM, which can create signature collisions with the case class companion methods.
+    * Note: We intentionally avoid `apply` overloads here because `IArray[Double]` erases to `Array[Double]` on the JVM,
+    * which can create signature collisions with the case class companion methods.
     */
   inline def equalWeights(values: Array[Double]): Empirical =
     Empirical(
@@ -183,4 +194,4 @@ object Empirical:
   /** Construct a weighted empirical distribution from arrays. */
   inline def weighted(values: Array[Double], weights: Array[Double]): Empirical =
     Empirical(IArray.from(values), IArray.from(weights))
-
+end Empirical

vecxt_re/src-jvm/dist/NegativeBinomial.scala

@@ -124,9 +124,14 @@ object NegativeBinomial:
     * @param tol
     *   convergence tolerance for parameter 'a'
     * @return
-    *   Named tuple with `dist`: the fitted NegativeBinomial distribution, and `converged`: whether the optimizer converged within maxIter
+    *   Named tuple with `dist`: the fitted NegativeBinomial distribution, and `converged`: whether the optimizer
+    *   converged within maxIter
     */
-  def mle(observations: Array[Int], maxIter: Int = 500, tol: Double = 1e-8): (dist: NegativeBinomial, converged: Boolean) =
+  def mle(
+      observations: Array[Int],
+      maxIter: Int = 500,
+      tol: Double = 1e-8
+  ): (dist: NegativeBinomial, converged: Boolean) =
     require(observations.nonEmpty, "observations must not be empty")
     require(observations.forall(_ >= 0), "all observations must be non-negative")
 
@@ -209,6 +214,7 @@ object NegativeBinomial:
 
         if aNew <= 0 || !aNew.isFinite then a = a / 2.0
         else a = aNew
+        end if
 
         converged = math.abs(step * delta) < tol * math.abs(a)
         iter += 1
@@ -221,25 +227,21 @@ object NegativeBinomial:
 
   /** Maximum likelihood estimation for the volume-adjusted Negative Binomial.
     *
-    * We observe pairs $(n_j, v_j)$ where $n_j$ is the count and $v_j$ is the volume ratio (historical volume / modeled volume).
-    * With parameters $(r, \beta)$ and $p = 1/(1+\beta v_j)$ the likelihood is
-    * $$
-    * L(r,\beta) = \prod_j \frac{\Gamma(r+n_j)}{\Gamma(r)\,\Gamma(n_j+1)} \left(\frac{\beta v_j}{1+\beta v_j}\right)^{n_j} \left(\frac{1}{1+\beta v_j}\right)^r.
-    * $$
-    * The log-likelihood is
-    * $$
-    * \ell(r,\beta) = \sum_j \big[\log\Gamma(r+n_j) - \log\Gamma(r) - \log\Gamma(n_j+1) + n_j(\log(\beta v_j) - \log(1+\beta v_j)) - r\,\log(1+\beta v_j)\big].
-    * $$
-    * Gradient components:
-    * $$\partial_\beta \ell = \sum_j \Big( \frac{n_j}{\beta(1+\beta v_j)} - \frac{r v_j}{1+\beta v_j} \Big),\quad
-    *   \partial_r \ell = \sum_j \big[\psi(r+n_j) - \psi(r) - \log(1+\beta v_j)\big],$$
-    * and Hessian entries:
-    * $$\partial^2_{\beta\beta} \ell = \sum_j \Big( \frac{r v_j}{(1+\beta v_j)^2} - \frac{n_j(1+2\beta v_j)}{\beta^2(1+\beta v_j)^2} \Big),$$
-    * $$\partial^2_{rr} \ell = \sum_j \big[\psi'(r+n_j) - \psi'(r)\big],\quad \partial^2_{\beta r} \ell = -\sum_j \frac{v_j}{1+\beta v_j}.$$
+    * We observe pairs $(n_j, v_j)$ where $n_j$ is the count and $v_j$ is the volume ratio (historical volume / modeled
+    * volume). With parameters $(r, \beta)$ and $p = 1/(1+\beta v_j)$ the likelihood is $$ L(r,\beta) = \prod_j
+    * \frac{\Gamma(r+n_j)}{\Gamma(r)\,\Gamma(n_j+1)} \left(\frac{\beta v_j}{1+\beta v_j}\right)^{n_j}
+    * \left(\frac{1}{1+\beta v_j}\right)^r. $$ The log-likelihood is $$ \ell(r,\beta) = \sum_j \big[\log\Gamma(r+n_j) -
+    * \log\Gamma(r) - \log\Gamma(n_j+1) + n_j(\log(\beta v_j) - \log(1+\beta v_j)) - r\,\log(1+\beta v_j)\big]. $$
+    * Gradient components: $$\partial_\beta \ell = \sum_j \Big( \frac{n_j}{\beta(1+\beta v_j)} - \frac{r v_j}{1+\beta
+    * v_j} \Big),\quad \partial_r \ell = \sum_j \big[\psi(r+n_j) - \psi(r) - \log(1+\beta v_j)\big],$$ and Hessian
+    * entries: $$\partial^2_{\beta\beta} \ell = \sum_j \Big( \frac{r v_j}{(1+\beta v_j)^2} - \frac{n_j(1+2\beta
+    * v_j)}{\beta^2(1+\beta v_j)^2} \Big),$$ $$\partial^2_{rr} \ell = \sum_j \big[\psi'(r+n_j) - \psi'(r)\big],\quad
+    * \partial^2_{\beta r} \ell = -\sum_j \frac{v_j}{1+\beta v_j}.$$
     *
     * Implementation details:
     *   - Initialize from method of moments on rates $n_j / v_j$; if underdispersed, start at a small $\beta$.
-    *   - Newton updates solve the $2\times2$ system from the gradient/Hessian; a tiny ridge is added to keep the Hessian invertible.
+    *   - Newton updates solve the $2\times2$ system from the gradient/Hessian; a tiny ridge is added to keep the
+    *     Hessian invertible.
     *   - Step halving is applied to enforce positivity of $r$ and $\beta$.
     *
     * @param observations
@@ -253,7 +255,12 @@ object NegativeBinomial:
     * @return
     *   tuple of fitted `NegativeBinomial(r, beta)` and a convergence flag
     */
-  def volweightedMle(observations: Array[Int], volumes: Array[Double], maxIter: Int = 500, tol: Double = 1e-8): (dist: NegativeBinomial, converged: Boolean) =
+  def volweightedMle(
+      observations: Array[Int],
+      volumes: Array[Double],
+      maxIter: Int = 500,
+      tol: Double = 1e-8
+  ): (dist: NegativeBinomial, converged: Boolean) =
     require(observations.nonEmpty, "observations must not be empty")
     require(observations.length == volumes.length, "observations and volumes must have the same length")
     require(observations.forall(_ >= 0), "all observations must be non-negative")
@@ -324,8 +331,7 @@ object NegativeBinomial:
       val hrrAdj = hrr + ridge
       val det = hbbAdj * hrrAdj - hbr * hbr
 
-      if det.isNaN || det.isInfinite || math.abs(det) < 1e-18 then
-        iter = maxIter
+      if det.isNaN || det.isInfinite || math.abs(det) < 1e-18 then iter = maxIter
       else
         val deltaBeta = (gBeta * hrrAdj - gR * hbr) / det
         val deltaR = (hbbAdj * gR - hbr * gBeta) / det
@@ -343,17 +349,23 @@ object NegativeBinomial:
         if newBeta > 0 && newR > 0 && newBeta.isFinite && newR.isFinite then
           beta = newBeta
           r = newR
-          converged =
-            math.abs(step * deltaBeta) <= tol * math.abs(beta) &&
-              math.abs(step * deltaR) <= tol * math.abs(r)
+          converged = math.abs(step * deltaBeta) <= tol * math.abs(beta) &&
+            math.abs(step * deltaR) <= tol * math.abs(r)
         else iter = maxIter
+        end if
       end if
 
       iter += 1
     end while
 
     (NegativeBinomial(r, beta), converged)
-  
-  inline def mleVolumeWeighted(observations: Array[Int], volumes: Array[Double], maxIter: Int = 100, tol: Double = 1e-8): (dist: NegativeBinomial, converged: Boolean) = volweightedMle(observations, volumes, maxIter, tol)
+  end volweightedMle
+
+  inline def mleVolumeWeighted(
+      observations: Array[Int],
+      volumes: Array[Double],
+      maxIter: Int = 100,
+      tol: Double = 1e-8
+  ): (dist: NegativeBinomial, converged: Boolean) = volweightedMle(observations, volumes, maxIter, tol)
 
 end NegativeBinomial

vecxt_re/src-jvm/dist/Pareto.scala

@@ -66,9 +66,9 @@ case class Pareto(scale: Double, shape: Double)
   def variance: Double = distribution.getVariance()
 
   private def guessMaxXForPlot = shape match
-      case s if s > 2 => mean + 4 * math.sqrt(variance) // mean and variance are defined
-      case s if s > 1 => mean + 20 * scale // no well defined variance
-      case _          => scale * 10 // no well defined mean
+    case s if s > 2 => mean + 4 * math.sqrt(variance) // mean and variance are defined
+    case s if s > 1 => mean + 20 * scale // no well defined variance
+    case _          => scale * 10 // no well defined mean
 
   def plot(using viz.LowPriorityPlotTarget) =
 
@@ -93,7 +93,6 @@ case class Pareto(scale: Double, shape: Double)
       _.layer._1.data.sequence.stop := maxX,
       _.layer._1.data.sequence.step := (maxX - scale) / 200,
       _.layer._1.transform.head.calculate := pdfExpr,
-      
       _ += (title = s"Pareto Distribution PDF (scale=$scale, shape=$shape)").asJson
     )
   end plot