Package newExamples

Source Code of newExamples.ModelFunction

package newExamples;

import gaia.cu1.params.GaiaParam;
import gaia.cu1.tools.exception.GaiaException;
import gaia.cu1.tools.numeric.algebra.GVector2d;

import gaiasimu.SimuException;
import gaiasimu.gaia.spacecraft.Transit;
import gaiasimu.universe.source.stellar.ExoPlanet;
import gaiasimu.universe.source.stellar.ExoPlanetPhysicalParameters;
import gaiasimu.universe.source.stellar.OrbitalParams;
import gaiasimu.universe.source.stellar.ComponentAstrometry;
import gaiasimu.universe.source.stellar.StarSystemSimuException;
import gaiasimu.universe.source.stellar.StellarAstrometry;

import java.io.IOException;
import java.util.ArrayList;
import java.util.TreeMap;

import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
import org.apache.commons.math.analysis.MultivariateMatrixFunction;


/**
* Model function used to fit the planets
*
* @author ebopp
*/
public class ModelFunction implements DifferentiableMultivariateVectorialFunction, Logging {

  /** Underlying model */
  private Model model;

  /** Setup of model parameters */
  private TreeMap<FitQuantity, FitVariable> paramSetup;

  /** The use AC data boolean switch. */
  private boolean useAc;
 
  /** The dimension of the vector space of the function values */
  private int valueDimension;
 
  /** Array of past parameters */
  private ArrayList<TreeMap<FitQuantity, Double> > pastParams;
 
  /** Array of past evaluation results */
  private ArrayList<double[]> pastResults;
 
  /** List of reference transits */
  private ArrayList<Transit> refTransits;


  /**
   * Constructor
   *
   * @param model      Fitting model of simulated system
   * @param useAc      Include across scan?
   * @param refTransits  List of reference transits to be fitted
   *
   * @throws SimuException
   * @throws GaiaException
   * @throws IOException
   */
  public ModelFunction(Model model, boolean useAc, ArrayList<Transit> refTransits)
  throws SimuException, GaiaException, IOException {

    // Initiate private variables
    this.model      = model;
    this.useAc      = useAc;
    this.refTransits  = refTransits;
    this.valueDimension  = refTransits.size() * (useAc? 2: 1);
   
    // Initialize stuff
    paramSetup = new TreeMap<FitQuantity, FitVariable>();
    pastParams = new ArrayList<TreeMap<FitQuantity, Double> >();
    pastResults = new ArrayList<double[]>();
   
    // Default settings for fitting parameters
    setVariable(FitQuantity.INCLINATION, false, Math.PI/4, 1e-5);
    setVariable(FitQuantity.NODEANGLE, true, 0.0, 1e-5);
    setVariable(FitQuantity.TIMEPERI, true, 0.0, 1e-2);
    setVariable(FitQuantity.OMEGA2, true, 0.0, 1e-5);
    setVariable(FitQuantity.ALPHA_OFFSET, true, 0.0, 1e-5);
    setVariable(FitQuantity.DELTA_OFFSET, true, 0.0, 1e-5);
    setVariable(FitQuantity.PARALLAX, true, model.getAstrometry().getVarpi(), 1e-5);
    setVariable(FitQuantity.MU_ALPHA, true, model.getAstrometry().getMuAlphaStar(), 1e-5);
    setVariable(FitQuantity.MU_DELTA, true, model.getAstrometry().getMuDelta(), 1e-5);
   
  }
 
 
  /**
   * Sets a fit variable.
   *
   * @param param      The quantity associated with this variable
   * @param fixed      Whether the variable should be fixed
   * @param defaultValue  The default value [unit depends on quantity]
   * @param diffH      Differential increment used for derivates w.r.t. this variable
   */
  public void setVariable ( FitQuantity param, boolean fixed, double defaultValue, double diffH ) {
   
    FitVariable var = new FitVariable ( fixed, defaultValue, diffH );
    paramSetup.put( param, var );
   
  }
 
 
  /**
   * Gets the variable associated to a certain physical quantity.
   *
   * @param param  The associated physical quantity
   * @return  The fit variable
   */
  public FitVariable getVariable ( FitQuantity param ) {
   
    return paramSetup.get ( param );
   
  }
 

  /**
   * Counts free parameters, effectively degrees of freedom.
   *
   * @return  Number of free parameters of this model function.
   */
  public int getNVarParams () {
   
    // Count free variables
    int varParams = 0;
    for ( FitQuantity qua: FitQuantity.values() )
      if ( ! paramSetup.get(qua).isFixed() )
        varParams++;
   
    // Return
    return varParams;
   
  }
 
 
  /**
   * Gets an array of initial values for the fitting parameters.
   *
   * @return  Array of initial values
   */
  public double[] getInitialValues () {
   
    // Build up array of initial values
    double[] initials = new double[getNVarParams()];
    int i = 0;
    for ( FitQuantity qua: FitQuantity.values() )
      if ( ! paramSetup.get(qua).isFixed() ) {
        initials[i] = paramSetup.get(qua).getDefaultValue();
        i++;
      }
   
    // Return
    return initials;
   
  }
 

  /**
   * Overloaded jacobian function, realized with two evaluations separated by small angle increment
   *
   * @see    org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction#jacobian()
   * @return  Jacobian matrix
   */
  @Override
  public MultivariateMatrixFunction jacobian() {
   
    return new MultivariateMatrixFunction() {
     
      public double[][] value ( double[] params ) {

        // Set up jacobian and zero offset values
        double[][] jacobianMatrix = new double [valueDimension][params.length];
        TreeMap<FitQuantity, Double> zeroMap = getMap(params);
        double[] zeroValue = evaluateFunction ( zeroMap );
       
        // Loop through fit variables
        int i = 0;
        for ( FitQuantity qua: FitQuantity.values() )
          if ( ! paramSetup.get(qua).isFixed() ) {
            double[] offset = evaluateWithOffset ( zeroMap, qua );
            double dx = paramSetup.get(qua).getDiffH();
            for ( int j = 0; j < valueDimension; j++ )
              jacobianMatrix[j][i] = (offset[j] - zeroValue[j]) / dx;
            i++;
          }

        // Return
        return jacobianMatrix;

      }

    };

  }


  /**
   * Overridden value method of multivariate vectorial function.
   *
   * @see org.apache.commons.math.analysis.MultivariateVectorialFunction#value(double[])
   *
   * @param params  The array containing supplied parameters
   * @return  The array of values
   */
  @Override
  public double[] value(double[] params)
  throws FunctionEvaluationException, IllegalArgumentException {
   
    // Map entries
    TreeMap<FitQuantity, Double> inputMap = getMap(params);
    if ( inputMap == null )
      return new double[valueDimension];

    // Call evaluation method
    return evaluateFunction(inputMap);

  }
 
 
  /**
   * Converts array of parameters to tree map.
   *
   * @param params  Parameter array
   * @return  Tree map
   */
  private TreeMap<FitQuantity, Double> getMap(double[] params) {

    // Check for correct array size
    if ( params.length != getNVarParams() )
      return null;
   
    // Map entries
    TreeMap<FitQuantity, Double> inputMap = new TreeMap<FitQuantity, Double>();
    int i = 0;
    for(FitQuantity qua: FitQuantity.values())
      if(!paramSetup.get(qua).isFixed()) {
        inputMap.put(qua, params[i]);
        i++;
      }
   
    // Return
    return inputMap;
   
  }


  /**
   * Internal function evaluation routine used for value and jacobian
   *
   * @param params  map of supplied quantities
   *
   * INCLINATION  Inclination of planetary orbit [rad]
   * NODEANGLE  Longitude of ascending node; node angle [rad]
   * ALPHA_OFFSET  Right ascension offset of star [mas]
   * DELTA_OFFSET  Declination offset of star [mas]
   * PARALLAX    Parallax [mas]
   * MU_ALPHA    Proper motion (alpha) [mas/yr]
   * MU_DELTA    Proper motion (delta) [mas/yr]
   *
   * @return    Value of function
   */
  protected double[] evaluateFunction(TreeMap<FitQuantity, Double> paramMap) {

    // In case of past evaluation return earlier result
    int index = 0;
    for ( TreeMap<FitQuantity, Double> pastParam: pastParams ) {
      boolean matches = true;
      for(FitQuantity qua: FitQuantity.values())
        if(paramMap.containsKey(qua) && (Math.abs(pastParam.get(qua) - paramMap.get(qua)) >= paramSetup.get(qua).getDiffH() * 0.5))
          matches = false;
      if(matches)
        return pastResults.get(index);
      index++;
    }
   
    // Create complete map including default values for fixed parameters not included in the map
    TreeMap<FitQuantity, Double> completeMap = new TreeMap<FitQuantity, Double>();
    for ( FitQuantity qua: FitQuantity.values() )
      if ( paramMap.containsKey ( qua ) )
        completeMap.put ( qua, paramMap.get(qua) );
      else
        completeMap.put ( qua, paramSetup.get(qua).getDefaultValue() );
   
    // Set up variables
    double newIncl = completeMap.get ( FitQuantity.INCLINATION );
    double nodeAngle = completeMap.get ( FitQuantity.NODEANGLE );
    double timePeri = completeMap.get ( FitQuantity.TIMEPERI );
    double omega2 = completeMap.get ( FitQuantity.OMEGA2 );
    double alphaOff = completeMap.get ( FitQuantity.ALPHA_OFFSET );
    double deltaOff = completeMap.get ( FitQuantity.DELTA_OFFSET );
    double parallax = completeMap.get ( FitQuantity.PARALLAX );
    double muAlphaStar = completeMap.get ( FitQuantity.MU_ALPHA );
    double muDelta = completeMap.get ( FitQuantity.MU_DELTA );
   
    // Debug output of parameters
    logger.debug ( String.format ( "%20.12e %20.12e %20.12e %20.12e %20.12e %20.12e %20.12e", newIncl, nodeAngle, alphaOff, deltaOff, parallax, muAlphaStar, muDelta ) );

    // Calculate right ascension and declination
    double alpha = (model.getRefAstro().getAlpha() + GaiaParam.Nature.MILLIARCSECOND_RADIAN * alphaOff) * GaiaParam.Nature.RADIAN_DEGREE;
    double delta = (model.getRefAstro().getDelta() + GaiaParam.Nature.MILLIARCSECOND_RADIAN * deltaOff) * GaiaParam.Nature.RADIAN_DEGREE;
   
    try {

      // Some shortcuts
      ExoPlanet planet = model.getPlanets().get(0);
      ExoPlanetPhysicalParameters phys = (ExoPlanetPhysicalParameters) planet.getPhysParam();
      OrbitalParams orb = ((ComponentAstrometry) planet.getAstrometry()).getOrbitalParams();
      StellarAstrometry astro = model.getAstrometry();
     
      // Change astrometry
      model.setAstrometry(alpha, delta, parallax, muAlphaStar, muDelta, astro.getRadialVelocity());
     
      // Modify planet
      double mSinI = phys.getMass() * Math.sin(orb.inclination);
      phys.setMass(mSinI / Math.sin(newIncl));
      orb.inclination = newIncl;
      orb.nodeangle = nodeAngle;
      orb.omega2 = omega2;
      orb.timeperiastron = timePeri;
     
      // Finalize the new model's system
      model.finalizeSystem();
     
      // Obtain field angles
      double[] vector;

      vector = new double[valueDimension];
      for ( int i=0; i < refTransits.size(); i++ ) {
       
        GVector2d fieldAngles = model.getFieldAngles(refTransits.get(i));
       
        if ( useAc ) {
          vector[2*i]   = fieldAngles.getX();
          vector[2*i+1] = fieldAngles.getY();
        } else {
          vector[i]   = fieldAngles.getX();
        }
       
      }
     
      // Save evaluation
      pastParams.add(paramMap);
      pastResults.add(vector);
     
      // Return
      return vector;

    } catch (SimuException e) {
      logger.error("SimuException occured while evaluating function.");
    } catch (StarSystemSimuException e) {
      logger.error("StarSystemSimuException occured while evaluating function.");
    }

    // Default zero in case of errors
    double[] defResult = new double[valueDimension];
    for ( int i = 0; i < defResult.length; i++ )
      defResult[i] = 0.0;
    return defResult;

  }
 
 
  /**
   * Returns the value of the model function offset by the differential increment of the given quantity
   *
   * @param params  The map of supplied quantities
   * @param diffQua  The offset quantity
   *
   * INCLINATION  Inclination of planetary orbit [rad]
   * ALPHA_OFFSET  Right ascension offset of star [mas]
   * DELTA_OFFSET  Declination offset of star [mas]
   * PARALLAX    Parallax [mas]
   * MU_ALPHA    Proper motion (alpha) [mas/yr]
   * MU_DELTA    Proper motion (delta) [mas/yr]
   *
   * @return    Value of partial differential w.r.t.
   */
  @SuppressWarnings("unchecked")
  protected double[] evaluateWithOffset ( TreeMap<FitQuantity, Double> paramMap, FitQuantity diffQua ) {
   
    // Set up differentiation variable
    FitVariable diffVar = paramSetup.get ( diffQua );
    double dx = diffVar.getDiffH();
   
    // Prepare inputs
    TreeMap<FitQuantity, Double> inputMap = (TreeMap<FitQuantity, Double>) paramMap.clone();
    double xVal = paramMap.get(diffQua);
    inputMap.put ( diffQua, xVal + dx );
   
    // Call evaluation function twice
    double[] result = evaluateFunction ( inputMap );
   
    // Return
    return result;
   
  }


}
TOP

Related Classes of newExamples.ModelFunction

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.