/*
    * Copyright (C) 2018 Alberto Irurueta Carro (alberto@irurueta.com)
    *
    * Licensed under the Apache License, Version 2.0 (the "License");
    * you may not use this file except in compliance with the License.
    * You may obtain a copy of the License at
    *
    *         http://www.apache.org/licenses/LICENSE-2.0
    *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package com.irurueta.navigation.indoor.radiosource;
  
  import com.irurueta.algebra.AlgebraException;
  import com.irurueta.algebra.Matrix;
  import com.irurueta.geometry.Point;
  import com.irurueta.navigation.LockedException;
  import com.irurueta.navigation.NotReadyException;
  import com.irurueta.navigation.geodesic.Constants;
  import com.irurueta.navigation.indoor.RadioSource;
  import com.irurueta.navigation.indoor.RssiReadingLocated;
  import com.irurueta.navigation.indoor.Utils;
  import com.irurueta.numerical.NumericalException;
  import com.irurueta.numerical.fitting.FittingException;
  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFitter;
  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFunctionEvaluator;
  
  import java.util.List;
  
  /**
   * Estimates position, transmitted power and path loss exponent of a
   * radio source (e.g. Wi-Fi access point or bluetooth beacon) assuming
   * that the radio source emits isotropically following the expression
   * below:
   * Pr = Pt*Gt*Gr*lambda^2 / (4*pi*d)^2,
   * where Pr is the received power (expressed in mW),
   * Gt is the Gain of the transmission antenna
   * Gr is the Gain of the receiver antenna
   * d is the distance between emitter and receiver
   * and lambda is the wavelength and is equal to: lambda = c / f,
   * where c is the speed of light
   * and f is the carrier frequency of the radio signal.
   * Because usually information about the antenna of the radio source cannot be
   * retrieved (because many measurements are made on unknown devices where
   * physical access is not possible), this implementation will estimate the
   * equivalent transmitted power as: Pte = Pt * Gt * Gr.
   * If Readings contain RSSI standard deviations, those values will be used,
   * otherwise it will be assumed an RSSI standard deviation of 1 dB.
   * <p>
   * IMPORTANT: Implementations of this class can choose to estimate a
   * combination of radio source position, transmitted power and path loss
   * exponent. However enabling all three estimations usually achieves
   * inaccurate results. When using this class, estimation must be of at least
   * one parameter (position, transmitted power or path loss exponent) when
   * initial values are provided for the other two, and at most it should consist
   * of two parameters (either position and transmitted power, position and
   * path loss exponent or transmitted power and path loss exponent), providing an
   * initial value for the remaining parameter.
   *
   * @param <S> a {@link RadioSource} type.
   * @param <P> a {@link Point} type.
   */
  @SuppressWarnings("DuplicatedCode")
  public abstract class RssiRadioSourceEstimator<S extends RadioSource, P extends Point<P>>
          extends RadioSourceEstimator<P, RssiReadingLocated<S, P>, RssiRadioSourceEstimatorListener<S, P>> {
  
      /**
       * Speed of light expressed in meters per second (m/s).
       */
      public static final double SPEED_OF_LIGHT = Constants.SPEED_OF_LIGHT;
  
      /**
       * Default standard deviations assumed for RSSI readings being fitted.
       */
      public static final double DEFAULT_POWER_STANDARD_DEVIATION = 1.0;
  
      /**
       * Default exponent typically used on free space for path loss propagation in
       * terms of distance. This value is used for free space environments.
       */
      public static final double DEFAULT_PATH_LOSS_EXPONENT = 2.0;
  
      /**
       * Indicates whether radio source position estimation is enabled or not by default.
       */
      public static final boolean DEFAULT_POSITION_ESTIMATION_ENABLED = true;
  
      /**
       * Indicates whether radio source transmitted power estimation is enabled or not by
       * default. Typically, this data is required for Wi-Fi Access points, but it is already
       * provided for Beacons (and hence its estimation is not needed).
       */
      public static final boolean DEFAULT_TRANSMITTED_POWER_ESTIMATION_ENABLED = true;
  
      /**
      * Indicates whether path loss estimation is enabled or not by default.
      */
     public static final boolean DEFAULT_PATHLOSS_ESTIMATION_ENABLED = false;
 
 
     /**
      * Indicates whether radio source position estimation is enabled or not.
      */
     private boolean positionEstimationEnabled = DEFAULT_POSITION_ESTIMATION_ENABLED;
 
     /**
      * Estimated transmitted power expressed in dBm's.
      */
     private double estimatedTransmittedPowerdBm;
 
     /**
      * Indicates whether transmitted power estimation is enabled or not.
      */
     private boolean transmittedPowerEstimationEnabled = DEFAULT_TRANSMITTED_POWER_ESTIMATION_ENABLED;
 
     /**
      * Estimated exponent typically used on free space for path loss propagation in
      * terms of distance.
      * On different environments path loss exponent might have different values:
      * - Free space: 2.0
      * - Urban Area: 2.7 to 3.5
      * - Suburban Area: 3 to 5
      * - Indoor (line-of-sight): 1.6 to 1.8
      * If path loss exponent estimation is not enabled, this value will always be equal to
      * {@link #DEFAULT_PATH_LOSS_EXPONENT}
      */
     private double estimatedPathLossExponent = DEFAULT_PATH_LOSS_EXPONENT;
 
     /**
      * Variance of estimated transmitted power.
      * This value will only be available when transmitted power
      * estimation is enabled.
      */
     private Double estimatedTransmittedPowerVariance;
 
     /**
      * Variance of estimated path loss exponent.
      * This value will only be available when path-loss
      * exponent estimation is enabled.
      */
     private Double estimatedPathLossExponentVariance;
 
     /**
      * Estimated chi square value.
      */
     private double chiSq;
 
     /**
      * Initial transmitted power to start the estimation of radio source
      * transmitted power.
      * If not defined, average value of received power readings will be used.
      */
     private Double initialTransmittedPowerdBm;
 
     /**
      * Initial position to start the estimation of radio source position.
      * If not defined, centroid of provided readings will be used.
      */
     private P initialPosition;
 
     /**
      * Initial exponent typically used on free space for path loss propagation in
      * terms of distance.
      * On different environments path loss exponent might have different values:
      * - Free space: 2.0
      * - Urban Area: 2.7 to 3.5
      * - Suburban Area: 3 to 5
      * - Indoor (line-of-sight): 1.6 to 1.8
      * <p>
      * If path loss exponent estimation is enabled, estimation will start at this
      * value and will converge to the most appropriate value.
      * If path loss exponent estimation is disabled, this value will be assumed
      * to be exact and the estimated path loss exponent will be equal to this
      * value.
      */
     private double initialPathLossExponent = DEFAULT_PATH_LOSS_EXPONENT;
 
     /**
      * Indicates whether path loss estimation is enabled or not.
      */
     private boolean pathLossEstimationEnabled = DEFAULT_PATHLOSS_ESTIMATION_ENABLED;
 
     /**
      * Levenberg-Marquardt fitter to find a solution.
      */
     private final LevenbergMarquardtMultiDimensionFitter fitter = new LevenbergMarquardtMultiDimensionFitter();
 
     /**
      * Constructor.
      */
     protected RssiRadioSourceEstimator() {
         super();
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings radio signal readings belonging to the same
      *                 radio source.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(final List<? extends RssiReadingLocated<S, P>> readings) {
         super(readings);
     }
 
     /**
      * Constructor.
      *
      * @param listener listener in charge of attending events raised by this instance.
      */
     protected RssiRadioSourceEstimator(final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(listener);
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings radio signal readings belonging to the same radio source.
      * @param listener listener in charge of attending events raised by this instance.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings,
             final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(readings, listener);
     }
 
     /**
      * Constructor.
      *
      * @param initialPosition initial position to start the estimation of radio
      *                        source position.
      */
     protected RssiRadioSourceEstimator(final P initialPosition) {
         this.initialPosition = initialPosition;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings        radio signal readings belonging to the same radio source.
      * @param initialPosition initial position to start the estimation of radio
      *                        source position.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final P initialPosition) {
         super(readings);
         this.initialPosition = initialPosition;
     }
 
     /**
      * Constructor.
      *
      * @param initialPosition initial position to start the estimation of radio
      *                        source position.
      * @param listener        listener in charge of attending events raised by this instance.
      */
     protected RssiRadioSourceEstimator(final P initialPosition, final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(listener);
         this.initialPosition = initialPosition;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings        radio signal readings belonging to the same radio source.
      * @param initialPosition initial position to start the estimation of radio
      *                        source position.
      * @param listener        listener in charge of attending events raised by this instance.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final P initialPosition,
             final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(readings, listener);
         this.initialPosition = initialPosition;
     }
 
     /**
      * Constructor.
      *
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's).
      */
     protected RssiRadioSourceEstimator(final Double initialTransmittedPowerdBm) {
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings                   radio signal readings belonging to the same radio source.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's).
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final Double initialTransmittedPowerdBm) {
         super(readings);
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      *
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      * @param listener                   listener in charge of attending events raised by this instance.
      */
     protected RssiRadioSourceEstimator(
             final Double initialTransmittedPowerdBm, final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(listener);
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings                   radio signal readings belonging to the same radio source.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      * @param listener                   listener in charge of attending events raised by this instance.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final Double initialTransmittedPowerdBm,
             final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(readings, listener);
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings                   radio signal readings belonging to the same radio source.
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final P initialPosition,
             final Double initialTransmittedPowerdBm) {
         super(readings);
         this.initialPosition = initialPosition;
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      *
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      */
     protected RssiRadioSourceEstimator(final P initialPosition, final Double initialTransmittedPowerdBm) {
         this.initialPosition = initialPosition;
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      *
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      * @param listener                   listener in charge of attending events raised by this instance.
      */
     protected RssiRadioSourceEstimator(
             final P initialPosition, final Double initialTransmittedPowerdBm,
             final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(listener);
         this.initialPosition = initialPosition;
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings                   radio signal readings belonging to the same radio source.
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      * @param listener                   listener in charge of attending events raised by this instance.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final P initialPosition,
             final Double initialTransmittedPowerdBm, final RssiRadioSourceEstimatorListener<S, P> listener) {
         super(readings, listener);
         this.initialPosition = initialPosition;
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings                   radio signal readings belonging to the same radio source.
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's).
      * @param initialPathLossExponent    initial path loss exponent. A typical value is 2.0.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final P initialPosition,
             final Double initialTransmittedPowerdBm, final double initialPathLossExponent) {
         this(readings, initialPosition, initialTransmittedPowerdBm);
         this.initialPathLossExponent = initialPathLossExponent;
     }
 
     /**
      * Constructor.
      *
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's).
      * @param initialPathLossExponent    initial path loss exponent. A typical value is 2.0.
      */
     protected RssiRadioSourceEstimator(
             final P initialPosition, final Double initialTransmittedPowerdBm, final double initialPathLossExponent) {
         this(initialPosition, initialTransmittedPowerdBm);
         this.initialPathLossExponent = initialPathLossExponent;
     }
 
     /**
      * Constructor.
      *
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      * @param initialPathLossExponent    initial path loss exponent. A typical value is 2.0.
      * @param listener                   listener in charge of attending events raised by this instance.
      */
     protected RssiRadioSourceEstimator(
             final P initialPosition, final Double initialTransmittedPowerdBm, final double initialPathLossExponent,
             final RssiRadioSourceEstimatorListener<S, P> listener) {
         this(initialPosition, initialTransmittedPowerdBm, listener);
         this.initialPathLossExponent = initialPathLossExponent;
     }
 
     /**
      * Constructor.
      * Sets radio signal readings belonging to the same radio source.
      *
      * @param readings                   radio signal readings belonging to the same radio source.
      * @param initialPosition            initial position to start the estimation of radio
      *                                   source position.
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted power
      *                                   (expressed in dBm's)
      * @param initialPathLossExponent    initial path loss exponent. A typical value is 2.0.
      * @param listener                   listener in charge of attending events raised by this instance.
      * @throws IllegalArgumentException if readings are not valid.
      */
     protected RssiRadioSourceEstimator(
             final List<? extends RssiReadingLocated<S, P>> readings, final P initialPosition,
             final Double initialTransmittedPowerdBm, final double initialPathLossExponent,
             final RssiRadioSourceEstimatorListener<S, P> listener) {
         this(readings, initialPosition, initialTransmittedPowerdBm, listener);
         this.initialPathLossExponent = initialPathLossExponent;
     }
 
     /**
      * Gets initial transmitted power to start the estimation of radio source
      * transmitted power (expressed in dBm's).
      * If not defined, average value of received power readings will be used.
      * <p>
      * If transmitted power estimation is enabled, estimation will start at this
      * value and will be converted to the most appropriate value.
      * If transmitted power estimation is disabled, this value will be assumed to be
      * exact and the estimated transmitted power will be equal to this value
      * (converted to dBm's).
      *
      * @return initial transmitted power to start the estimation of radio source
      * transmitted power.
      */
     public Double getInitialTransmittedPowerdBm() {
         return initialTransmittedPowerdBm;
     }
 
     /**
      * Sets initial transmitted power to start the estimation of radio source
      * transmitted power (expressed in dBm's).
      * If not defined, average value of received power readings will be used.
      * <p>
      * If transmitted power estimation is enabled, estimation will start at this
      * value and will be converted to the most appropriate value.
      * If transmitted power estimation is disabled, this value will be assumed to be
      * exact and the estimated transmitted power will be equal to this value
      * (converted to dBm's).
      *
      * @param initialTransmittedPowerdBm initial transmitted power to start the
      *                                   estimation of radio source transmitted
      *                                   power.
      * @throws LockedException if estimator is locked.
      */
     public void setInitialTransmittedPowerdBm(final Double initialTransmittedPowerdBm) throws LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         this.initialTransmittedPowerdBm = initialTransmittedPowerdBm;
     }
 
     /**
      * Gets initial transmitted power to start the estimation of radio source
      * transmitted power (expressed in mW).
      * If not defined, average value of received power readings will be used.
      * <p>
      * If transmitted power estimation is enabled, estimation will start at this
      * value and will be converted to the most appropriate value.
      * If transmitted power estimation is disabled, this value will be assumed to be
      * exact and the estimated transmitted power will be equal to this value
      * (converted to dBm's).
      *
      * @return initial transmitted power to start the estimation of radio source
      * transmitted power.
      */
     public Double getInitialTransmittedPower() {
         return initialTransmittedPowerdBm != null ? Utils.dBmToPower(initialTransmittedPowerdBm) : null;
     }
 
     /**
      * Sets initial transmitted power to start the estimation of radio source
      * transmitted power (expressed in mW).
      * If not defined, average value of received power readings will be used.
      * <p>
      * If transmitted power estimation is enabled, estimation will start at this
      * value and will be converted to the most appropriate value.
      * If transmitted power estimation is disabled, this value will be assumed to be
      * exact and the estimated transmitted power will be equal to this value
      * (converted to dBm's).
      *
      * @param initialTransmittedPower initial transmitted power to start the
      *                                estimation of radio source transmitted power.
      * @throws LockedException          if estimator is locked.
      * @throws IllegalArgumentException if provided value is negative.
      */
     public void setInitialTransmittedPower(final Double initialTransmittedPower) throws LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         if (initialTransmittedPower != null) {
             if (initialTransmittedPower < 0.0) {
                 throw new IllegalArgumentException();
             }
             initialTransmittedPowerdBm = Utils.powerTodBm(initialTransmittedPower);
         } else {
             initialTransmittedPowerdBm = null;
         }
     }
 
     /**
      * Indicates whether transmitted power estimation is enabled or not.
      *
      * @return true if transmitted power estimation is enabled, false otherwise.
      */
     public boolean isTransmittedPowerEstimationEnabled() {
         return transmittedPowerEstimationEnabled;
     }
 
     /**
      * Specifies whether transmitted power estimation is enabled or not.
      *
      * @param transmittedPowerEstimationEnabled true if transmitted power estimation is enabled,
      *                                          false otherwise.
      * @throws LockedException if estimator is locked.
      */
     public void setTransmittedPowerEstimationEnabled(final boolean transmittedPowerEstimationEnabled)
             throws LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         this.transmittedPowerEstimationEnabled = transmittedPowerEstimationEnabled;
     }
 
     /**
      * Gets initial position to start the estimation of radio source position.
      * If not defined, centroid of provided readings will be used.
      * <p>
      * If position estimation is enabled, estimation will start at this value
      * and will converge to the most appropriate value.
      * If position estimation is disabled, this value will be assumed to
      * be exact and the estimated position will be equal to this value.
      *
      * @return initial position to start the estimation of radio source position.
      */
     public P getInitialPosition() {
         return initialPosition;
     }
 
     /**
      * Sets initial position to start the estimation of radio source position.
      * If not defined, centroid of provided fingerprints will be used.
      * <p>
      * If position estimation is enabled, estimation will start at this value
      * and will converge to the most appropriate value.
      * If position estimation is disabled, this value will be assumed to
      * be exact and the estimated position will be equal to this value.
      *
      * @param initialPosition initial position to start the estimation of radio
      *                        source position.
      * @throws LockedException if estimator is locked.
      */
     public void setInitialPosition(final P initialPosition) throws LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         this.initialPosition = initialPosition;
     }
 
     /**
      * Indicates whether radio source position estimation is enabled or not.
      *
      * @return true if position estimation is enabled, false otherwise.
      */
     public boolean isPositionEstimationEnabled() {
         return positionEstimationEnabled;
     }
 
     /**
      * Specifies whether radio source position estimation is enabled or not.
      *
      * @param positionEstimationEnabled true if position estimation is enabled,
      *                                  false otherwise.
      * @throws LockedException if estimator is locked.
      */
     public void setPositionEstimationEnabled(final boolean positionEstimationEnabled) throws LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         this.positionEstimationEnabled = positionEstimationEnabled;
     }
 
     /**
      * Gets initial exponent typically used on free space for path loss propagation
      * in terms of distance.
      * On different environments path loss exponent might have different value:
      * - Free space: 2.0
      * - Urban Area: 2.7 to 3.5
      * - Suburban Area: 3 to 5
      * - Indoor (line-of-sight): 1.6 to 1.8
      * <p>
      * If path loss exponent estimation is enabled, estimation will start at this
      * value and will converge to the most appropriate value.
      * If path loss exponent estimation is disabled, this value will be assumed
      * to be exact and the estimated path loss exponent will be equal to this
      * value.
      *
      * @return initial path loss exponent.
      */
     public double getInitialPathLossExponent() {
         return initialPathLossExponent;
     }
 
     /**
      * Sets initial exponent typically used on free space for path loss propagation
      * in terms of distance.
      * On different environments path loss exponent might have different value:
      * - Free space: 2.0
      * - Urban Area: 2.7 to 3.5
      * - Suburban Area: 3 to 5
      * - Indoor (line-of-sight): 1.6 to 1.8
      * <p>
      * If path loss exponent estimation is enabled, estimation will start at this
      * value and will converge to the most appropriate value.
      * If path loss exponent estimation is disabled, this value will be assumed
      * to be exact and the estimated path loss exponent will be equal to this
      * value.
      *
      * @param initialPathLossExponent initial path loss exponent.
      * @throws LockedException if estimator is locked.
      */
     public void setInitialPathLossExponent(final double initialPathLossExponent) throws LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         this.initialPathLossExponent = initialPathLossExponent;
     }
 
     /**
      * Indicates whether path loss estimation is enabled or not.
      *
      * @return true if path loss estimation is enabled, false otherwise.
      */
     public boolean isPathLossEstimationEnabled() {
         return pathLossEstimationEnabled;
     }
 
     /**
      * Specifies whether path loss estimation is enabled or not.
      *
      * @param pathLossEstimationEnabled true if path loss estimation is enabled,
      *                                  false otherwise.
      * @throws LockedException if estimator is locked.
      */
     public void setPathLossEstimationEnabled(final boolean pathLossEstimationEnabled) throws LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         this.pathLossEstimationEnabled = pathLossEstimationEnabled;
     }
 
     /**
      * Indicates whether this instance is ready to start the estimation.
      *
      * @return true if this instance is ready, false otherwise.
      */
     @Override
     public boolean isReady() {
         // at least one parameter estimation must be enabled
         return (positionEstimationEnabled || transmittedPowerEstimationEnabled || pathLossEstimationEnabled)
                 // if position estimation is disabled, an initial position must be provided
                 && !(!positionEstimationEnabled && initialPosition == null)
                 // if transmitted power estimation is disabled, an initial transmitted power must be provided
                 && !(!transmittedPowerEstimationEnabled && initialTransmittedPowerdBm == null)
                 // readings must also be valid
                 && areValidReadings(readings);
     }
 
     /**
      * Estimate position, transmitted power and path loss exponent.
      *
      * @throws RadioSourceEstimationException if estimation fails.
      * @throws NotReadyException              if estimator is not ready.
      * @throws LockedException                if estimator is locked.
      */
     @SuppressWarnings("all")
     @Override
     public void estimate() throws RadioSourceEstimationException, NotReadyException, LockedException {
         if (isLocked()) {
             throw new LockedException();
         }
         if (!isReady()) {
             throw new NotReadyException();
         }
 
         try {
             locked = true;
 
             if (listener != null) {
                 listener.onEstimateStart(this);
             }
 
             if (positionEstimationEnabled && !transmittedPowerEstimationEnabled && !pathLossEstimationEnabled) {
                 // only position estimation is enabled
                 setupFitterPosition();
             } else if (!positionEstimationEnabled && transmittedPowerEstimationEnabled && !pathLossEstimationEnabled) {
                 // only transmitted power estimation is enabled
                 setupFitterTransmittedPower();
             } else if (!positionEstimationEnabled && !transmittedPowerEstimationEnabled && pathLossEstimationEnabled) {
                 // only pathloss estimation is enabled
                 setupFitterPathLossExponent();
             } else if (positionEstimationEnabled && transmittedPowerEstimationEnabled && !pathLossEstimationEnabled) {
                 // position and transmitted power enabled
                 setupFitterPositionAndTransmittedPower();
             } else if (positionEstimationEnabled && !transmittedPowerEstimationEnabled && pathLossEstimationEnabled) {
                 // position and pathloss enabled
                 setupFitterPositionAndPathLossExponent();
             } else if (!positionEstimationEnabled && transmittedPowerEstimationEnabled && pathLossEstimationEnabled) {
                 // transmitted power and pathloss enabled
                 setupFitterTransmittedPowerAndPathLossExponent();
             } else {
                 // position, transmitted power and pathloss enabled
                 setupFitterPositionTransmittedPowerAndPathLossExponent();
             }
 
             fitter.fit();
 
             // estimated position and transmitted power
             final var a = fitter.getA();
             final var dims = getNumberOfDimensions();
 
             estimatedCovariance = fitter.getCovar();
             chiSq = fitter.getChisq();
 
             var pos = 0;
             estimatedPositionCoordinates = new double[dims];
             if (positionEstimationEnabled) {
                 // position estimation enabled
                 System.arraycopy(a, 0, estimatedPositionCoordinates, 0, dims);
 
                 if (estimatedCovariance != null) {
                     final var d = dims - 1;
                     if (estimatedPositionCovariance == null) {
                         estimatedPositionCovariance = estimatedCovariance.getSubmatrix(0, 0, d,
                                 d);
                     } else {
                         estimatedCovariance.getSubmatrix(0, 0, d, d,
                                 estimatedPositionCovariance);
                     }
                 }
                 pos += dims;
             } else {
                 // position estimation disabled
                 if (initialPosition != null) {
                     for (var i = 0; i < dims; i++) {
                         estimatedPositionCoordinates[i] = initialPosition.getInhomogeneousCoordinate(i);
                     }
                 }
 
                 estimatedPositionCovariance = null;
             }
 
             if (transmittedPowerEstimationEnabled) {
                 // transmitted power estimation enabled
                 estimatedTransmittedPowerdBm = a[pos];
 
                 if (estimatedCovariance != null) {
                     estimatedTransmittedPowerVariance = estimatedCovariance.getElementAt(pos, pos);
                 }
                 pos++;
             } else {
                 // transmitted power estimation disabled
                 if (initialTransmittedPowerdBm != null) {
                     estimatedTransmittedPowerdBm = initialTransmittedPowerdBm;
                 }
                 estimatedTransmittedPowerVariance = null;
             }
 
             if (pathLossEstimationEnabled) {
                 // pathloss exponent estimation enabled
                 estimatedPathLossExponent = a[pos];
 
                 if (estimatedCovariance != null) {
                     estimatedPathLossExponentVariance = estimatedCovariance.getElementAt(pos, pos);
                 }
             } else {
                 // pathloss exponent estimation disabled
                 estimatedPathLossExponent = initialPathLossExponent;
                 estimatedPathLossExponentVariance = null;
             }
 
 
             if (listener != null) {
                 listener.onEstimateEnd(this);
             }
         } catch (final NumericalException e) {
             throw new RadioSourceEstimationException(e);
         } finally {
             locked = false;
         }
     }
 
     /**
      * Gets estimated transmitted power expressed in milli watts (mW).
      *
      * @return estimated transmitted power expressed in milli watts.
      */
     public double getEstimatedTransmittedPower() {
         return Utils.dBmToPower(estimatedTransmittedPowerdBm);
     }
 
     /**
      * Gets estimated transmitted power expressed in dBm's.
      *
      * @return estimated transmitted power expressed in dBm's.
      */
     public double getEstimatedTransmittedPowerdBm() {
         return estimatedTransmittedPowerdBm;
     }
 
     /**
      * Gets estimated exponent typically used on free space for path loss propagation in
      * terms of distance.
      * On different environments path loss exponent might have different values:
      * - Free space: 2.0
      * - Urban Area: 2.7 to 3.5
      * - Suburban Area: 3 to 5
      * - Indoor (line-of-sight): 1.6 to 1.8
      * If path loss exponent estimation is not enabled, this value will always be equal to
      * {@link #DEFAULT_PATH_LOSS_EXPONENT}
      *
      * @return estimated path loss exponent.
      */
     public double getEstimatedPathLossExponent() {
         return estimatedPathLossExponent;
     }
 
     /**
      * Gets estimated transmitted power variance.
      * This value will only be available when transmitted power
      * estimation is enabled.
      *
      * @return estimated transmitted power variance or null.
      */
     public Double getEstimatedTransmittedPowerVariance() {
         return estimatedTransmittedPowerVariance;
     }
 
     /**
      * Gets estimated path loss exponent variance.
      * This value will only be available when path-loss
      * exponent estimation is enabled.
      *
      * @return estimated path loss exponent variance or null.
      */
     public Double getEstimatedPathLossExponentVariance() {
         return estimatedPathLossExponentVariance;
     }
 
     /**
      * Gets estimated chi square value.
      *
      * @return estimated chi square value.
      */
     public double getChiSq() {
         return chiSq;
     }
 
     /**
      * Setups fitter to estimated position.
      *
      * @throws FittingException if Levenberg-Marquardt fitting fails.
      */
     private void setupFitterPosition() throws FittingException {
         // because all readings must belong to the same radio source, we
         // obtain the frequency of the first radio source on the first reading
         var reading = readings.get(0);
         final var frequency = reading.getSource().getFrequency();
 
         // n = 2.0, is the path loss exponent (which is typically 2.0)
 
         // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
         // lambda = c/f, where lambda is wavelength,
         // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the received Gain
         // Pr = Pte*c^n/((4*pi*f)^n * d^n)
 
 
         // compute k as the constant part of the isotropic received power formula
         // so that: Pr = Pte*k/d^n
         final var k = Math.pow(SPEED_OF_LIGHT / (4.0 * Math.PI * frequency), initialPathLossExponent);
         final var kdB = 10.0 * Math.log10(k);
 
         final var dims = getNumberOfDimensions();
         final var initialTransmittedPowerdBm = computeInitialTransmittedPowerdBm();
 
         // for numerical accuracy reasons, a logarithmic version of the previous
         // formula will be used instead
         // Pr (dBm) = 10 * log(Pte * k / d^n) = 10*log(k) + 10*log(Pte) - 10*n*log(d)
 
         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
             @Override
             public int getNumberOfDimensions() {
                 return dims;
             }
 
             @Override
             public double[] createInitialParametersArray() {
                 final var initial = new double[dims];
                 computeInitialPosition(initial, dims);
                 return initial;
             }
 
             @Override
             public double evaluate(
                     final int i, final double[] point, final double[] params, final double[] derivatives) {
                 var sqrDistance = 0.0;
                 for (var j = 0; j < dims; j++) {
                     final var diff = params[j] - point[j];
                     sqrDistance += diff * diff;
 
                     // n is mInitialPathLossExponent, which is typically 2.0
                     derivatives[j] = -10.0 * initialPathLossExponent * diff;
                 }
 
                 // derivatives respect position coordinates are (2D case):
                 // f(x,y) = -5*n*log((x - xap)^2 + (y - yap)^2)
                 // df/dx = -5*n*2*(x - xap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffX/(ln(10)*sqrDistance)
                 // df/dy = -5*n*2*(y - yap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffY/(ln(10)*sqrDistance)
                 final var ln10PerSqrDistance = Math.log(10.0) * sqrDistance;
                 if (ln10PerSqrDistance != 0.0) {
                     for (var j = 0; j < dims; j++) {
                         derivatives[j] /= ln10PerSqrDistance;
                     }
                 }
 
                 // d^2 = (x - xap)^2 + (y - yap)^2
                 // d^n = (d^2)^n/2
 
                 // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
                 // n is the path loss exponent
                 // lambda = c/f, where lambda is wavelength,
                 // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the
                 //           received Gain
                 // Pr = Pte*c^n/((4*pi*f)^n * d^n)
                 // Pr (dBm) = 10*log(k) + 10*log(Pte) - 10*log(d^n) =
                 // 10*log(k) + 10*log(Pte) - 10*log((d^2)^n/2) =
                 // 10*log(k) + 10*log(Pte) - 10*n/2*log(d^2) =
                 // 10*log(k) + 10*log(Pte) - 5*n*log(d^2) =
                 return kdB + initialTransmittedPowerdBm - 5.0 * initialPathLossExponent * Math.log10(sqrDistance);
             }
         });
 
         final var numReadings = readings.size();
         try {
             final var x = new Matrix(numReadings, dims);
             final var y = new double[numReadings];
             final var standardDeviations = new double[numReadings];
             for (var i = 0; i < numReadings; i++) {
                 reading = readings.get(i);
                 final var position = reading.getPosition();
 
                 for (var j = 0; j < dims; j++) {
                     x.setElementAt(i, j, position.getInhomogeneousCoordinate(j));
                 }
 
                 standardDeviations[i] = reading.getRssiStandardDeviation() != null ? reading.getRssiStandardDeviation()
                         : DEFAULT_POWER_STANDARD_DEVIATION;
                 y[i] = reading.getRssi();
             }
 
             fitter.setInputData(x, y, standardDeviations);
         } catch (final AlgebraException ignore) {
             // never happens
         }
     }
 
     /**
      * Setups fitter to estimate transmitted power.
      *
      * @throws FittingException if Levenberg-Marquardt fitting fails.
      */
     private void setupFitterTransmittedPower() throws FittingException {
         // because all readings must belong to the same radio source, we
         // obtain the frequency of the first radio source on the first reading
         var reading = readings.get(0);
         final var frequency = reading.getSource().getFrequency();
 
         // n = 2.0, is the path loss exponent (which is typically 2.0)
 
         // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
         // lambda = c/f, where lambda is wavelength,
         // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the received Gain
         // Pr = Pte*c^n/((4*pi*f)^n * d^n)
 
 
         // compute k as the constant part of the isotropic received power formula
         // so that: Pr = Pte*k/d^n
         final var k = Math.pow(SPEED_OF_LIGHT / (4.0 * Math.PI * frequency), initialPathLossExponent);
         final var kdB = 10.0 * Math.log10(k);
         final var initialTransmittedPowerdBm = computeInitialTransmittedPowerdBm();
 
         // for numerical accuracy reasons, a logarithmic version of the previous
         // formula will be used instead
         // Pr (dBm) = 10 * log(Pte * k / d^n) = 10*log(k) + 10*log(Pte) - 10*n*log(d)
 
         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
             @Override
             public int getNumberOfDimensions() {
                 return 1;
             }
 
             @Override
             public double[] createInitialParametersArray() {
                 final var initial = new double[1];
                 initial[0] = initialTransmittedPowerdBm;
                 return initial;
             }
 
             @Override
             public double evaluate(
                     final int i, final double[] point, final double[] params, final double[] derivatives) {
                 final var sqrDistance = initialPosition.sqrDistanceTo(readings.get(i).getPosition());
 
                 final var transmittedPowerdBm = params[0];
 
                 // derivative respect transmitted power Pt (dBm) = 10*log(Pte)
                 derivatives[0] = 1.0;
 
                 // d^2 = (x - xap)^2 + (y - yap)^2
                 // d^n = (d^2)^n/2
 
                 // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
                 // n is the path loss exponent
                 // lambda = c/f, where lambda is wavelength,
                 // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the
                 //          received Gain
                 // Pr = Pte*c^n/((4*pi*f)^n * d^n)
                 // Pr (dBm) = 10*log(k) + 10*log(Pte) - 10*log(d^n) =
                 // 10*log(k) + 10*log(Pte) - 10*log((d^2)^n/2) =
                 // 10*log(k) + 10*log(Pte) - 10*n/2*log(d^2) =
                 // 10*log(k) + 10*log(Pte) - 5*n*log(d^2) =
                 return kdB + transmittedPowerdBm - 5.0 * initialPathLossExponent * Math.log10(sqrDistance);
             }
         });
 
         final var numReadings = readings.size();
         try {
             final var x = new Matrix(numReadings, 1);
             final var y = new double[numReadings];
             final var standardDeviations = new double[numReadings];
 
             x.initialize(initialTransmittedPowerdBm);
             for (var i = 0; i < numReadings; i++) {
                 reading = readings.get(i);
 
                 standardDeviations[i] = reading.getRssiStandardDeviation() != null ? reading.getRssiStandardDeviation()
                         : DEFAULT_POWER_STANDARD_DEVIATION;
                 y[i] = reading.getRssi();
             }
 
             fitter.setInputData(x, y, standardDeviations);
         } catch (final AlgebraException ignore) {
             // never happens
         }
     }
 
     /**
      * Setups fitter to estimated path loss exponent.
      *
      * @throws FittingException if Levenberg-Marquardt fitting fails.
      */
     private void setupFitterPathLossExponent() throws FittingException {
         // because all readings must belong to the same radio source, we
         // obtain the frequency of the first radio source on the first reading
         var reading = readings.get(0);
         final var frequency = reading.getSource().getFrequency();
 
         // n = 2.0, is the path loss exponent
 
         // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
         // lambda = c/f, where lambda is wavelength,
         // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the received Gain
         // Pr = Pte*c^n/((4*pi*f)^n * d^n)
 
 
         // k is defined so that: Pr = Pte * k^n / d^n so that
         // k = (c/(4*pi*f))
         final var k = SPEED_OF_LIGHT / (4.0 * Math.PI * frequency);
         final var kdB = 10.0 * Math.log10(k);
         final var initialTransmittedPowerdBm = computeInitialTransmittedPowerdBm();
 
         // for numerical accuracy reasons, a logarithmic version of the previous
         // formula will be used instead
         // Pr (dBm) = 10 * log(Pte * k^n / d^n) = 10*n*log(k) + 10*log(Pte) - 10*n*log(d)
 
         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
             @Override
             public int getNumberOfDimensions() {
                 return 1;
             }
 
             @Override
             public double[] createInitialParametersArray() {
                 final var initial = new double[1];
                 initial[0] = initialPathLossExponent;
                 return initial;
             }
 
             @Override
             public double evaluate(
                     final int i, final double[] point, final double[] params, final double[] derivatives) {
                 final var sqrDistance = initialPosition.sqrDistanceTo(readings.get(i).getPosition());
                 final var pathLossExponent = params[0];
 
                 // derivative respect to path loss exponent
                 // f(x,y,n) = n*kdB -5*n*log((x - xap)^2 + (y - yap)^2)
                 // df/dn = kdB -5*log((x - xap)^2 + (y - yap)^2) = kdB - 5*log(sqrDistance)
                 final var logSqrDistance = Math.log10(sqrDistance);
                 derivatives[0] = kdB - 5 * logSqrDistance;
 
                 // d^2 = (x - xap)^2 + (y - yap)^2
                 // d^n = (d^2)^n/2
 
                 // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
                 // n is the path loss exponent
                 // lambda = c/f, where lambda is wavelength,
                 // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the
                 //           received Gain
                 // Pr = Pte*c^n/((4*pi*f)^n * d^n)
                 // Pr (dBm) = 10*log(k^n) + 10*log(Pte) - 10*log(d^n) =
                 // 10*log(k^n) + 10*log(Pte) - 10*log((d^2)^n/2) =
                 // 10*n*log(k) + 10*log(Pte) - 10*n/2*log(d^2) =
                 // 10*n*log(k) + 10*log(Pte) - 5*n*log(d^2) =
                 return pathLossExponent * kdB + initialTransmittedPowerdBm - 5.0 * pathLossExponent * logSqrDistance;
             }
         });
 
         final var numReadings = readings.size();
         try {
             final var x = new Matrix(numReadings, 1);
             final var y = new double[numReadings];
             final var standardDeviations = new double[numReadings];
 
             x.initialize(initialPathLossExponent);
             for (var i = 0; i < numReadings; i++) {
                 reading = readings.get(i);
 
                 standardDeviations[i] = reading.getRssiStandardDeviation() != null ? reading.getRssiStandardDeviation()
                         : DEFAULT_POWER_STANDARD_DEVIATION;
                 y[i] = reading.getRssi();
             }
 
             fitter.setInputData(x, y, standardDeviations);
         } catch (final AlgebraException ignore) {
             // never happens
         }
     }
 
     /**
      * Setups fitter to estimate transmitted power and position.
      *
      * @throws FittingException if Levenberg-Marquardt fitting fails.
      */
     private void setupFitterPositionAndTransmittedPower() throws FittingException {
         // because all readings must belong to the same radio source, we
         // obtain the frequency of the first radio source on the first reading
         var reading = readings.get(0);
         final var frequency = reading.getSource().getFrequency();
 
         // n = 2.0, is the path loss exponent (which is typically 2.0)
 
         // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
         // lambda = c/f, where lambda is wavelength,
         // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the received Gain
         // Pr = Pte*c^n/((4*pi*f)^n * d^n)
 
 
         // compute k as the constant part of the isotropic received power formula
         // so that: Pr = Pte*k/d^n
         final var k = Math.pow(SPEED_OF_LIGHT / (4.0 * Math.PI * frequency), initialPathLossExponent);
         final var kdB = 10.0 * Math.log10(k);
 
         final var dims = getNumberOfDimensions();
         final var dimsPlus1 = dims + 1;
 
         final var initialTransmittedPowerdBm = computeInitialTransmittedPowerdBm();
 
         // for numerical accuracy reasons, a logarithmic version of the previous
         // formula will be used instead
         // Pr (dBm) = 10 * log(Pte * k / d^n) = 10*log(k) + 10*log(Pte) - 10*n*log(d)
 
         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
 
             @Override
             public int getNumberOfDimensions() {
                 return dimsPlus1;
             }
 
             @Override
             public double[] createInitialParametersArray() {
                 final var initial = new double[dimsPlus1];
 
                 // initial position
                 computeInitialPosition(initial, dims);
 
                 // initial transmitted power
                 initial[dims] = initialTransmittedPowerdBm;
 
                 return initial;
             }
 
             @Override
             public double evaluate(
                     final int i, final double[] point, final double[] params, final double[] derivatives) {
                 var sqrDistance = 0.0;
                 for (var j = 0; j < dims; j++) {
                     final var diff = params[j] - point[j];
                     sqrDistance += diff * diff;
 
                     // n is mInitialPathLossExponent, which is typically 2.0
                     derivatives[j] = -10.0 * initialPathLossExponent * diff;
                 }
 
                 final var transmittedPowerdBm = params[dims];
 
                 // derivatives respect position coordinates are (2D case):
                 // f(x,y) = -5*n*log((x - xap)^2 + (y - yap)^2)
                 // df/dx = -5*n*2*(x - xap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffX/(ln(10)*sqrDistance)
                 // df/dy = -5*n*2*(y - yap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffY/(ln(10)*sqrDistance)
                 final var ln10PerSqrDistance = Math.log(10.0) * sqrDistance;
                 if (ln10PerSqrDistance != 0.0) {
                     for (var j = 0; j < dims; j++) {
                         derivatives[j] /= ln10PerSqrDistance;
                     }
                 }
 
                 // derivative respect transmitted power
                 derivatives[dims] = 1.0;
 
                 // d^2 = (x - xap)^2 + (y - yap)^2
                 // d^n = (d^2)^n/2
 
                 // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
                 // n is the path loss exponent
                 // lambda = c/f, where lambda is wavelength,
                 // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the
                 //               received Gain
                 // Pr = Pte*c^n/((4*pi*f)^n * d^n)
                 // Pr (dBm) = 10*log(k) + 10*log(Pte) - 10*log(d^n) =
                 // 10*log(k) + 10*log(Pte) - 10*log((d^2)^n/2) =
                 // 10*log(k) + 10*log(Pte) - 10*n/2*log(d^2) =
                 // 10*log(k) + 10*log(Pte) - 5*n*log(d^2) =
                 return kdB + transmittedPowerdBm - 5.0 * initialPathLossExponent * Math.log10(sqrDistance);
             }
         });
 
         final var numReadings = readings.size();
         try {
             final var x = new Matrix(numReadings, dimsPlus1);
             final var y = new double[numReadings];
             final var standardDeviations = new double[numReadings];
             for (var i = 0; i < numReadings; i++) {
                 reading = readings.get(i);
                 final var position = reading.getPosition();
 
                 for (var j = 0; j < dims; j++) {
                     x.setElementAt(i, j, position.getInhomogeneousCoordinate(j));
                 }
                 x.setElementAt(i, dims, initialTransmittedPowerdBm);
 
                 standardDeviations[i] = reading.getRssiStandardDeviation() != null ? reading.getRssiStandardDeviation()
                         : DEFAULT_POWER_STANDARD_DEVIATION;
                 y[i] = reading.getRssi();
             }
 
             fitter.setInputData(x, y, standardDeviations);
         } catch (final AlgebraException ignore) {
             // never happens
         }
     }
 
     /**
      * Setups fitter to estimate position and path loss exponent.
      *
      * @throws FittingException if Levenberg-Marquardt fitting fails.
      */
     private void setupFitterPositionAndPathLossExponent() throws FittingException {
         // because all readings must belong to the same radio source, we
         // obtain the frequency of the first radio source on the first reading
         var reading = readings.get(0);
         final var frequency = reading.getSource().getFrequency();
 
         // n = 2.0, is the path loss exponent (which is typically 2.0)
 
         // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
         // lambda = c/f, where lambda is wavelength,
         // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the received Gain
         // Pr = Pte*c^n/((4*pi*f)^n * d^n)
 
 
         // compute k as the constant part of the isotropic received power formula
         // so that: Pr = Pte*k/d^n
         final var k = SPEED_OF_LIGHT / (4.0 * Math.PI * frequency);
         final var kdB = 10.0 * Math.log10(k);
 
         final var dims = getNumberOfDimensions();
         final var dimsPlus1 = dims + 1;
 
         final var initialTransmittedPowerdBm = computeInitialTransmittedPowerdBm();
 
         // for numerical accuracy reasons, a logarithmic version of the previous
         // formula will be used instead
         // Pr (dBm) = 10 * log(Pte * k / d^n) = 10*log(k) + 10*log(Pte) - 10*n*log(d)
 
         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
             @Override
             public int getNumberOfDimensions() {
                 return dimsPlus1;
             }
 
             @Override
             public double[] createInitialParametersArray() {
                 final var initial = new double[dimsPlus1];
 
                 // initial position
                 computeInitialPosition(initial, dims);
 
                 // initial path loss exponent
                 initial[dims] = initialPathLossExponent;
                 return initial;
             }
 
             @Override
             public double evaluate(
                     final int i, final double[] point, final double[] params, final double[] derivatives) {
                 var sqrDistance = 0.0;
                 final var pathLossExponent = params[dims];
                 for (var j = 0; j < dims; j++) {
                     final var diff = params[j] - point[j];
                     sqrDistance += diff * diff;
 
                     // n is mInitialPathLossExponent, which is typically 2.0
                     derivatives[j] = -10.0 * pathLossExponent * diff;
                 }
 
                 // derivatives respect position coordinates are (2D case):
                 // f(x,y,n) = n*kdB -5*n*log((x - xap)^2 + (y - yap)^2)
                 // df/dx = -5*n*2*(x - xap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffX/(ln(10)*sqrDistance)
                 // df/dy = -5*n*2*(y - yap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffY/(ln(10)*sqrDistance)
                 // df/dn = kdB -5*log((x - xap)^2 + (y - yap)^2) = kdB - 5*log(sqrDistance)
                 final var ln10PerSqrDistance = Math.log(10.0) * sqrDistance;
                 if (ln10PerSqrDistance != 0.0) {
                     for (var j = 0; j < dims; j++) {
                         derivatives[j] /= ln10PerSqrDistance;
                     }
                 }
 
                 // derivative respect to path loss exponent
                 final var logSqrDistance = Math.log10(sqrDistance);
                 derivatives[dims] = kdB - 5 * logSqrDistance;
 
                 // d^2 = (x - xap)^2 + (y - yap)^2
                 // d^n = (d^2)^n/2
                 // k = (c/(4*pi*f))^n
 
                 // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
                 // n is the path loss exponent
                 // lambda = c/f, where lambda is wavelength,
                 // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the
                 //               received Gain
                 // Pr = Pte*c^n/((4*pi*f)^n * d^n)
                 // Pr (dBm) = 10*log(k^n) + 10*log(Pte) - 10*log(d^n) =
                 // 10*log(k^n) + 10*log(Pte) - 10*log((d^2)^n/2) =
                 // 10*n*log(k) + 10*log(Pte) - 10*n/2*log(d^2) =
                 // 10*n*log(k) + 10*log(Pte) - 5*n*log(d^2) =
 
                 return pathLossExponent * kdB + initialTransmittedPowerdBm - 5.0 * pathLossExponent * logSqrDistance;
             }
         });
 
         final var numReadings = readings.size();
         try {
             final var x = new Matrix(numReadings, dimsPlus1);
             final var y = new double[numReadings];
             final var standardDeviations = new double[numReadings];
             for (var i = 0; i < numReadings; i++) {
                 reading = readings.get(i);
                 final var position = reading.getPosition();
 
                 for (var j = 0; j < dims; j++) {
                     x.setElementAt(i, j, position.getInhomogeneousCoordinate(j));
                 }
                 x.setElementAt(i, dims, initialPathLossExponent);
 
                 standardDeviations[i] = reading.getRssiStandardDeviation() != null ? reading.getRssiStandardDeviation()
                         : DEFAULT_POWER_STANDARD_DEVIATION;
                 y[i] = reading.getRssi();
             }
 
             fitter.setInputData(x, y, standardDeviations);
         } catch (final AlgebraException ignore) {
             // never happens
         }
     }
 
     /**
      * Setups fitter to estimate transmitted power and path loss exponent.
      *
      * @throws FittingException if Levenberg-Marquardt fitting fails.
      */
     private void setupFitterTransmittedPowerAndPathLossExponent() throws FittingException {
         // because all readings must belong to the same radio source, we
         // obtain the frequency of the first radio source on the first reading
         var reading = readings.get(0);
         final var frequency = reading.getSource().getFrequency();
 
         // n = 2.0, is the path loss exponent
 
         // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
         // lambda = c/f, where lambda is wavelength,
         // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the received Gain
         // Pr = Pte*c^n/((4*pi*f)^n * d^n)
 
 
         // k is defined so that: Pr = Pte * k^n / d^n so that
         // k = (c/(4*pi*f))
         final var k = SPEED_OF_LIGHT / (4.0 * Math.PI * frequency);
         final var kdB = 10.0 * Math.log10(k);
 
         final var initialTransmittedPowerdBm = computeInitialTransmittedPowerdBm();
 
         // for numerical accuracy reasons, a logarithmic version of the previous
         // formula will be used instead
         // Pr (dBm) = 10 * log(Pte * k^n / d^n) = 10*n*log(k) + 10*log(Pte) - 10*n*log(d)
 
         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
             @Override
             public int getNumberOfDimensions() {
                 return 2;
             }
 
             @Override
             public double[] createInitialParametersArray() {
                 final var initial = new double[2];
 
                 // initial transmitted power
                 initial[0] = initialTransmittedPowerdBm;
 
                 // initial path loss exponent
                 initial[1] = initialPathLossExponent;
 
                 return initial;
             }
 
             @Override
             public double evaluate(
                     final int i, final double[] point, final double[] params, final double[] derivatives) {
                 final var sqrDistance = initialPosition.sqrDistanceTo(readings.get(i).getPosition());
                 final var transmittedPowerdBm = params[0];
                 final var pathLossExponent = params[1];
 
                 // derivative respect transmitted power
                 derivatives[0] = 1.0;
 
                 // derivative respect to path loss exponent
                 final var logSqrDistance = Math.log10(sqrDistance);
                 derivatives[1] = kdB - 5 * logSqrDistance;
 
                 // d^2 = (x - xap)^2 + (y - yap)^2
                 // d^n = (d^2)^n/2
                 // k = (c/(4*pi*f))^n
 
                 // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
                 // n is the path loss exponent
                 // lambda = c/f, where lambda is wavelength,
                 // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the
                 //                  received Gain
                 // Pr = Pte*c^n/((4*pi*f)^n * d^n)
                 // Pr (dBm) = 10*log(k^n) + 10*log(Pte) - 10*log(d^n) =
                 // 10*log(k^n) + 10*log(Pte) - 10*log((d^2)^n/2) =
                 // 10*n*log(k) + 10*log(Pte) - 10*n/2*log(d^2) =
                 // 10*n*log(k) + 10*log(Pte) - 5*n*log(d^2) =
 
                 return pathLossExponent * kdB + transmittedPowerdBm - 5.0 * pathLossExponent * logSqrDistance;
             }
         });
 
         final var numReadings = readings.size();
         try {
             final var x = new Matrix(numReadings, 2);
             final var y = new double[numReadings];
             final var standardDeviations = new double[numReadings];
             for (var i = 0; i < numReadings; i++) {
                 reading = readings.get(i);
 
                 x.setElementAt(i, 0, initialTransmittedPowerdBm);
                 x.setElementAt(i, 1, initialPathLossExponent);
 
                 standardDeviations[i] = reading.getRssiStandardDeviation() != null ? reading.getRssiStandardDeviation()
                         : DEFAULT_POWER_STANDARD_DEVIATION;
                 y[i] = reading.getRssi();
             }
 
             fitter.setInputData(x, y, standardDeviations);
         } catch (AlgebraException ignore) {
             // never happens
         }
     }
 
     /**
      * Setups fitter to estimate transmitted power, position and path
      * loss exponent.
      *
      * @throws FittingException if Levenberg-Marquardt fitting fails.
      */
     private void setupFitterPositionTransmittedPowerAndPathLossExponent() throws FittingException {
         // because all readings must belong to the same radio source, we
         // obtain the frequency of the first radio source on the first reading
         var reading = readings.get(0);
         final var frequency = reading.getSource().getFrequency();
 
         // n = 2.0, is the path loss exponent
 
         // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
         // lambda = c/f, where lambda is wavelength,
         // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the received Gain
         // Pr = Pte*c^n/((4*pi*f)^n * d^n)
 
         // k is defined so that: Pr = Pte * k^n / d^n so that
         // k = (c/(4*pi*f))
         final var k = SPEED_OF_LIGHT / (4.0 * Math.PI * frequency);
         final var kdB = 10.0 * Math.log10(k);
 
         final var dims = getNumberOfDimensions();
         final var dimsPlus1 = dims + 1;
         final var dimsPlus2 = dims + 2;
 
         final var initialTransmittedPowerdBm = computeInitialTransmittedPowerdBm();
 
         // for numerical accuracy reasons, a logarithmic version of the previous
         // formula will be used instead
         // Pr (dBm) = 10 * log(Pte * k^n / d^n) = 10*n*log(k) + 10*log(Pte) - 10*n*log(d)
 
         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
 
             @Override
             public int getNumberOfDimensions() {
                 return dimsPlus2;
             }
 
             @Override
             public double[] createInitialParametersArray() {
                 final var initial = new double[dimsPlus2];
 
                 // initial position
                 computeInitialPosition(initial, dims);
 
                 // initial transmitted power
                 initial[dims] = initialTransmittedPowerdBm;
 
                 // initial path loss exponent
                 initial[dimsPlus1] = initialPathLossExponent;
 
                 return initial;
             }
 
             @Override
             public double evaluate(
                     final int i, final double[] point, final double[] params, final double[] derivatives) {
                 var sqrDistance = 0.0;
                 final var pathLossExponent = params[dimsPlus1];
                 for (var j = 0; j < dims; j++) {
                     final var diff = params[j] - point[j];
                     sqrDistance += diff * diff;
 
                     // n is mInitialPathLossExponent, which is typically 2.0
                     derivatives[j] = -10.0 * pathLossExponent * diff;
                 }
 
                 final var transmittedPowerdBm = params[dims];
 
                 // derivatives respect position coordinates are (2D case):
                 // f(x,y,n) = n*kdB -5*n*log((x - xap)^2 + (y - yap)^2)
                 // df/dx = -5*n*2*(x - xap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffX/(ln(10)*sqrDistance)
                 // df/dy = -5*n*2*(y - yap)/(ln(10)*((x - xap)^2 + (y - yap)^2)) = -10*n*diffY/(ln(10)*sqrDistance)
                 // df/dn = kdB -5*log((x - xap)^2 + (y - yap)^2) = kdB - 5*log(sqrDistance)
                 final var ln10PerSqrDistance = Math.log(10.0) * sqrDistance;
                 if (ln10PerSqrDistance != 0.0) {
                     for (var j = 0; j < dims; j++) {
                         derivatives[j] /= ln10PerSqrDistance;
                     }
                 }
 
                 // derivative respect transmitted power
                 derivatives[dims] = 1.0;
 
                 // derivative respect to path loss exponent
                 final var logSqrDistance = Math.log10(sqrDistance);
                 derivatives[dimsPlus1] = kdB - 5 * logSqrDistance;
 
                 // d^2 = (x - xap)^2 + (y - yap)^2
                 // d^n = (d^2)^n/2
                 // k = (c/(4*pi*f))^n
 
                 // Pr = Pt*Gt*Gr*lambda^n/(4*pi*d)^n,    where Pr is the received power
                 // n is the path loss exponent
                 // lambda = c/f, where lambda is wavelength,
                 // Pte = Pt*Gt*Gr, is the equivalent transmitted power, Gt is the transmitted Gain and Gr is the
                 //           received Gain
                 // Pr = Pte*c^n/((4*pi*f)^n * d^n)
                 // Pr (dBm) = 10*log(k^n) + 10*log(Pte) - 10*log(d^n) =
                 // 10*log(k^n) + 10*log(Pte) - 10*log((d^2)^n/2) =
                 // 10*n*log(k) + 10*log(Pte) - 10*n/2*log(d^2) =
                 // 10*n*log(k) + 10*log(Pte) - 5*n*log(d^2) =
 
                 return pathLossExponent * kdB + transmittedPowerdBm - 5.0 * pathLossExponent * logSqrDistance;
             }
         });
 
         final var numReadings = readings.size();
         try {
             final var x = new Matrix(numReadings, dimsPlus2);
             final var y = new double[numReadings];
             final var standardDeviations = new double[numReadings];
             for (var i = 0; i < numReadings; i++) {
                 reading = readings.get(i);
                 final var position = reading.getPosition();
 
                 for (var j = 0; j < dims; j++) {
                     x.setElementAt(i, j, position.getInhomogeneousCoordinate(j));
                 }
                 x.setElementAt(i, dims, initialTransmittedPowerdBm);
                 x.setElementAt(i, dimsPlus1, initialPathLossExponent);
 
                 standardDeviations[i] = reading.getRssiStandardDeviation() != null ? reading.getRssiStandardDeviation()
                         : DEFAULT_POWER_STANDARD_DEVIATION;
                 y[i] = reading.getRssi();
             }
 
             fitter.setInputData(x, y, standardDeviations);
         } catch (final AlgebraException ignore) {
             // never happens
         }
     }
 
     /**
      * Computes initial transmitted power expressed in dBm's.
      * If no initial transmitted power is provided, the average of all measures
      * is used, otherwise provided value is used.
      *
      * @return initial transmitted power.
      */
     private double computeInitialTransmittedPowerdBm() {
         if (initialTransmittedPowerdBm == null) {
             // compute average transmitted power (in mW)
             final var num = readings.size();
             var result = 0.0;
             for (final var reading : readings) {
                 final var rssi = reading.getRssi();
                 result += rssi / num;
             }
             return result;
         } else {
             // convert initial value
             return initialTransmittedPowerdBm;
         }
     }
 
     /**
      * Computes initial position.
      *
      * @param result array where result will be stored.
      * @param dims   number of dimensions of position coordinates (either 2 or 3).
      */
     private void computeInitialPosition(final double[] result, final int dims) {
         if (initialPosition == null) {
             final var num = readings.size();
 
             // compute average centroid of fingerprint positions
             for (final var reading : readings) {
                 final var position = reading.getPosition();
                 for (var i = 0; i < dims; i++) {
                     result[i] += position.getInhomogeneousCoordinate(i) / num;
                 }
             }
         } else {
             // copy initial position
             for (var i = 0; i < dims; i++) {
                 result[i] = initialPosition.getInhomogeneousCoordinate(i);
             }
         }
     }
 }