class DupireLocalVolatilityCalculator implements LocalVolatilityCalculator {

  @Override
  public DeformedSurface localVolatilityFromPrice(
      Surface callPriceSurface,
      double spot,
      Function<Double, Double> interestRate,
      Function<Double, Double> dividendRate) {

    Function<DoublesPair, ValueDerivatives> func = new Function<DoublesPair, ValueDerivatives>() {
      @Override
      public ValueDerivatives apply(DoublesPair x) {
        double t = x.getFirst();
        double k = x.getSecond();
        double r = interestRate.apply(t);
        double q = dividendRate.apply(t);
        double price = callPriceSurface.zValue(t, k);
        DoubleArray priceSensi = callPriceSurface.zValueParameterSensitivity(t, k).getSensitivity();
        double divT = FIRST_DERIV.differentiate(u -> callPriceSurface.zValue(u, k)).apply(t);
        DoubleArray divTSensi = FIRST_DERIV_SENSI.differentiate(
            u -> callPriceSurface.zValueParameterSensitivity(u.get(0), k).getSensitivity())
            .apply(DoubleArray.of(t)).column(0);
        double divK = FIRST_DERIV.differentiate(l -> callPriceSurface.zValue(t, l)).apply(k);
        DoubleArray divKSensi = FIRST_DERIV_SENSI.differentiate(
            l -> callPriceSurface.zValueParameterSensitivity(t, l.get(0)).getSensitivity())
            .apply(DoubleArray.of(k)).column(0);
        double divK2 = SECOND_DERIV.differentiate(l -> callPriceSurface.zValue(t, l)).apply(k);
        DoubleArray divK2Sensi = SECOND_DERIV_SENSI.differentiateNoCross(
            l -> callPriceSurface.zValueParameterSensitivity(t, l.get(0)).getSensitivity())
            .apply(DoubleArray.of(k)).column(0);
        double var = 2d * (divT + q * price + (r - q) * k * divK) / (k * k * divK2);
        if (var < 0d) {
          throw new IllegalArgumentException("Negative variance");
        }
        double localVol = Math.sqrt(var);
        double factor = 1d / (localVol * k * k * divK2);
        DoubleArray localVolSensi = divTSensi.multipliedBy(factor)
            .plus(divKSensi.multipliedBy((r - q) * k * factor))
            .plus(priceSensi.multipliedBy(q * factor))
            .plus(divK2Sensi.multipliedBy(-0.5 * localVol / divK2));
        return ValueDerivatives.of(localVol, localVolSensi);
      }
    };
    SurfaceMetadata metadata = DefaultSurfaceMetadata.builder()
        .xValueType(ValueType.YEAR_FRACTION)
        .yValueType(ValueType.STRIKE)
        .zValueType(ValueType.LOCAL_VOLATILITY)
        .surfaceName(SurfaceName.of("localVol_" + callPriceSurface.getName()))
        .build();
    return DeformedSurface.of(metadata, callPriceSurface, func);
  }

}