Skip to content

Commit

Permalink
wip
Browse files Browse the repository at this point in the history
  • Loading branch information
Peva Blanchard committed Jul 20, 2023
1 parent 3736120 commit 3c2b192
Show file tree
Hide file tree
Showing 10 changed files with 644 additions and 91 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -67,10 +67,10 @@ class Assessment(
}

fun inventory(): Inventory {
val impactFactorMatrix = this.observableMatrix.matrix.matDiv(this.controllableMatrix.matrix.negate())
val impactFactorMatrix = this.controllableMatrix.matrix.negate().matDiv(this.observableMatrix.matrix)
?.let { ImpactFactorMatrix(observablePorts, controllablePorts, it) }
?: throw EvaluatorException("The system cannot be solved")
val supplyMatrix = this.observableMatrix.matrix.transpose().matDiv(demandMatrix.matrix.transpose())
val supplyMatrix = demandMatrix.matrix.transpose().matDiv(this.observableMatrix.matrix.transpose())
?.transpose()
?.let { SupplyMatrix(observablePorts, it) }
?: throw EvaluatorException("The system cannot be solved")
Expand Down
68 changes: 9 additions & 59 deletions src/main/kotlin/ch/kleis/lcaplugin/core/math/DualMatrix.kt
Original file line number Diff line number Diff line change
@@ -1,74 +1,19 @@
package ch.kleis.lcaplugin.core.math

import ch.kleis.lcaplugin.core.math.ejml.EJMLMatrix
import ch.kleis.lcaplugin.core.math.ejml.EJMLMatrixFactory
import org.ejml.data.DMatrixSparseCSC
import org.jetbrains.kotlinx.multik.api.mk
import org.jetbrains.kotlinx.multik.api.zeros
import org.jetbrains.kotlinx.multik.ndarray.data.*
import org.jetbrains.kotlinx.multik.ndarray.operations.forEachMultiIndexed
import org.jetbrains.kotlinx.multik.ndarray.operations.minus
import org.jetbrains.kotlinx.multik.ndarray.operations.plus
import org.jetbrains.kotlinx.multik.ndarray.operations.unaryMinus

/*
TODO: TEST ME !!!
*/

private fun d2Ejml(m: D2Array<Double>): EJMLMatrix {
val (rows, cols) = m.shape
val r = EJMLMatrixFactory.INSTANCE.zero(rows, cols)
m.forEachMultiIndexed { index, d ->
val (i, j) = index
r.set(i, j, d)
}
return r
}

private fun d3Ejml(m: D3Array<Double>): EJMLMatrix {
val (rows, cols, nps) = m.shape
val nParams = if (nps == 0) 1 else nps
val r = EJMLMatrixFactory.INSTANCE.zero(rows, cols * nParams)
m.forEachMultiIndexed { index, d ->
val (i, j, k) = index
r.set(i, j * nParams + k, d)
}
return r
}

private fun ejmlD2(m: EJMLMatrix): D2Array<Double> {
val (rows, cols) = Pair(m.rowDim(), m.colDim())
val r = mk.zeros<Double>(rows, cols)
val iterator = m.matrix.getMatrix<DMatrixSparseCSC>().createCoordinateIterator()
while (iterator.hasNext()) {
val v = iterator.next()
r[v.row, v.col] = v.value
}
return r
}

private fun ejmlD3(m: EJMLMatrix, nParams: Int): D3Array<Double> {
val (rows, cols) = Pair(m.rowDim(), m.colDim())
val r = mk.zeros<Double>(rows, cols, nParams)
if (nParams == 0) return r
val iterator = m.matrix.getMatrix<DMatrixSparseCSC>().createCoordinateIterator()
while (iterator.hasNext()) {
val v = iterator.next()
val i = v.row
val j = v.col / nParams
val k = v.col % nParams
r[i, j, k] = v.value
}
return r
}

open class DualMatrix(
data class DualMatrix(
private val zeroth: D2Array<Double>,
private val first: D3Array<Double>,
) {
private val nParams = first.shape[2]


companion object {
fun zeros(rows: Int, cols: Int, nParams: Int): DualMatrix {
return DualMatrix(
Expand All @@ -78,10 +23,10 @@ open class DualMatrix(
}
}

open fun value(row: Int, col: Int): DualNumber {
fun value(row: Int, col: Int): DualNumber {
return DualNumber(
zeroth[row, col],
first[row, col].asDNArray().asD1Array(),
first[row, col].flatten() as D1Array<Double>,
)
}

Expand Down Expand Up @@ -114,6 +59,11 @@ open class DualMatrix(
)
}

/*
X = R \ L
L @ X = R
*/

fun matDiv(other: DualMatrix): DualMatrix? {
val leftEjml0 = d2Ejml(this.zeroth)
val rightEjml0 = d2Ejml(other.zeroth)
Expand Down Expand Up @@ -157,7 +107,7 @@ open class DualMatrix(

fun set(row: Int, col: Int, value: DualNumber) {
zeroth[row, col] = value.zeroth
first[row, col] = value.first
first[row, col] = if (value.first.isEmpty()) mk.zeros(nParams) else value.first
}

fun add(row: Int, col: Int, value: DualNumber) {
Expand Down
18 changes: 4 additions & 14 deletions src/main/kotlin/ch/kleis/lcaplugin/core/math/DualNumber.kt
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
package ch.kleis.lcaplugin.core.math

import ch.kleis.lcaplugin.core.lang.evaluator.EvaluatorException
import org.jetbrains.kotlinx.multik.api.d1array
import org.jetbrains.kotlinx.multik.api.mk
import org.jetbrains.kotlinx.multik.api.zeros
Expand All @@ -17,15 +16,6 @@ data class DualNumber(
return "$zeroth"
}

private fun broadcast(l: D1Array<Double>, r: D1Array<Double>): Pair<D1Array<Double>, D1Array<Double>> {
return when {
l.isEmpty() -> mk.zeros<Double>(r.size) to r
r.isEmpty() -> l to mk.zeros(l.size)
l.size != r.size -> throw EvaluatorException("d1arrays cannot be broadcast")
else -> l to r
}
}

companion object {
fun constant(c: Double): DualNumber {
return constant(c, 0)
Expand Down Expand Up @@ -64,31 +54,31 @@ data class DualNumber(
}

operator fun plus(other: DualNumber): DualNumber {
val (thisFirst, otherFirst) = broadcast(this.first, other.first)
val (thisFirst, otherFirst) = d1MatchDimensions(this.first, other.first)
return DualNumber(
this.zeroth + other.zeroth,
thisFirst + otherFirst,
)
}

operator fun minus(other: DualNumber): DualNumber {
val (thisFirst, otherFirst) = broadcast(this.first, other.first)
val (thisFirst, otherFirst) = d1MatchDimensions(this.first, other.first)
return DualNumber(
this.zeroth - other.zeroth,
thisFirst - otherFirst,
)
}

operator fun times(other: DualNumber): DualNumber {
val (thisFirst, otherFirst) = broadcast(this.first, other.first)
val (thisFirst, otherFirst) = d1MatchDimensions(this.first, other.first)
return DualNumber(
this.zeroth * other.zeroth,
thisFirst * other.zeroth + otherFirst * zeroth,
)
}

operator fun div(other: DualNumber): DualNumber {
val (thisFirst, otherFirst) = broadcast(this.first, other.first)
val (thisFirst, otherFirst) = d1MatchDimensions(this.first, other.first)
return DualNumber(
this.zeroth / other.zeroth,
thisFirst / other.zeroth - otherFirst * this.zeroth / other.zeroth.pow(2)
Expand Down
69 changes: 69 additions & 0 deletions src/main/kotlin/ch/kleis/lcaplugin/core/math/MathUtils.kt
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
package ch.kleis.lcaplugin.core.math

import ch.kleis.lcaplugin.core.lang.evaluator.EvaluatorException
import ch.kleis.lcaplugin.core.math.ejml.EJMLMatrix
import ch.kleis.lcaplugin.core.math.ejml.EJMLMatrixFactory
import org.ejml.data.DMatrixSparseCSC
import org.jetbrains.kotlinx.multik.api.mk
import org.jetbrains.kotlinx.multik.api.zeros
import org.jetbrains.kotlinx.multik.ndarray.data.D1Array
import org.jetbrains.kotlinx.multik.ndarray.data.D2Array
import org.jetbrains.kotlinx.multik.ndarray.data.D3Array
import org.jetbrains.kotlinx.multik.ndarray.data.set
import org.jetbrains.kotlinx.multik.ndarray.operations.forEachMultiIndexed

fun d1MatchDimensions(l: D1Array<Double>, r: D1Array<Double>): Pair<D1Array<Double>, D1Array<Double>> {
return when {
l.isEmpty() -> mk.zeros<Double>(r.size) to r
r.isEmpty() -> l to mk.zeros(l.size)
l.size != r.size -> throw EvaluatorException("d1arrays cannot be broadcast")
else -> l to r
}
}

fun d2Ejml(m: D2Array<Double>): EJMLMatrix {
val (rows, cols) = m.shape
val r = EJMLMatrixFactory.INSTANCE.zero(rows, cols)
m.forEachMultiIndexed { index, d ->
val (i, j) = index
r.set(i, j, d)
}
return r
}

fun d3Ejml(m: D3Array<Double>): EJMLMatrix {
val (rows, cols, nps) = m.shape
val nParams = if (nps == 0) 1 else nps
val r = EJMLMatrixFactory.INSTANCE.zero(rows, cols * nParams)
m.forEachMultiIndexed { index, d ->
val (i, j, k) = index
r.set(i, j * nParams + k, d)
}
return r
}

fun ejmlD2(m: EJMLMatrix): D2Array<Double> {
val (rows, cols) = Pair(m.rowDim(), m.colDim())
val r = mk.zeros<Double>(rows, cols)
val iterator = m.matrix.getMatrix<DMatrixSparseCSC>().createCoordinateIterator()
while (iterator.hasNext()) {
val v = iterator.next()
r[v.row, v.col] = v.value
}
return r
}

fun ejmlD3(m: EJMLMatrix, nParams: Int): D3Array<Double> {
val (rows, cols) = Pair(m.rowDim(), m.colDim())
val r = mk.zeros<Double>(rows, cols, nParams)
if (nParams == 0) return r
val iterator = m.matrix.getMatrix<DMatrixSparseCSC>().createCoordinateIterator()
while (iterator.hasNext()) {
val v = iterator.next()
val i = v.row
val j = v.col / nParams
val k = v.col % nParams
r[i, j, k] = v.value
}
return r
}
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ class EJMLMatrix(internal val matrix: SimpleMatrix) {

fun matDiv(other: EJMLMatrix): EJMLMatrix? {
val solver = EJMLSolver.INSTANCE
return solver.solve(this, other)
return solver.solve(other, this)
}

fun matMul(other: EJMLMatrix): EJMLMatrix {
Expand Down
92 changes: 92 additions & 0 deletions src/test/kotlin/ch/kleis/lcaplugin/core/math/DualMatrixTest.kt
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
package ch.kleis.lcaplugin.core.math

import ch.kleis.lcaplugin.core.math.MatrixFixture.Companion.makeDualMatrix
import org.junit.Test
import kotlin.test.assertEquals


class DualMatrixTest {
private val dx = DualNumber.basis(2, 0)
private val dy = DualNumber.basis(2, 1)
private fun c(d: Double): DualNumber {
return DualNumber.constant(d)
}

@Test
fun value() {
// given
val m = makeDualMatrix(2, 3, 2,
arrayOf(
c(1.0), c(2.0) + dx, c(3.0),
c(4.0), c(5.0), c(6.0),
)
)

// when
val actual = m.value(0, 1, )

// then
val expected = c(2.0) + dx
assertEquals(expected, actual)
}

@Test
fun set() {
// given
val m = makeDualMatrix(2, 3, 2,
arrayOf(
c(1.0), c(2.0), c(3.0),
c(4.0), c(5.0), c(6.0),
)
)

// when
m.set(0, 1, c(3.0) + dx)

// then
val expected = makeDualMatrix(2, 3, 2,
arrayOf(
c(1.0), c(3.0) + dx, c(3.0),
c(4.0), c(5.0), c(6.0),
)
)
assertEquals(expected, m)
}


@Test
fun minus() {
}

@Test
fun plus() {
}

@Test
fun matMul() {
}

@Test
fun matDiv() {
}

@Test
fun negate() {
}

@Test
fun transpose() {
}

@Test
fun rowDim() {
}

@Test
fun colDim() {
}

@Test
fun add() {
}
}
Loading

0 comments on commit 3c2b192

Please sign in to comment.