// // FlatField.java // /* VisAD system for interactive analysis and visualization of numerical data. Copyright (C) 1996 - 2017 Bill Hibbard, Curtis Rueden, Tom Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and Tommy Jasmin. This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA */ package visad; import java.rmi.RemoteException; import java.util.Arrays; import java.util.Random; /** FlatField is the VisAD class for finite samplings of functions whose range type and range coordinate systems are simple enough to allow efficient representation. The DomainSet, DomainCoordinateSystem, RangeSet, RangeCoordinateSystem and RangeCoordinateSystems variables of FlatField are immutable.and will be of * length 1. * @throws TypeException The type of the range is neither RealType nor * RealTupleType. */ public CoordinateSystem[] getRangeCoordinateSystem() throws TypeException { MathType RangeType = ((FunctionType) Type).getRange(); if (!((FunctionType) Type).getReal()) { throw new TypeException("FlatField.getRangeCoordinateSystem: " + "Range is not Real, need DefaultCoordinateSystem index"); } CoordinateSystem[] cs = {RangeCoordinateSystem}; return cs; } /** * Returns the sampling set of each flat component. * @return The sampling set of each component in the flat range. */ public Set[] getRangeSets() { Set[] sets = new Set[RangeSet.length]; System.arraycopy(RangeSet, 0, sets, 0, sets.length); return sets; } /** * Returns the CoordinateSystem of a component of the range. The MathType * of the range shall be a TupleType. * * @param i The index of the component. The value shall be * greater than or equal to zero and less that the * number of components in the TupleType of the * range. * @return The CoordinateSystem of theA FlatField range type may be either a RealType (for a function with range = R), a RealTupleType (for a function with range = R^n for n > 0), or a TupleType of RealType-s and RealTupleType-s..
VisAD avoids invoking methods once per datum through the use of FlatField's. These are logically Field's of Tuple's of RealType's and RealTupleType's. Internally FlatField's are stored as arrays of numerical values, rather than the arrays of data objects stored in Field's. Many of the methods in the FlatField class and in other classes (e.g., CoordinateTransform, Set, Unit) process data in the form double[Dimension][Length] where Length is the number of samples in a Field and Dimension is the number of Tuple elements in the Field range values. Note that the order of the Length and Dimension indices are reversed as array indices. This allows efficient processing of long columns of Field value components. For example, if Latitude is one component of Field values, then any computation involving Latitude can be applied in a tight loop to all Latitude's in the Field.
FlatField's support range types more general than RealTuple's. To understand the motive, consider a set of observations that include Latitude, Longitude, Altitude, Pressure, Temperature, etc. We can organize these as a Field whose range values have the Tuple type:
(Latitude, Longitude, Altitude, Pressure, Temperature, ...)However, in order to declare that (Latitude, Longitude, Altitude) is a coordinate system with coordinate transform functions to other spatial coordinate systems, we need to organize:
(Latitude, Longitude, Altitude)as a RealTupleType. Hence the range type of the Field of observations must be:
((Latitude, Longitude, Altitude), Pressure, Temperature, ...)which is not a RealTupleType (since one of its components is a RealTupleType). In order to process such data efficiently, FlatField's must support range types that are Tuple's of RealType's and RealTupleType's.*/ public class FlatField extends FieldImpl implements FlatFieldIface { protected int TupleDimension; // dimension of Type.getFlatRange() private Set RangeSet[]; // one 1-D Set per range components private int RangeMode[]; // DOUBLE, FLOAT, INT, SHORT or BYTE // coordinate system of the function range R^n protected CoordinateSystem RangeCoordinateSystem; // used if Type.Real protected CoordinateSystem[] RangeCoordinateSystems; // used for Flat, not for Real private boolean MissingFlag; protected Unit[] RangeUnits; // RangeErrors, like range values, are not immutable private ErrorEstimate[] RangeErrors; // for tuple memeber i, only one of these is not null, depending on RangeSet[i] private double[][] DoubleRange; private float[][] FloatRange; private long[][] LongRange; // not currently used because array // indices are int's private int[][] IntRange; private short[][] ShortRange; private byte[][] ByteRange; private static final int MISSING1 = Byte.MIN_VALUE; // least byte private static final int MISSING2 = Short.MIN_VALUE; // least short private static final int MISSING4 = Integer.MIN_VALUE; // least int // private static final int MISSING8 = Long.MIN_VALUE; // least long private static final int DOUBLE = 1; private static final int FLOAT = 2; private static final int LONG = 3; // not currently used because array // indices are int's private static final int INT = 4; private static final int SHORT = 5; private static final int BYTE = 6; /** * Constructs a FlatField from a function type. The domain Set is the * default Set of the function domain. * @param type The type of the function. * @throws VisADException Couldn't create necessary VisAD object. */ public FlatField(FunctionType type) throws VisADException { this(type, type.getDomain().getDefaultSet(), null, null, null, null); } /** * Constructs a FlatField from a function type and a (non-default) domain * Set. * @param type The type of the function. * @param domain_set The sampling set of the domain. May be
* null, in which case the sampling set is * the default Set of the function domain. * @throws VisADException Couldn't create necessary VisAD object. */ public FlatField(FunctionType type, Set domain_set) throws VisADException { this(type, domain_set, null, null, null, null); } /** * Constructs a FlatField from a function type, a sampling set of the domain, * a coordinate system for the range, sampling sets for the range components, * and units for the range components. * @param type The type of the function. * @param domain_set The sampling set of the domain. May be *null, in which case the * sampling set is set to the default sampling * set of the function domain (i.e.* type.getDomain().getDefaultSet()). * @param range_coord_sys Optional coordinate system for the range. * May benull. If non-null *, then the range of the function shall * be a RealTuple (i.e.type.getRange() * instanceof RealTupleshall be * true) and the reference coordinate system * of the range coordinate system shall be * the same as the reference coordiinate * system of the function range (i.e. *range_coord_sys().getReference().equals( * ((RealTuple)type.getRange()).getReference()) *shall be true). * @param range_sets The sampling sets of the (flat) range * components. May benull, in which * case the default sampling sets of the range * component RealType-s are used (i.e.* ((RealType)type.getFlatRange().getComponent(i)) * .getDefaultSet()for alli* in the flat range). If non-null, * then theith flat range * component values are stored in bytes * ifrange_sets[i].getLength() < * 256, stored in shorts if* range_sets[i].getLength() < 65536, etc. * @param units The units of the (flat) range components. May * benull, in which case the default * units of the flat range RealType-s are used. * @throws VisADException Couldn't create necessary VisAD object. */ public FlatField(FunctionType type, Set domain_set, CoordinateSystem range_coord_sys, Set[] range_sets, Unit[] units) throws VisADException { this(type, domain_set, range_coord_sys, null, range_sets, units); } /** * Constructs a FlatField from a function type, a sampling set of the domain, * coordinate systems for the range components, sampling sets for the range * components, and units for the range components. * @param type The type of the function. * @param domain_set The sampling set of the domain. May be null, in * which case the sampling set is set to the * default sampling set of the function domain * (i.e.* type.getDomain().getDefaultSet()). * @param range_coord_syses Optional coordinate systems for the range * components. May benull. * If theith component * of the range is a RealTuple, then *range_coord_syses[i]may be * non-nulland, if so, is the * coordinate system for that component and shall * have the same reference coordinate system as * the corresponding component in the range of the * function type. * @param range_sets The sampling sets of the (flat) range * components. May benull, in which * case the default sampling sets of the range * component RealType-s are used (i.e.* ((RealType)type.getFlatRange().getComponent(i)) * .getDefaultSet()for alli* in the flat range). If non-null, * then theith flat range * component values are stored in bytes * ifrange_sets[i].getLength() < * 256, stored in shorts if* range_sets[i].getLength() < 65536, etc. * @param units The units of the (flat) range components. May * benull, in which case the default * units of the flat range RealType-s are used. * @throws VisADException Couldn't create necessary VisAD object. */ public FlatField(FunctionType type, Set domain_set, CoordinateSystem[] range_coord_syses, Set[] range_sets, Unit[] units) throws VisADException { this(type, domain_set, null, range_coord_syses, range_sets, units); } /** * Constructs a FlatField from a function type, a sampling set of the domain, * a coordinate system for the range, coordinate systems for the range * components, sampling sets for the range components, and units for the range * components. This is the most general constructor. * @param type The type of the function. * @param domain_set The sampling set of the domain. May be null, in * which case the sampling set is set to the * default sampling set of the function domain * (i.e.* type.getDomain().getDefaultSet()). * @param range_coord_sys Optional coordinate system for the range. * May benull. If non-null *, then the range of the function shall * be a RealTuple (i.e.type.getRange() * instanceof RealTupleshall be * true) and the reference coordinate system * of the range coordinate system shall be * the same as the reference coordiinate * system of the function range (i.e. *range_coord_sys().getReference().equals( * ((RealTuple)type.getRange()).getReference()) *shall be true). * @param range_coord_syses Optional coordinate systems for the range * components. May benull. * If theith component * of the range is a RealTuple, then *range_coord_syses[i]may be * non-nulland, if so, is the * coordinate system for that component and shall * have the same reference coordinate system as * the corresponding component in the range of the * function type. * @param range_sets The sampling sets of the (flat) range * components. May benull, in which * case the default sampling sets of the range * component RealType-s are used (i.e.* ((RealType)type.getFlatRange().getComponent(i)) * .getDefaultSet()for alli* in the flat range). If non-null, * then theith flat range * component values are stored in bytes * ifrange_sets[i].getLength() < * 256, stored in shorts if* range_sets[i].getLength() < 65536, etc. * @param units The units of the (flat) range components. May * benull, in which case the default * units of the flat range RealType-s are used. * @throws SetException ifrange_setsis non-null *and either one of its elements is *nullor an element's dimension is * not one. * @throws VisADException Couldn't create necessary VisAD object. */ public FlatField(FunctionType type, Set domain_set, CoordinateSystem range_coord_sys, CoordinateSystem[] range_coord_syses, Set[] range_sets, Unit[] units) throws VisADException { //Call FieldImpl.ctor passing false to tell it not to create the Range array super(type, domain_set, false); pr ("ctor"); if (!type.getFlat()) { throw new FieldException("FlatField: FunctionType must be Flat"); } MathType RangeType = type.getRange(); RealTupleType FlatRange = type.getFlatRange(); TupleDimension = FlatRange.getDimension(); DoubleRange = new double[TupleDimension][]; // set RangeSet RangeSet = new Set[TupleDimension]; RangeMode = new int[TupleDimension]; if (range_sets == null) { // set the default range sampling; if no default, use double for (int i=0; inull i-th * component. Won't benulland will * be of length 1. The single element might be *null. * @throws TypeException The type of the range is either RealType or * RealTupleType. */ public CoordinateSystem[] getRangeCoordinateSystem(int i) throws TypeException { if (((FunctionType) Type).getReal()) { throw new TypeException("FlatField.getRangeCoordinateSystem: " + "Range is Real, cannot specify CoordinateSystem index"); } CoordinateSystem[] cs = {RangeCoordinateSystems[i]}; return cs; } /** return array of Units associated with each RealType component of range; these may differ from default Units of range RealTypes, but must be convertable; the second index has length = 1 (since all samples have the same Units) */ public Unit[][] getRangeUnits() { Unit[][] units = new Unit[RangeUnits.length][1]; for (int i=0; ierrors[i] is the {@link ErrorEstimate} for the * i-th {@link RealType} component. * * @param errors The error estimates for the range values. * @throws FieldException iferrorsis non-nulland *errors.length != getRangeDimension(). */ public void setRangeErrors(ErrorEstimate[] errors) throws FieldException { synchronized (RangeErrors) { if (errors == null) { for (int i=0; igetLength()) { throw new FieldException("setSamples: bad array length"); } } boolean copy = false; int length = range[0].length; //-- packValues starts here: // NOTE INVERTED ORDER OF range ARRAY INDICES !!! int[] index; synchronized (DoubleRange) { for (int i=0; i getLength() ) { throw new FieldException("setSamples: indices array too long"); } for (int i=0; i [i][j] is * the jth sample value of the *ith component of the range. * @throws VisADException if a VisAD object couldn't be created. */ public float[][] getFloats() throws VisADException { return getFloats(true); } /** * Returns the range values in their default units as floats. * * @param copy Whether or not the returned array might be the * actual range array. Iftrue, then * the returned array will not be the actual range * array. * @return The range values in their default units * as determined by the {@link MathType} of * the range. Element[i][j]is * thejth sample value of the *ith component of the range. * @throws VisADException if a VisAD object couldn't be created. */ public float[][] getFloats (boolean copy) throws VisADException { pr ("getFloats(" + copy + ")"); float[][] values = unpackFloats(copy); Unit[] units_out = ((FunctionType) Type).getFlatRange().getDefaultUnits(); //Only converty the data if the output units not equal to the input units if (!Arrays.equals (units_out, RangeUnits)) { values = Unit.convertTuple(values, RangeUnits, units_out); } return values; } /** * Returns the range values in their default units as doubles. * * @return The range values in their default units * as determined by the {@link MathType} of * the range. Element[i][j]is * thejth sample value of the *ith component of the range. * @throws VisADException if a VisAD object couldn't be created. */ public double[][] getValues() throws VisADException { return getValues(true); } /** * Returns the range values in their default units as doubles. * * @param copy Whether or not the returned array might be the * actual range array. Iftrue, then * the returned array will not be the actual range * array. * @return The range values in their default units * as determined by the {@link MathType} of * the range. Element[i][j]is * thejth sample value of the *ith component of the range. * @throws VisADException if a VisAD object couldn't be created. */ public double[][] getValues(boolean copy) throws VisADException { pr ("getValues(" + copy + ")"); double[][] values = unpackValues(copy); Unit[] units_out = ((FunctionType) Type).getFlatRange().getDefaultUnits(); double[][]result = Unit.convertTuple(values, RangeUnits, units_out); return result; } /** * Get String values for Text components * @return null (there are none for FlatFields) * @throws VisADException but doesn't happen * @throws RemoteException but doesn't happen */ public String[][] getStringValues() throws VisADException, RemoteException { return null; } /** * Get values for 'Flat' components in default range Unit-s. * @return The range values in their default units * as determined by the {@link MathType} of * the range. Element[i]is * the value of theith component * of the flattened range. * @throws VisADException if a VisAD object couldn't be created. */ /*- TDR June 1998 */ public double[] getValues( int s_index ) throws VisADException { double[] values = new double[TupleDimension]; values = unpackValues( s_index ); double[][] n_values = new double[TupleDimension][1]; for ( int ii = 0; ii < TupleDimension; ii++ ) { n_values[ii][0] = values[ii]; } Unit[] units_out = ((FunctionType) Type).getFlatRange().getDefaultUnits(); double[][] r_values = Unit.convertTuple(n_values, RangeUnits, units_out); for ( int ii = 0; ii < values.length; ii++ ) { values[ii] = r_values[ii][0]; } return values; } /** * Get default range Unit-s for 'Flat' components. * @return array of the default Units for each of the RealTypes * components in the flattened range. */ public Unit[] getDefaultRangeUnits() { return ((FunctionType) Type).getFlatRange().getDefaultUnits(); } /** * Get the range value at the index-th sample. * FlatField does not override evaluate, but the correctness * of FlatField.evaluate depends on overriding getSample * @return Data object (Real, RealTuple, or Tuple) corresponding to * the range at the index-th sample. * @throws VisADException problem getting data * @throws RemoteException problem getting data from remote object */ public Data getSample(int index) throws VisADException, RemoteException { int[] inds; if (isMissing() || index < 0 || index >= getLength()) { return ((FunctionType) Type).getRange().missingData(); } double[][] range = new double[TupleDimension][1]; double[][] range1; synchronized (DoubleRange) { for (int i=0; i= getLength()) return; // disect range into doubles double[] vals = new double[TupleDimension]; // holder for errors of transformed values; ErrorEstimate[] errors_out = new ErrorEstimate[TupleDimension]; if (range == null) { for (int j=0; j values[i] will have been replaced * if necessary. * @param units The units for the values and error estimates. * It's length shall be values.length. * On output,units[i]will have been * replaced with the absolute form of the input * unit if necessary. * @param errors The error estimates. It's length shall be *values.length. On output, *errors[i]will have been replaced * if necessary. * @throws UnitException Unit conversion error. */ protected static void makeRational(double[][] values, Unit[] units, ErrorEstimate[] errors) throws UnitException { for (int j=0; jn_comps ) { throw new VisADException("extract: component index too large"); } } MathType new_range = ((TupleType)rangeType).getComponent( component ); FunctionType new_type = new FunctionType( domainType, new_range); int cnt = 0; int t_cnt = 0; for ( ii = 0; ii < component; ii++ ) { m_type = ((TupleType)rangeType).getComponent(ii); if ( m_type instanceof RealType ) { cnt++; } else { cnt += ((RealTupleType)m_type).getDimension(); t_cnt++; } } if ( new_range instanceof RealType ) { compSize = 1; flat_indices = new int[compSize]; flat_indices[0] = cnt; coord_sys = null; } else { compSize = ((RealTupleType)new_range).getDimension(); flat_indices = new int[ compSize ]; for ( jj = 0; jj < compSize; jj++ ) { flat_indices[jj] = cnt++; } coord_sys = RangeCoordinateSystems[t_cnt]; } ErrorEstimate[] errors_out = new ErrorEstimate[ compSize ]; Unit[] units_out = new Unit[ compSize ]; Set[] rangeSet_out = new Set[ compSize ]; boolean needDoubles = false; for ( ii = 0; ii < compSize; ii++ ) { units_out[ii] = RangeUnits[ flat_indices[ii] ]; errors_out[ii] = RangeErrors[ flat_indices[ii] ]; rangeSet_out[ii] = RangeSet[ flat_indices[ii] ]; if (rangeSet_out[ii] instanceof DoubleSet) needDoubles = true; } FlatField new_field = new FlatField( new_type, domainSet, coord_sys, null, rangeSet_out, units_out ); new_field.setRangeErrors( errors_out ); if (needDoubles) { double[][] new_values = new double[ compSize ][ n_samples ]; for ( ii = 0; ii < compSize; ii++ ) { new_values[ii] = unpackOneRangeComp( flat_indices[ii], copy); } new_field.setSamples( new_values, false ); } else { float[][] new_floats = new float[ compSize ][ n_samples ]; for ( ii = 0; ii < compSize; ii++ ) { new_floats[ii] = unpackOneFloatRangeComp( flat_indices[ii], copy); } new_field.setSamples( new_floats, false ); } return new_field; } public Data derivative( RealTuple location, RealType[] d_partial_s, MathType[] derivType_s, int error_mode ) throws VisADException, RemoteException { int ii, jj, kk, dd, rr, tt, pp, ss; Set domainSet = getDomainSet(); int domainDim = domainSet.getDimension(); int manifoldDimension = domainSet.getManifoldDimension(); int n_samples = domainSet.getLength(); CoordinateSystem d_coordsys = getDomainCoordinateSystem(); RealTupleType d_reference = (d_coordsys == null) ? null : d_coordsys.getReference(); MathType m_type = null; MathType[] m_types = null; RealType r_type = null; RealType[] r_types = null; TupleType t_type = null; boolean thisDomainFlag = true; if ( manifoldDimension != domainDim ) { throw new SetException("derivative: manifoldDimension must equal "+ "domain dimension" ); } error_mode = Data.NO_ERRORS; if ( location != null ) { thisDomainFlag = false; } RealTupleType domainType = ((FunctionType)Type).getDomain(); RealType[] r_comps = domainType.getRealComponents(); RealType[] r_compsRange = (((FunctionType)Type).getFlatRange()).getRealComponents(); RealType[] r_compsRef = (d_reference == null) ? null : d_reference.getRealComponents(); MathType RangeType = ((FunctionType)Type).getRange(); int n_partials; // number of partial derivatives to compute -* //- get all components for this function's domain -* if ( d_partial_s == null ) { n_partials = domainDim; d_partial_s = r_comps; } else { n_partials = d_partial_s.length; if ( n_partials > domainDim ) { throw new VisADException("derivative: too many d_partial components"); } } int[] u_index = new int[n_partials]; double[][] u_vectors = new double[n_partials][domainDim]; //- verify that input RealType-s match the Function's domain -* //- create unit vectors for the d_partial RealTypes -* int found = 0; int foundRef = 0; for ( ii = 0; ii < n_partials; ii++ ) { for ( jj = 0; jj < domainDim; jj++ ) { u_vectors[ii][jj] = 0d; if ( r_comps[jj].equals(d_partial_s[ii]) ) { u_index[ii] = jj; u_vectors[ii][jj] = 1d; found++; } else if ( d_reference != null ) { if ( r_compsRef[jj].equals(d_partial_s[ii]) ) { u_index[ii] = jj; u_vectors[ii][jj] = 1d; foundRef++; } } } } boolean transform; //- flag indicating coordinate transform is required --* if ( found == 0 ) { if ( foundRef == 0 ) { throw new VisADException("derivative: d_partial_s not in domain or reference"); } else if ( 0 < foundRef && foundRef < n_partials ) { throw new VisADException("derivative: d_partial_s must ALL be in function's "+ "domain or ALL in domain's reference"); } else { transform = true; } } else if ( 0 < found && found < n_partials ) { throw new VisADException("derivative: d_partial_s must ALL be in function's "+ "domain or ALL in domain's reference"); } else { transform = false; } Unit[] D_units = null; Unit[][] R_units = getRangeUnits(); String[][] derivNames = null; Unit[][] derivUnits = new Unit[ n_partials ][ TupleDimension ]; MathType[] new_range = new MathType[ n_partials ]; MathType[] new_types = new MathType[ n_partials ]; if ( !transform ) { D_units = domainSet.getSetUnits(); } else { D_units = d_reference.getDefaultUnits(); } //- Create derivative Units array -* for ( ii = 0; ii < n_partials; ii++ ) { for ( jj = 0; jj < TupleDimension; jj++ ) { if (( R_units == null)||( D_units == null )) { derivUnits[ii][jj] = null; } else if (( R_units[jj][0] == null )||( D_units[u_index[ii]] == null )) { derivUnits[ii][jj] = null; } else { derivUnits[ii][jj] = R_units[jj][0].divide( D_units[ u_index[ii] ] ); } } } if ( derivType_s == null ) { for ( ii = 0; ii < n_partials; ii++ ) { MathType M_type = Type.cloneDerivative( d_partial_s[ii] ); if ( thisDomainFlag ) { new_types[ii] = M_type; } else { new_types[ii] = ((FunctionType)M_type).getRange(); } } derivType_s = new_types; } else //- check supplied derivType-s for compatibility -* { if ( derivType_s.length != n_partials ) { throw new VisADException("derivative: must be a single MathType "+ "for each domain RealType"); } for ( ii = 0; ii < n_partials; ii++ ) { if ( thisDomainFlag ) { if ( !Type.equalsExceptName(derivType_s[ii]) ) { throw new TypeException("derivative: incompatible with function range"); } } else { if ( !((((FunctionType)Type).getRange()).equalsExceptName(derivType_s[ii])) ) { throw new TypeException("derivative: incompatible with function range"); } } } } double[][][] p_derivatives = null; DataImpl[] datums = new DataImpl[ n_partials ]; ErrorEstimate[][] rangeErrors_out = new ErrorEstimate[ n_partials ][ TupleDimension ]; if ( thisDomainFlag ) { p_derivatives = new double[ n_partials ][ TupleDimension ][ n_samples ]; for ( ii = 0; ii < n_partials; ii++ ) { datums[ii] = (DataImpl) cloneDouble( derivType_s[ii], derivUnits[ii], null ); } if ( isMissing() ) { if ( n_partials == 1 ) { return datums[0]; } else { return new Tuple( datums ); } } } else { p_derivatives = new double[ n_partials ][ TupleDimension ][ 1 ]; if ( isMissing() ) { for ( ii = 0; ii < n_partials; ii++ ) { for ( jj = 0; jj < TupleDimension; jj++ ) { p_derivatives[ii][jj][0] = Double.NaN; rangeErrors_out[ii][jj] = null; } } } } if ( !isMissing() ) //- skip computations. if thisDomainFlag is also set, then //- method would have already returned. { //- compute derivative-s, return FlatField or Tuple of FlatFields, or Data --* double[][] rangeValues = null; int[][] neighbors = null; int n_points; int n_index; int m_index; int index; float distance; float step; float f_sum; double d_sum; ErrorEstimate[] domainErrors = domainSet.getSetErrors(); //- Handle LinearSet case separately for efficiency -* if(( domainSet instanceof LinearSet )&&( thisDomainFlag )) { rangeValues = getValues(); //- each partial derivative -* for ( kk = 0; kk < n_partials; kk++ ) { //- get manifoldDimension index for this real axis ( LinearSet only ) -* m_index = u_index[kk]; //- get neigbors and separation along this axis -* neighbors = domainSet.getNeighbors( m_index ); step = (float) (((LinearSet)domainSet).getLinear1DComponent(m_index)).getStep(); //- compute derivative for each sample and each range component -* for ( ii = 0; ii < n_samples; ii++ ) { if ( neighbors[ii][0] == -1) { distance = step; n_index = neighbors[ii][1]; index = ii; } else if ( neighbors[ii][1] == -1 ) { distance = step; n_index = ii; index = neighbors[ii][0]; } else { distance = 2.f*step; n_index = neighbors[ii][1]; index = neighbors[ii][0]; } for ( rr = 0; rr < TupleDimension; rr++ ) { p_derivatives[kk][rr][ii] = ( rangeValues[rr][ n_index ] - rangeValues[rr][ index ] )/distance; } } if ( error_mode != Data.NO_ERRORS ) { for ( rr = 0; rr < TupleDimension; rr++ ) { double[] d_values = p_derivatives[kk][rr]; rangeErrors_out[kk][rr] = new ErrorEstimate( d_values, derivUnits[kk][rr], Data.DIVIDE, RangeErrors[rr], domainErrors[m_index], error_mode ); } } } neighbors = null; rangeValues = null; } else //- GriddedSet, IrregularSet --* { float dotproduct; float inv_dotproduct; float[][] weights = null; float sum_weights; float[][] Samples; //- compute derivative at this Set's sample locations --* if ( thisDomainFlag ) { rangeValues = getValues(); neighbors = new int[n_samples][]; weights = new float[n_samples][]; domainSet.getNeighbors( neighbors, weights ); if ( transform ) { Samples = domainSet.getSamples(true); Samples = CoordinateSystem.transformCoordinates( d_reference, null, null, null, domainType, d_coordsys, null, null, Samples ); } else { Samples = domainSet.getSamples(false); } } //- compute derivative at selected ( probably interpolated locations ) --* else { double[][] new_rangeValues; double[] d_array; int[][] new_neighbors; n_samples = 1; FlatField new_field; float[][] new_Samples; float[][] evalSamples; float[][] org_Samples = domainSet.getSamples(false); new_field = (FlatField) resample( new SingletonSet(location, null, null, null ), Data.WEIGHTED_AVERAGE, error_mode ); evalSamples = (new_field.getDomainSet()).getSamples(false); neighbors = new int[n_samples][]; weights = new float[n_samples][]; ((SimpleSet)domainSet).valueToInterp( evalSamples, neighbors, weights ); n_points = neighbors[0].length; new_neighbors = new int[n_samples][ n_points ]; new_rangeValues = new double[ TupleDimension ][ n_points + 1 ]; new_Samples = new float[ domainDim ][ n_points + 1 ]; for ( ii = 0; ii < domainDim; ii++ ) { new_Samples[ii][0] = evalSamples[ii][0]; } d_array = new_field.unpackValues(0); for ( ii = 0; ii < TupleDimension; ii++ ) { new_rangeValues[ii][0] = d_array[ii]; } for ( kk = 0; kk < n_points; kk++ ) { d_array = unpackValues( neighbors[0][kk] ); new_neighbors[0][kk] = kk + 1; for ( ii = 0; ii < TupleDimension; ii++ ) { new_rangeValues[ii][kk+1] = d_array[ii]; } for ( ii = 0; ii < domainDim; ii++ ) { new_Samples[ii][kk+1] = org_Samples[ii][ neighbors[0][kk] ]; } } neighbors = new_neighbors; rangeValues = new_rangeValues; Samples = new_Samples; if ( transform ) { Samples = CoordinateSystem.transformCoordinates( d_reference, null, null, null, domainType, d_coordsys, null, null, Samples ); } } //- compute derivatives for each sample --* for ( ii = 0; ii < n_samples; ii++ ) { n_points = neighbors[ii].length; double[] distances = new double[ n_points ]; double[][] derivatives = new double[ n_points ][ TupleDimension ]; double[][] uvecPoint = new double[ n_points ][ domainDim ]; //- neighbors loop -* for ( kk = 0; kk < n_points; kk++ ) { f_sum = 0; for ( dd = 0; dd < domainDim; dd++ ) { f_sum += (Samples[dd][ neighbors[ii][kk] ] - Samples[dd][ii])* (Samples[dd][ neighbors[ii][kk] ] - Samples[dd][ii]); uvecPoint[kk][dd] = Samples[dd][ neighbors[ii][kk] ] - Samples[dd][ii]; } distances[kk] = Math.sqrt((double)f_sum); for ( rr = 0; rr < TupleDimension; rr++ ) { derivatives[kk][rr] = rangeValues[rr][ neighbors[ii][kk] ] - rangeValues[rr][ii]; } } //- Interpolate for each partial derivative -* for ( pp = 0; pp < n_partials; pp++ ) { sum_weights = 0f; for ( kk = 0; kk < n_points; kk++ ) { dotproduct = 0; for ( dd = 0; dd < domainDim; dd++ ) { dotproduct += uvecPoint[kk][dd]*u_vectors[pp][dd]; } inv_dotproduct = 1f/dotproduct; if ( ! Float.isInfinite(inv_dotproduct) ) { sum_weights += weights[ii][kk]; for ( rr = 0; rr < TupleDimension; rr++ ) { p_derivatives[pp][rr][ii] += derivatives[kk][rr]*inv_dotproduct*weights[ii][kk]; } } } for ( rr = 0; rr < TupleDimension; rr++ ) { p_derivatives[pp][rr][ii] /= sum_weights; } } } } } //-not missing branch -* double[][] s_samples = null; for ( pp = 0; pp < n_partials; pp++ ) { s_samples = p_derivatives[pp]; if ( thisDomainFlag ) { ((FlatField)datums[pp]).setSamples( s_samples ); ((FlatField)datums[pp]).setRangeErrors( rangeErrors_out[pp] ); } else { MathType M_type = derivType_s[pp]; if ( M_type instanceof RealType ) { datums[pp] = new Real( (RealType) M_type, s_samples[0][0], derivUnits[pp][0], rangeErrors_out[pp][0]); } else if ( M_type instanceof RealTupleType ) { Real[] reals = new Real[ TupleDimension ]; for ( ii = 0; ii < TupleDimension; ii++ ) { reals[ii] = new Real( (RealType)((TupleType)M_type).getComponent(ii), s_samples[ii][0], derivUnits[pp][ii], rangeErrors_out[pp][0] ); } datums[pp] = new RealTuple( (RealTupleType) M_type, reals, RangeCoordinateSystem ); } else if ( M_type instanceof TupleType ) { ss = 0; tt = 0; int n_comps = ((TupleType)M_type).getDimension(); Data[] s_datums = new Data[ n_comps ]; for ( ii = 0; ii < n_comps; ii++ ) { m_type = ((TupleType)M_type).getComponent(ii); if ( m_type instanceof RealType ) { s_datums[ii] = new Real( (RealType)m_type, s_samples[tt][0], derivUnits[pp][ii], rangeErrors_out[pp][ii] ); tt++; } else if ( m_type instanceof RealTupleType ) { int n_compsI = ((TupleType)m_type).getDimension(); Real[] reals = new Real[ n_compsI ]; for ( jj = 0; jj < n_compsI; jj++ ) { reals[jj] = new Real( (RealType)((TupleType)m_type).getComponent(jj), s_samples[tt][0], derivUnits[pp][jj], rangeErrors_out[pp][jj] ); tt++; } s_datums[ii] = new RealTuple( (RealTupleType) m_type, reals, RangeCoordinateSystems[ss] ); ss++; } } datums[pp] = new Tuple( (TupleType)M_type, s_datums ); } } } if ( n_partials == 1 ) { return datums[0]; } else { return new Tuple( datums ); } } public Data derivative( int error_mode ) throws VisADException, RemoteException { MathType[] derivType_s = null; RealType[] d_partial_s = null; return derivative( null, d_partial_s, derivType_s, error_mode ); } public Data derivative( MathType[] derivType_s, int error_mode ) throws VisADException, RemoteException { return derivative( null, null, derivType_s, error_mode ); } public Function derivative( RealType d_partial, int error_mode ) throws VisADException, RemoteException { MathType[] derivType_s = null; RealType[] d_partial_s = new RealType[1]; d_partial_s[0] = d_partial; return (Function) derivative( null, d_partial_s, derivType_s, error_mode ); } public Function derivative( RealType d_partial, MathType derivType, int error_mode ) throws VisADException, RemoteException { MathType[] derivType_s = new MathType[1]; RealType[] d_partial_s = new RealType[1]; derivType_s[0] = derivType; d_partial_s[0] = d_partial; return (Function) derivative( null, d_partial_s, derivType_s, error_mode ); } /** * Resamples the range to domain samples of a given set. Resampling is either * by nearest neighbor or mulit-linear interpolation. NOTE: This code is very * similar to FieldImpl.resample(Set,int,int). * @param set The set of points at which to resample this * field. * @param sampling_mode Resampling mode: Data.NEAREST_NEIGHBOR or * Data.WEIGHTED_AVERAGE * @param error_mode Error estimation mode: Data.DEPENDENT, * Data.INDEPENDENT, or Data.NO_ERRORS. * @return Field of resampled data. RangeSet objects * in result are set to DoubleSet. NOTE: May * return this (i.e., not a copy). */ public Field resample (Set set, int sampling_mode, int error_mode) throws VisADException, RemoteException { /* NB: resampling is done in this method for a float domain. If * you make changes to this method, make the corresponding changes * in resampleDouble. */ Set domainSet = getDomainSet(); visad.util.Trace.call1("FlatField.resample"); if (domainSet.equals(set)) { // nothing to do visad.util.Trace.call2("FlatField.resample", "sampling set==domain set"); return this; } int dim = domainSet.getDimension(); if (dim != set.getDimension()) { throw new SetException("FlatField.resample: bad Set Dimension"); } if (domainSet instanceof GriddedDoubleSet) { return resampleDouble(set, sampling_mode, error_mode); } CoordinateSystem coord_sys = set.getCoordinateSystem(); Unit[] units = set.getSetUnits(); ErrorEstimate[] errors = (error_mode == NO_ERRORS) ? new ErrorEstimate[dim] : set.getSetErrors(); // create (initially missing) FlatField for return Set[] sets = new Set[TupleDimension]; for (int i=0; i getLength()) || !shouldBeDouble() || sampling_mode == WEIGHTED_AVERAGE) { values = unpackFloats(false); // values = Set.doubleToFloat(unpackValues()); } float[][] new_values = new float[TupleDimension][length]; float[] new_valuesJ; float[] valuesJ; if (sampling_mode == WEIGHTED_AVERAGE && domainSet instanceof SimpleSet) { // resample by interpolation int[][] indices = new int[length][]; float[][] coefs = new float[length][]; ((SimpleSet) domainSet).valueToInterp(vals, indices, coefs); /* DEBUG // System.out.println("DomainSet = " + domainSet); // System.out.println("set = " + set); // for (i=0; i 0) ii = 1; if (indices == null) ii = 0; for (i=0; i 0) { v = valuesJ[indices[i][0]] * coefs[i][0]; for (k=1; k 0) { float[][] xvals = new float[len][]; for (k = 0; k 0) { a = values[j][error_indices[2*i][0]] * error_coefs[2*i][0]; for (k=1; k 0) { b = values[j][error_indices[2*i+1][0]] * error_coefs[2*i+1][0]; for (k=1; k 0) ii = 1; if (indices == null) ii = 0; for (i=0; i = 0) { pr = true; if (i == 0) { System.out.println("DomainSet = " + domainSet); System.out.println("set = " + set); } System.out.println("NEAREST_NEIGHBOR indices[" + i + "] = " + indices[i]); } } */ // WLH 20 July 2000 if (values != null) { for (j=0; j = 0) ? valuesJ[indices[i]]: Float.NaN); } } } else { for (i=0; i = 0) { float[] xvals = unpackFloats(indices[i]); for (j=0; j = 0) ? values[j][error_indices[2*i]]: Double.NaN); float b = (float) ((error_indices[2*i+1] >= 0) ? values[j][error_indices[2*i+1]]: Double.NaN); sampling_partials[j][i] = Math.abs(b - a); } } } } if (sampling_errors) { for (j=0; j getLength()) || !shouldBeDouble() || sampling_mode == WEIGHTED_AVERAGE) { values = unpackFloats(false); // values = Set.doubleToFloat(unpackValues()); } double[][] new_values = new double[TupleDimension][length]; double[] new_valuesJ; float[] valuesJ; if (sampling_mode == WEIGHTED_AVERAGE) { // resample by interpolation int[][] indices = new int[length][]; double[][] coefs = new double[length][]; ((GriddedDoubleSet) domainSet).doubleToInterp(vals, indices, coefs); /* DEBUG // System.out.println("DomainSet = " + domainSet); // System.out.println("set = " + set); // for (i=0; i 0) ii = 1; if (indices == null) ii = 0; for (i=0; i 0) { v = valuesJ[indices[i][0]] * coefs[i][0]; for (k=1; k 0) { float[][] xvals = new float[len][]; for (k = 0; k 0) { a = values[j][error_indices[2*i][0]] * error_coefs[2*i][0]; for (k=1; k 0) { b = values[j][error_indices[2*i+1][0]] * error_coefs[2*i+1][0]; for (k=1; k 0) ii = 1; if (indices == null) ii = 0; for (i=0; i = 0) { pr = true; if (i == 0) { System.out.println("DomainSet = " + domainSet); System.out.println("set = " + set); } System.out.println("NEAREST_NEIGHBOR indices[" + i + "] = " + indices[i]); } } */ // WLH 20 July 2000 if (values != null) { for (j=0; j = 0) ? valuesJ[indices[i]]: Float.NaN); } } } else { for (i=0; i = 0) { float[] xvals = unpackFloats(indices[i]); for (j=0; j = 0) ? values[j][error_indices[2*i]]: Double.NaN); float b = (float) ((error_indices[2*i+1] >= 0) ? values[j][error_indices[2*i+1]]: Double.NaN); sampling_partials[j][i] = Math.abs(b - a); } } } } if (sampling_errors) { for (j=0; j = 0) any_mapped = true; } if (!any_mapped) return shadow; // check for any range coordinate systems boolean anyRangeRef = (RangeCoordinateSystem != null); if (RangeCoordinateSystems != null) { for (int i=0; i = 0 || anyRangeRef) { double min = Double.MAX_VALUE; double max = -Double.MAX_VALUE; if (valuesDI!=null) { for (int j=0; j = 0 && k < shadow.ranges[0].length) { shadow.ranges[0][k] = Math.min(shadow.ranges[0][k], min); shadow.ranges[1][k] = Math.max(shadow.ranges[1][k], max); } } } if (RangeCoordinateSystem != null) { // computeRanges for Reference (relative to range) RealTypes ShadowRealTupleType range_type = (ShadowRealTupleType) ((ShadowFunctionType) type).getRange(); shad_ref = range_type.getReference(); shadow = computeReferenceRanges(range_type, RangeCoordinateSystem, RangeUnits, shadow, shad_ref, ranges); } else if (RangeCoordinateSystems != null) { TupleType RangeType = (TupleType) ((FunctionType) Type).getRange(); int j = 0; for (int i=0; i Note that it is possible to simultaneously modify the domain-set of * both this instance and the clone by modifying the values in the array * returned by invoking getSamples(false)on the domain-set of * either this instance or the clone. Don't do this unless you enjoy * debugging. * * @return A clone of this instance. * @throws RuntimeException if a {@link VisADException} occurs. */ public Object clone() { FlatField clone; try { clone = (FlatField)super.clone(); } catch (CloneNotSupportedException ex) { throw new Error("Assertion failure"); // can't happen } synchronized(DoubleRange) { if (!MissingFlag) { try { // DRM 3 Oct 2006 - we need to recreate this otherwise // it's shared and it's mutable! clone.DoubleRange = new double[TupleDimension][]; double[][] values = unpackValues(true); clone.packValues(values, false); } catch (VisADException ex) { throw new RuntimeException(ex.toString()); } try { clone.setRangeErrors(RangeErrors); } catch (FieldException ex) { throw new Error("Assertion failure"); // can't happen } } } return clone; } String valuesString() throws VisADException { int rowlength; StringBuffer s = new StringBuffer(""); int ncolumns = 8 / TupleDimension; if (ncolumns < 1) ncolumns = 1; if (getDomainSet() instanceof GriddedSet) { // && ((GriddedSet) getDomainSet()).ManifoldDimension == 2) { rowlength = ((GriddedSet) getDomainSet()).getLength(0); } else { rowlength = getLength(); } RealTupleType range = ((FunctionType) Type).getFlatRange(); RealType[] types = range.getRealComponents(); double[][] values = unpackValues(); int rl = rowlength; int i = 0; while (ii) s.append(", "); if (TupleDimension == 1) { s.append(new Real(types[0], values[0][k], RangeUnits[0]).toString()); } else if (((FunctionType) Type).getReal()) { String t = "(" + new Real(types[0], values[0][k], RangeUnits[0]); for (int j=1; j 0) t = t + ", "; MathType type = RangeType.getComponent(l); if (type instanceof RealType) { t = t + new Real(types[j], values[j][k], RangeUnits[j]); j++; } else { int mm = ((TupleType) type).getDimension(); t = t + "(" + new Real(types[j], values[j][k], RangeUnits[j]); j++; for (int kk=1; kk Ads help maintain this website.