Commit 829ef182 authored by Eva Darulova's avatar Eva Darulova

Refactor absolute error computation (step 1)

The following has changed:
- absolute error computation previously present in ErrorFunctions
  is now split into range computation and error computation.
- it does not annotate the tree anymore, the information is passed via maps in the context
- error computation supports mixed-precision (without a front-end)

What is probably broken or needs fixing:
- fixed-point code generation
- RangePhase and AbsErrorPhases cleanup and deduplication
parent 2f848c4b
......@@ -28,15 +28,17 @@ case class Context(
// but don't want to pollute the nice and clean trees.
// If these get too many, move to their own "Summary".
// indexed by FunDef.id
inputRanges: Map[Identifier, Map[Identifier, Interval]] = Map(),
inputErrors: Map[Identifier, Map[Identifier, Rational]] = Map(),
specInputRanges: Map[Identifier, Map[Identifier, Interval]] = Map(),
specInputErrors: Map[Identifier, Map[Identifier, Rational]] = Map(),
// for now we only support a single result value, i.e. no tuples
// this map is indexed by fnc.id -> potentially partial interval bound of result
// and similar for the errors
resultRangeBounds: Map[Identifier, PartialInterval] = Map(),
resultErrorBounds: Map[Identifier, Rational] = Map()
//requiredOutputRanges: Map[Identifier, Map[Identifier, PartialInterval]] = Map(),
//requiredOutputErrors: Map[Identifier, Map[Identifier, Rational]] = Map()
specResultRangeBounds: Map[Identifier, PartialInterval] = Map(),
specResultErrorBounds: Map[Identifier, Rational] = Map(),
// the analysed/computed roundoff errors for each function
resultAbsoluteErrors: Map[Identifier, Rational] = Map(),
resultRealRanges: Map[Identifier, Interval] = Map()
) {
// on the first creation of a context, we also update the context variable
......
......@@ -50,23 +50,16 @@ object InfoPhase extends DaisyPhase {
reporter.info(fnc.id)
val infoString: String = getLastExpression(fnc.body.get) match {
case x: NumAnnotation if x.hasError =>
val absError = x.absError
val range = x.interval
val relError: Double = if (range.xlo <= zero && zero <= range.xhi) {
Double.NaN
} else {
max(abs(absError / range.xlo), abs(absError / range.xhi)).toDouble
}
s"abs-error: ${x.absError}, range: ${x.interval},\nrel-error: ${relError}"
case x: NumAnnotation =>
"final error: - "
case _ => ""
}
val absError = ctx.resultAbsoluteErrors(fnc.id)
val range = ctx.resultRealRanges(fnc.id)
val relError: Double = if (range.xlo <= zero && zero <= range.xhi) {
Double.NaN
} else {
max(abs(absError / range.xlo), abs(absError / range.xhi)).toDouble
}
val infoString: String =
s"abs-error: ${absError}, real range: ${range},\nrel-error: ${relError}"
reporter.info(infoString)
if (out != null) {
......
......@@ -106,7 +106,7 @@ object DynamicPhase extends DaisyPhase {
val body = fnc.body.get
reporter.info("evaluating " + id + "...")
val inputRanges: Map[Identifier, Interval] = ctx.inputRanges(id).map({
val inputRanges: Map[Identifier, Interval] = ctx.specInputRanges(id).map({
case (id, i) =>
(id, Interval(i.mid - inputRangeFactor * i.radius,
i.mid + inputRangeFactor * i.radius))
......
......@@ -16,7 +16,7 @@ import lang.Identifiers._
Prerequisites:
- SpecsProcessingPhase
*/
object RangeErrorPhase extends DaisyPhase with ErrorFunctions {
object RangeErrorPhase extends DaisyPhase with RoundoffEvaluators with IntervalSubdivision {
override val name = "range-error phase"
override val description = "Computes ranges and absolute errors"
......@@ -37,10 +37,6 @@ object RangeErrorPhase extends DaisyPhase with ErrorFunctions {
var reporter: Reporter = null
// Var for error functions
// trackRoundoffErrs and uniformPrecision assigned below in run
attachToTree = true
override def run(ctx: Context, prg: Program): (Context, Program) = {
reporter = ctx.reporter
reporter.info(s"\nStarting $name")
......@@ -51,9 +47,8 @@ object RangeErrorPhase extends DaisyPhase with ErrorFunctions {
var errorMethod = "affine" // only one supported so far
var trackInitialErrs = true
var trackRoundoffErrs = true
// setting trait variables
trackRoundoffErrs = true
var uniformPrecision: Precision = Float64
// process relevant options
......@@ -92,11 +87,13 @@ object RangeErrorPhase extends DaisyPhase with ErrorFunctions {
val fncsToConsider: Seq[String] = functionsToConsider(ctx, prg)
for (fnc <- prg.defs)
if (!fnc.precondition.isEmpty && !fnc.body.isEmpty && fncsToConsider.contains(fnc.id.toString)){
val res: Map[Identifier, (Rational, Interval)] = prg.defs.filter(fnc =>
!fnc.precondition.isEmpty &&
!fnc.body.isEmpty &&
fncsToConsider.contains(fnc.id.toString)).map(fnc => {
reporter.info("analyzing fnc: " + fnc.id)
val inputValMap: Map[Identifier, Interval] = ctx.inputRanges(fnc.id)
val inputValMap: Map[Identifier, Interval] = ctx.specInputRanges(fnc.id)
// If we track both input and roundoff errors, then we pre-compute
// the roundoff errors for those variables that do not have a user-defined
......@@ -104,14 +101,14 @@ object RangeErrorPhase extends DaisyPhase with ErrorFunctions {
val inputErrorMap: Map[Identifier, Rational] =
if (trackInitialErrs && trackRoundoffErrs){
val inputErrs = ctx.inputErrors(fnc.id)
val inputErrs = ctx.specInputErrors(fnc.id)
val allIDs = fnc.params.map(_.id).toSet
val missingIDs = allIDs -- inputErrs.keySet
inputErrs ++ missingIDs.map( id => (id -> uniformPrecision.absRoundoff(inputValMap(id))))
} else if(trackInitialErrs) {
val inputErrs = ctx.inputErrors(fnc.id)
val inputErrs = ctx.specInputErrors(fnc.id)
val allIDs = fnc.params.map(_.id).toSet
val missingIDs = allIDs -- inputErrs.keySet
inputErrs ++ missingIDs.map( id => (id -> zero))
......@@ -131,35 +128,38 @@ object RangeErrorPhase extends DaisyPhase with ErrorFunctions {
// TODO: Interval-only based error estimation; should be a very quick fix
(rangeMethod, errorMethod) match {
val (resError: Rational, resRange: Interval) = (rangeMethod, errorMethod) match {
case ("interval", "affine") =>
errorIntervalAffine(fnc.body.get, inputValMap, inputErrorMap, uniformPrecision)
uniformRoundoff_IA_AA(fnc.body.get, inputValMap, inputErrorMap, uniformPrecision, trackRoundoffErrs)
case ("affine", "affine") =>
errorAffineAffine(fnc.body.get, inputValMap, inputErrorMap, uniformPrecision)
uniformRoundoff_AA_AA(fnc.body.get, inputValMap, inputErrorMap, uniformPrecision, trackRoundoffErrs)
case ("smt", "affine") =>
errorSMTAffine(fnc.body.get, inputValMap, inputErrorMap, uniformPrecision)
uniformRoundoff_SMT_AA(fnc.body.get, inputValMap, inputErrorMap, uniformPrecision, trackRoundoffErrs)
// default is to use the method that attaches the info to trees.
case ("subdiv", _) =>
evaluateSubdiv(
fnc.body.get,
val tmp = doIntervalSubdivision( //evaluateSubdiv(
fnc.body.get, lang.TreeOps.freeVariablesOf(fnc.body.get),
inputValMap,
inputErrorMap,
trackRoundoffErrs,
uniformPrecision)
(tmp._2, tmp._1)
case _ =>
reporter.fatalError(s"Your combination of $rangeMethod and $errorMethod" +
"for computing ranges and errors is not supported.")
null
}
}
(fnc.id -> (resError, resRange))
}).toMap
timer.stop
ctx.reporter.info(s"Finished $name")
(ctx, prg)
(ctx.copy(resultAbsoluteErrors = res.map(x => (x._1 -> x._2._1)),
resultRealRanges = res.map(x => (x._1 -> x._2._2))), prg)
}
......
......@@ -70,7 +70,7 @@ object RangePhase extends DaisyPhase {
case "smt" =>
evaluateSMT(fnc.body.get, Map.empty)
case "subdiv" =>
evaluateSubdiv(fnc.body.get, ctx.inputRanges(fnc.id), Map.empty)
evaluateSubdiv(fnc.body.get, ctx.specInputRanges(fnc.id), Map.empty)
}
......
......@@ -121,8 +121,8 @@ object SpecsProcessingPhase extends DaisyPhase {
reporter.debug("range bounds: " + resRanges.mkString("\n"))
reporter.debug("error bounds: " + resErrors.mkString("\n"))
(ctx.copy(inputRanges = allRanges, inputErrors = allErrors,
resultRangeBounds = resRanges, resultErrorBounds = resErrors), prg)
(ctx.copy(specInputRanges = allRanges, specInputErrors = allErrors,
specResultRangeBounds = resRanges, specResultErrorBounds = resErrors), prg)
}
......
......@@ -173,7 +173,7 @@ object TreeOps {
}
/** Returns the set of free variables in an expression */
def variablesOf(expr: Expr): Set[Identifier] = {
def freeVariablesOf(expr: Expr): Set[Identifier] = {
fold[Set[Identifier]] {
case (e, subs) =>
val subvs = subs.flatten.toSet
......@@ -186,6 +186,20 @@ object TreeOps {
}(expr)
}
/** Returns the set of all variables occuring in an expression */
def allVariablesOf(expr: Expr): Set[Identifier] = {
fold[Set[Identifier]] {
case (e, subs) =>
val subvs = subs.flatten.toSet
e match {
case Variable(i) => subvs + i
//case Let(i, _, _) => subvs - i
//case Lambda(args, _) => subvs -- args.map(_.id)
case _ => subvs
}
}(expr)
}
/**
A term is an expression which (for our purposes) does not contain
propositional logic connectives.
......
......@@ -22,7 +22,7 @@ import lang.Types._
import lang.Trees._
import lang.Identifiers._
import utils.Bijection
import lang.TreeOps.variablesOf
import lang.TreeOps.freeVariablesOf
import utils.Rational
import Rational._
import lang.Constructors._
......@@ -102,7 +102,7 @@ abstract class SMTLibSolver(val context: Context) {
/* Add a constraint */
def assertConstraint(expr: Expr): Unit = {
try {
variablesOf(expr).foreach(declareVariable)
freeVariablesOf(expr).foreach(declareVariable)
val term = toSMT(expr)(Map())
emit(SMTAssert(term))
......
package daisy
package utils
import analysis._
import analysis.DynamicPhase._
import lang.Trees._
import lang.Identifiers._
import lang.NumAnnotation
import Interval._
import Rational._
import FinitePrecision._
import daisy.analysis.Sampler._
import MPFRFloat.{abs => mpfr_abs, max => mpfr_max, min => mpfr_min}
import scala.collection.immutable._
trait DynamicEvaluators {
var dynamicSamplesDefault: Int = 100000
/**
* Estimates the (roundoff) error of an expression dynamically
* by executing the finite-precision computation side-by-side with a higher-
* precision one in MPFR.
* This version does not automatically add any roundoff to the inputs,
* you have to specify input errors manually, or map all variable identifiers
* to zero if you don't want any input errors.
* Errors added to each input are randomly taken from the interval
* [- input-error, + input-error]
*
* @param expr expression whose error is to be calculated
* @param inputValMap map from variable identifier to its input interval
* @param inputErrorMap map from variable identifier to its max. input error
* @param dynamicSamples how many samples to use
*/
def errorDynamic(
expr: Expr,
inputValMap: Map[Identifier, Interval],
inputErrorMap: Map[Identifier, Rational],
dynamicSamples: Int = dynamicSamplesDefault): Rational = {
val inputRanges: Map[Identifier, Interval] = inputValMap.map({
case (id, i) =>
(id, Interval(i.mid - inputRangeFactor * i.radius,
i.mid + inputRangeFactor * i.radius))
})
val rand = new util.Random(3691285)
val sampler = new Uniform(inputRanges, 485793) //no seed = System millis
// TODO: this should not be using the ErrorMeasurer as it only needs the abs error
val measurer = new ErrorMeasurerMPFR()
var currentMaxAbsMPFR = measurer.maxAbsError
var currentMaxAbs: Double = measurer.maxAbsError.doubleValue
var i = 0
while (i < dynamicSamples) {
i = i + 1
// no input errors
val dblInputsWithoutErrors: Map[Identifier, Double] = sampler.next
// baseline without errors
val mpfrInputs: Map[Identifier, MPFRFloat] = dblInputsWithoutErrors.map({
case (x, d) => (x -> MPFRFloat.fromDouble(d))
})
// add errors to lower precision version, if there are any
val dblInputs: Map[Identifier, Double] = dblInputsWithoutErrors.map({
case (x, d) =>
if (rand.nextBoolean) (x -> (d + inputErrorMap(x).toDouble * rand.nextDouble))
else (x -> (d - inputErrorMap(x).toDouble * rand.nextDouble))
})
val dblOutput: Double = evalDouble(expr, dblInputs)
val mpfrOutput: MPFRFloat = evalMPFR(expr, mpfrInputs)
measurer.nextValues(dblOutput, mpfrOutput)
// Invariant that absolute errors have to grow monotonically
assert(currentMaxAbsMPFR <= measurer.maxAbsError)
currentMaxAbsMPFR = measurer.maxAbsError
assert(currentMaxAbs <= measurer.maxAbsError.doubleValue)
currentMaxAbs = measurer.maxAbsError.doubleValue
}
Rational.fromString(measurer.maxAbsError.toString())
}
/**
* Estimates the (roundoff) error of an expression dynamically
* by executing the finite-precision computation side-by-side with a higher-
* precision one in MPFR.
* In this version roundoff errors on inputs are considered by converting
* a double-valued sample into a String and considering that string as the
* 'baseline'. If you don't want roundoff errors, use the other errorDynamic
* function.
*
* @param expr expression whose error is to be calculated
* @param inputValMap map from variable identifier to its input interval
* @param dynamicSamples how many samples to use
*/
def errorDynamicWithInputRoundoff(
expr: Expr,
inputValMap: Map[Identifier, Interval],
dynamicSamples: Int = dynamicSamplesDefault): Rational = {
val inputRanges: Map[Identifier, Interval] = inputValMap.map({
case (id, i) =>
(id, Interval(i.mid - inputRangeFactor * i.radius,
i.mid + inputRangeFactor * i.radius))
})
val sampler = new Uniform(inputRanges, 485793) //no seed = System millis
// TODO: no need for the ErrorMeasurer
val measurer = new ErrorMeasurerMPFR()
var currentMaxAbsMPFR = measurer.maxAbsError
var currentMaxAbs: Double = measurer.maxAbsError.doubleValue
var i = 0
while (i < dynamicSamples) {
i = i + 1
// roundoff errors on inputs
val dblInputs: Map[Identifier, Double] = sampler.next
val mpfrInputs: Map[Identifier, MPFRFloat] = dblInputs.map({
case (x, d) => (x -> MPFRFloat.fromString(d.toString))
})
val dblOutput: Double = evalDouble(expr, dblInputs)
val mpfrOutput: MPFRFloat = evalMPFR(expr, mpfrInputs)
measurer.nextValues(dblOutput, mpfrOutput)
// Invariant that absolute errors have to grow monotonically
assert(currentMaxAbsMPFR <= measurer.maxAbsError)
currentMaxAbsMPFR = measurer.maxAbsError
assert(currentMaxAbs <= measurer.maxAbsError.doubleValue)
currentMaxAbs = measurer.maxAbsError.doubleValue
}
Rational.fromString(measurer.maxAbsError.toString())
}
}
\ No newline at end of file
This diff is collapsed.
......@@ -51,7 +51,21 @@ object FinitePrecision {
}
}
sealed abstract class Precision {
private def allPrec: List[Precision] = List(Float32, Float64, DoubleDouble, QuadDouble)
def getUpperBound(lhs: Precision, rhs: Precision): Precision = (lhs, rhs) match {
case (Fixed(a), Fixed(b)) if (a == b) => lhs
case (Fixed(a), Fixed(b)) =>
throw new Exception("mixed-precision currently unsupported for fixed-points")
case _ =>
if (allPrec.indexOf(lhs) <= allPrec.indexOf(rhs)) {
rhs
} else {
lhs
}
}
sealed abstract class Precision extends Ordered[Precision] {
/* The range of values that are representable by this precision. */
def range: (Rational, Rational)
......@@ -81,6 +95,12 @@ object FinitePrecision {
// PERFORMANCE: this may not be the fastest way
Rational.fromDouble(math.ulp(Rational.abs(r).floatValue)/2)
}
def compare(that: Precision): Int = that match {
case Float32 => 0
case Float64 => -1
case DoubleDouble => -1
}
}
case object Float64 extends Precision {
......@@ -97,6 +117,12 @@ object FinitePrecision {
def absRoundoff(r: Rational): Rational = {
Rational.fromDouble(math.ulp(Rational.abs(r).doubleValue)/2)
}
def compare(that: Precision): Int = that match {
case Float32 => 1
case Float64 => 0
case DoubleDouble => -1
}
}
case object DoubleDouble extends Precision {
......@@ -115,6 +141,12 @@ object FinitePrecision {
def absRoundoff(r: Rational): Rational = {
doubleDoubleEps * Rational.abs(r)
}
def compare(that: Precision): Int = that match {
case Float32 => 1
case Float64 => 1
case DoubleDouble => 0
}
}
case object QuadDouble extends Precision {
......@@ -133,6 +165,7 @@ object FinitePrecision {
def absRoundoff(r: Rational): Rational = {
quadDoubleEps * Rational.abs(r)
}
def compare(that: Precision): Int = ???
}
/*
......@@ -170,6 +203,9 @@ object FinitePrecision {
// TODO: don't we have to also subtract 1 for the sign?
32 - Integer.numberOfLeadingZeros(value)
}
def compare(that: Precision): Int = that match {
case Fixed(x) => bitlength.compare(x)
}
}
}
\ No newline at end of file
package daisy
package utils
import analysis._
import analysis.DynamicPhase._
import lang.Trees._
import lang.Identifiers._
import lang.NumAnnotation
import Interval._
import Rational._
import FinitePrecision._
import daisy.analysis.Sampler._
import MPFRFloat.{abs => mpfr_abs, max => mpfr_max, min => mpfr_min}
import scala.collection.immutable._
trait IntervalSubdivision extends RoundoffEvaluators {
/**
* Performs an interval subdivision, recording only the result's range
* and error.
*
* The tree does not get annotated.
* The initErrorMap only contains initial errors, NOT roundoffs.
*
* @param expr expression to be analysed
* @param inputParams free variables of the expr
* @param initValMap ranges of free variables
* @param initErrorMap initial errors of free variables (if any)
* These do NOT include roundoff errors
* @param trackRoundoffErrs whether to track roundoff errors or to only
* propagate initial errors (as provided)
* @param uniformPrecision the default precision to use when doing absolute
* errors calculations
*/
def doIntervalSubdivision(
expr: Expr,
inputParams: Set[Identifier],
initValMap: Map[Identifier, Interval],
initErrorMap: Map[Identifier, Rational],
//trackInitialErrs: Boolean,
trackRoundoffErrs: Boolean,
uniformPrecision: Precision): (Interval, Rational) = {
// TODO: this should not be hardcoded
val numSplits = 10
val inputsSubdiv: Seq[Map[Identifier, Interval]] = initValMap.foldLeft(Seq(Map[Identifier, Interval]()))({
case (currSeq: Seq[Map[Identifier, Interval]], (id, intrvl)) =>
val xlo = intrvl.xlo
val splitWidth = intrvl.width / numSplits
val splits: Seq[Interval] = (0 until numSplits).map(i =>
Interval(xlo + i * splitWidth, xlo + (i+1) * splitWidth))
currSeq.flatMap( m =>
splits.map(i => m + (id -> i))
)
})
val allIDs = inputParams
val missingIDs = allIDs -- initErrorMap.keySet
// the initErrorMap will simply be empty, if initial errors are not being tracked
def getInputErrors(rangeMap: Map[Identifier, Interval]): Map[Identifier, Rational] =
if (trackRoundoffErrs){
initErrorMap ++ missingIDs.map(id => (id -> uniformPrecision.absRoundoff(rangeMap(id))))
} else {
initErrorMap ++ missingIDs.map( id => (id -> zero))
}
val (resError, resInterval) = uniformRoundoff_IA_AA(expr,
inputsSubdiv.head, getInputErrors(inputsSubdiv.head), uniformPrecision)
// val (resInterval, resError) = evaluate[Interval, AffineForm](
// expr, inputsSubdiv.head,
// getInputErrors(inputsSubdiv.head).map(x => (x._1 -> AffineForm.fromError(x._2))),
// Interval.apply, AffineForm.zero, AffineForm.fromError, AffineForm.apply,
// false, trackRoundoffErrs, uniformPrecision)
var maxInterval = resInterval
var maxError = resError
for (m <- inputsSubdiv.tail) {
val (tmpError, tmpRange) = uniformRoundoff_IA_AA(expr, m,
getInputErrors(m), uniformPrecision)
// evaluate[Interval, AffineForm](
// expr, m,
// getInputErrors(m).map(x => (x._1 -> AffineForm.fromError(x._2))),
// Interval.apply, AffineForm.zero, AffineForm.fromError, AffineForm.apply,
// false, trackRoundoffErrs, uniformPrecision)
maxInterval = maxInterval.union(tmpRange)
maxError = max(maxError, tmpError)
}
(maxInterval, maxError)
}
/**
* Uses interval subdivision to compute ranges and affine arithmetic for errors.
*
* Unlike the other subdivision method, this one annotates the trees with