diff --git a/math/src/main/scala/breeze/stats/regression/LeastSquares.scala b/math/src/main/scala/breeze/stats/regression/LeastSquares.scala index 727264b3c..08cffd4f4 100644 --- a/math/src/main/scala/breeze/stats/regression/LeastSquares.scala +++ b/math/src/main/scala/breeze/stats/regression/LeastSquares.scala @@ -1,16 +1,27 @@ package breeze.stats.regression + import breeze.generic.UFunc import breeze.linalg._ import org.netlib.util.intW import dev.ludovic.netlib.lapack.LAPACK.{getInstance => lapack} + import java.util.Arrays + +trait NumericType[T] + +object NumericType { + implicit object FloatIsNumeric extends NumericType[Float] + implicit object DoubleIsNumeric extends NumericType[Double] +} private object leastSquaresImplementation { - def doLeastSquares( + + + def doLeastSquaresDouble( data: DenseMatrix[Double], outputs: DenseVector[Double], - workArray: Array[Double]): LeastSquaresRegressionResult = { + workArray: Array[Double]): LeastSquaresRegressionResult[Double] = { require(data.rows == outputs.size) require(data.rows > data.cols + 1) require(workArray.length >= 2 * data.rows * data.cols) @@ -37,139 +48,201 @@ private object leastSquaresImplementation { for (i <- 0 until (data.rows - data.cols)) { r2 = r2 + math.pow(outputs.data(data.cols + i), 2) } - LeastSquaresRegressionResult(coefficients, r2) + LeastSquaresRegressionResult[Double](coefficients, r2) + } + + def doLeastSquaresFloat( + data: DenseMatrix[Float], + outputs: DenseVector[Float], + workArray: Array[Float]): LeastSquaresRegressionResult[Float] = { + require(data.rows == outputs.size) + require(data.rows > data.cols + 1) + require(workArray.length >= 2 * data.rows * data.cols) + + val info = new intW(0) + lapack.sgels( + "N", + data.rows, + data.cols, + 1, + data.data, + data.rows, + outputs.data, + data.rows, + workArray, + workArray.length, + info) + if (info.`val` < 0) { + throw new ArithmeticException("Least squares did not converge.") + } + + val coefficients = new DenseVector[Float](Arrays.copyOf(outputs.data, data.cols)) + var r2 = 0.0 + for (i <- 0 until (data.rows - data.cols)) { + r2 = r2 + math.pow(outputs.data(data.cols + i), 2) + } + LeastSquaresRegressionResult[Float](coefficients, r2.toFloat) } } -case class LeastSquaresRegressionResult(coefficients: DenseVector[Double], rSquared: Double) - extends RegressionResult[DenseVector[Double], Double] { - def apply(x: DenseVector[Double]): Double = coefficients.dot(x) - def apply(X: DenseMatrix[Double]): DenseVector[Double] = X * coefficients + + + +case class LeastSquaresRegressionResult[T](coefficients: DenseVector[T], rSquared: T)(implicit ev: NumericType[T]) + extends RegressionResult[DenseVector[T], T] { + + def apply(x: DenseVector[T]): T = ev match { + + case _: NumericType[Float] => + val coeffs = coefficients.asInstanceOf[DenseVector[Float]] + val x_ins = x.asInstanceOf[DenseVector[Float]] + (coeffs .dot(x_ins)).asInstanceOf[T] + case _: NumericType[Double] => + val coeffs = coefficients.asInstanceOf[DenseVector[Double]] + val x_ins = x.asInstanceOf[DenseVector[Double]] + (coeffs .dot(x_ins)).asInstanceOf[T] + case _ => throw new UnsupportedOperationException("Unsupported numeric type. Only Float and Double are supported") + + } + def apply(X: DenseMatrix[T]): DenseVector[T] = ev match { + case _: NumericType[Float] => + val coeffs = coefficients.asInstanceOf[DenseVector[Float]] + val mat = X.asInstanceOf[DenseMatrix[Float]] + (mat * coeffs).asInstanceOf[DenseVector[T]] + + case _: NumericType[Double] => + val coeffs = coefficients.asInstanceOf[DenseVector[Double]] + val mat = X.asInstanceOf[DenseMatrix[Double]] + (mat * coeffs).asInstanceOf[DenseVector[T]] + + case _ => throw new UnsupportedOperationException("Unsupported numeric type. Only Float and Double are supported") + } } + object leastSquares extends UFunc { implicit val matrixVectorWithWorkArrayDouble - : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] = + new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] { def apply( data: DenseMatrix[Double], outputs: DenseVector[Double], - workArray: Array[Double]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, workArray) + workArray: Array[Double]): LeastSquaresRegressionResult[Double] = + leastSquaresImplementation.doLeastSquaresDouble(data.copy, outputs.copy, workArray) } implicit val matrixVectorWithWorkArrayFloat - : Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] = + new Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] { def apply( data: DenseMatrix[Float], outputs: DenseVector[Float], - workArray: Array[Float]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, workArray) + workArray: Array[Float]): LeastSquaresRegressionResult[Float] = + leastSquaresImplementation.doLeastSquaresFloat(data.copy, outputs.copy, workArray) } implicit val matrixVectorSpecifiedWorkDouble - : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] = + new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] { def apply( data: DenseMatrix[Double], outputs: DenseVector[Double], - workSize: Int): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, new Array[Double](workSize)) + workSize: Int): LeastSquaresRegressionResult[Double] = + leastSquaresImplementation.doLeastSquaresDouble(data.copy, outputs.copy, new Array[Double](workSize)) } implicit val matrixVectorSpecifiedWorkFloat - : Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] = + new Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] { def apply( data: DenseMatrix[Float], outputs: DenseVector[Float], - workSize: Int): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, new Array[Float](workSize)) + workSize: Int): LeastSquaresRegressionResult[Float] = + leastSquaresImplementation.doLeastSquaresFloat(data.copy, outputs.copy, new Array[Float](workSize)) } - implicit val matrixVectorDouble: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] = - new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] { + implicit val matrixVectorDouble: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] = + new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] { def apply( data: DenseMatrix[Double], - outputs: DenseVector[Double]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares( + outputs: DenseVector[Double]): LeastSquaresRegressionResult[Double] = + leastSquaresImplementation.doLeastSquaresDouble( data.copy, outputs.copy, new Array[Double](math.max(1, data.rows * data.cols * 2))) } - implicit val matrixVectorFloat: Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult] = - new Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult] { + implicit val matrixVectorFloat: Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] = + new Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] { def apply( data: DenseMatrix[Float], - outputs: DenseVector[Float]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares( - data.copy, - outputs.copy, + outputs: DenseVector[Float]): LeastSquaresRegressionResult[Float] = + leastSquaresImplementation.doLeastSquaresFloat( + data, + outputs, new Array[Float](math.max(1, data.rows * data.cols * 2))) } } object leastSquaresDestructive extends UFunc { implicit val matrixVectorWithWorkArrayDouble - : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] = + new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] { def apply( data: DenseMatrix[Double], outputs: DenseVector[Double], - workArray: Array[Double]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data, outputs, workArray) + workArray: Array[Double]): LeastSquaresRegressionResult[Double] = + leastSquaresImplementation.doLeastSquaresDouble(data, outputs, workArray) } implicit val matrixVectorWithWorkArrayFloat - : Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] = + new Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] { def apply( data: DenseMatrix[Float], outputs: DenseVector[Float], - workArray: Array[Float]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data, outputs, workArray) + workArray: Array[Float]): LeastSquaresRegressionResult[Float] = + leastSquaresImplementation.doLeastSquaresFloat(data, outputs, workArray) } implicit val matrixVectorSpecifiedWorkDouble - : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] = + new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] { def apply( data: DenseMatrix[Double], outputs: DenseVector[Double], - workSize: Int): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data, outputs, new Array[Float](workSize)) + workSize: Int): LeastSquaresRegressionResult[Double] = + leastSquaresImplementation.doLeastSquaresDouble(data, outputs, new Array[Double](workSize)) } implicit val matrixVectorSpecifiedWorkFloat - : Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult] = - new Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult] { + : Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] = + new Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] { def apply( data: DenseMatrix[Float], outputs: DenseVector[Float], - workSize: Int): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares(data, outputs, new Array[Float](workSize)) + workSize: Int): LeastSquaresRegressionResult[Float] = + leastSquaresImplementation.doLeastSquaresFloat(data, outputs, new Array[Float](workSize)) } - implicit val matrixVectorDouble: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] = - new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] { + implicit val matrixVectorDouble: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] = + new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] { def apply( data: DenseMatrix[Double], - outputs: DenseVector[Double]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares( + outputs: DenseVector[Double]): LeastSquaresRegressionResult[Double] = + leastSquaresImplementation.doLeastSquaresDouble( data, outputs, new Array[Double](math.max(1, data.rows * data.cols * 2))) } - implicit val matrixVectorFloat: Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult] = - new Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult] { + implicit val matrixVectorFloat: Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] = + new Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] { def apply( data: DenseMatrix[Float], - outputs: DenseVector[Float]): LeastSquaresRegressionResult = - leastSquaresImplementation.doLeastSquares( + outputs: DenseVector[Float]): LeastSquaresRegressionResult[Float] = + leastSquaresImplementation.doLeastSquaresFloat( data, outputs, new Array[Float](math.max(1, data.rows * data.cols * 2))) diff --git a/math/src/main/scala/breeze/stats/regression/LeastSquaresDinesh.scala b/math/src/main/scala/breeze/stats/regression/LeastSquaresDinesh.scala deleted file mode 100644 index 17d97cec1..000000000 --- a/math/src/main/scala/breeze/stats/regression/LeastSquaresDinesh.scala +++ /dev/null @@ -1,161 +0,0 @@ -package breeze.stats.regression - -import breeze.generic.UFunc -import breeze.linalg._ -import org.netlib.util.intW -import dev.ludovic.netlib.lapack.LAPACK.{getInstance => lapack} - -import java.util.Arrays -import scala.reflect.ClassTag - -private object leastSquaresImplementation { - def doLeastSquares[T: ClassTag]( - data: DenseMatrix[T], - outputs: DenseVector[T], - workArray: Array[T]) - (implicit num: Numeric[T]): LeastSquaresRegressionResult[T] = { - import num._ - require(data.rows == outputs.size) - require(data.rows > data.cols + 1) - require(workArray.length >= 2 * data.rows * data.cols) - - val info = new intW(0) - lapack match { - case lapack_double if num == implicitly[Numeric[Double]] => - lapack_double.dgels( - "N", - data.rows, - data.cols, - 1, - data.data.asInstanceOf[Array[Double]], - data.rows, - outputs.data.asInstanceOf[Array[Double]], - data.rows, - workArray.asInstanceOf[Array[Double]], - workArray.length, - info) - - if (info.`val` < 0) { - throw new ArithmeticException("Least squares did not converge.") - } - - val coefficients = new DenseVector[Double](Arrays.copyOf(outputs.data.asInstanceOf[Array[Double]], data.cols)) - var r2 = 0.toDouble - for (i <- 0 until (data.rows - data.cols)) { - r2 = r2 + math.pow(outputs.data(data.cols + i).toDouble, num.fromInt(2).toDouble) - } - LeastSquaresRegressionResult(coefficients.asInstanceOf[DenseVector[T]], r2.asInstanceOf[T]) - - - case lapack_float if num == implicitly[Numeric[Float]] => - lapack_float.sgels( - "N", - data.rows, - data.cols, - 1, - data.data.asInstanceOf[Array[Float]], - data.rows, - outputs.data.asInstanceOf[Array[Float]], - data.rows, - workArray.asInstanceOf[Array[Float]], - workArray.length, - info) - - if (info.`val` < 0) { - throw new ArithmeticException("Least squares did not converge.") - } - - val coefficients = new DenseVector[Float](Arrays.copyOf(outputs.data.asInstanceOf[Array[Float]], data.cols)) - var r2 = 0.toDouble - for (i <- 0 until (data.rows - data.cols)) { - r2 = r2 + math.pow(outputs.data(data.cols + i).toDouble, num.fromInt(2).toDouble) - } - LeastSquaresRegressionResult(coefficients.asInstanceOf[DenseVector[T]], r2.toFloat.asInstanceOf[T]) - case _ => - throw new UnsupportedOperationException("Unsupported numeric type. Only Float and Double are supported") - } - - } -} - -case class LeastSquaresRegressionResult[T](coefficients: DenseVector[T], rSquared: T) - extends RegressionResult[DenseVector[T], T] { - def apply(x: DenseVector[T]): T = - - ( coefficients.asInstanceOf[DenseVector[Double]] .dot( x.asInstanceOf[DenseVector[Double]] )).asInstanceOf[T] - - def apply(X: DenseMatrix[T]): DenseVector[T] = - ( X .asInstanceOf[DenseMatrix[Double]] * (coefficients.asInstanceOf[DenseVector[Double]])).asInstanceOf[DenseVector[T]] -} - - -object leastSquares extends UFunc { - implicit def matrixVectorWithWorkArray[T: ClassTag]( - implicit num: Numeric[T]): Impl3[DenseMatrix[T], DenseVector[T], Array[T], LeastSquaresRegressionResult[T]] = - new Impl3[DenseMatrix[T], DenseVector[T], Array[T], LeastSquaresRegressionResult[T]] { - def apply( - data: DenseMatrix[T], - outputs: DenseVector[T], - workArray: Array[T]): LeastSquaresRegressionResult[T] = - leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, workArray) - } - - implicit def matrixVectorSpecifiedWork[T: ClassTag]( - implicit num: Numeric[T]): Impl3[DenseMatrix[T], DenseVector[T], Int, LeastSquaresRegressionResult[T]] = - new Impl3[DenseMatrix[T], DenseVector[T], Int, LeastSquaresRegressionResult[T]] { - def apply( - data: DenseMatrix[T], - outputs: DenseVector[T], - workSize: Int): LeastSquaresRegressionResult[T] = - leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, new Array[T](workSize)) - } - - implicit def matrixVector[T: ClassTag]( - implicit num: Numeric[T]): Impl2[DenseMatrix[T], DenseVector[T], LeastSquaresRegressionResult[T]] = - new Impl2[DenseMatrix[T], DenseVector[T], LeastSquaresRegressionResult[T]] { - def apply( - data: DenseMatrix[T], - outputs: DenseVector[T]): LeastSquaresRegressionResult[T] = - leastSquaresImplementation.doLeastSquares( - data.copy, - outputs.copy, - new Array[T](math.max(1, data.rows * data.cols * 2))) - } -} - - -object leastSquaresDestructive extends UFunc { - - implicit def matrixVectorWithWorkArray[T: ClassTag]( - implicit num: Numeric[T]): Impl3[DenseMatrix[T], DenseVector[T], Array[T], LeastSquaresRegressionResult[T]] = - new Impl3[DenseMatrix[T], DenseVector[T], Array[T], LeastSquaresRegressionResult[T]] { - def apply( - data: DenseMatrix[T], - outputs: DenseVector[T], - workArray: Array[T]): LeastSquaresRegressionResult[T] = - leastSquaresImplementation.doLeastSquares(data, outputs, workArray) - } - - implicit def matrixVectorSpecifiedWork[T: ClassTag]( - implicit num: Numeric[T]): Impl3[DenseMatrix[T], DenseVector[T], Int, LeastSquaresRegressionResult[T]] = - new Impl3[DenseMatrix[T], DenseVector[T], Int, LeastSquaresRegressionResult[T]] { - def apply( - data: DenseMatrix[T], - outputs: DenseVector[T], - workSize: Int): LeastSquaresRegressionResult[T] = - leastSquaresImplementation.doLeastSquares(data, outputs, new Array[T](workSize)) - } - - implicit def matrixVector[T: ClassTag]( - implicit num: Numeric[T]): Impl2[DenseMatrix[T], DenseVector[T], LeastSquaresRegressionResult[T]] = - new Impl2[DenseMatrix[T], DenseVector[T], LeastSquaresRegressionResult[T]] { - def apply( - data: DenseMatrix[T], - outputs: DenseVector[T]): LeastSquaresRegressionResult[T] = - leastSquaresImplementation.doLeastSquares( - data, - outputs, - new Array[T](math.max(1, data.rows * data.cols * 2))) - } - -}