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;
}
}