001    /*
002     * This file is part of McIDAS-V
003     *
004     * Copyright 2007-2013
005     * Space Science and Engineering Center (SSEC)
006     * University of Wisconsin - Madison
007     * 1225 W. Dayton Street, Madison, WI 53706, USA
008     * https://www.ssec.wisc.edu/mcidas
009     * 
010     * All Rights Reserved
011     * 
012     * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and
013     * some McIDAS-V source code is based on IDV and VisAD source code.  
014     * 
015     * McIDAS-V is free software; you can redistribute it and/or modify
016     * it under the terms of the GNU Lesser Public License as published by
017     * the Free Software Foundation; either version 3 of the License, or
018     * (at your option) any later version.
019     * 
020     * McIDAS-V is distributed in the hope that it will be useful,
021     * but WITHOUT ANY WARRANTY; without even the implied warranty of
022     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
023     * GNU Lesser Public License for more details.
024     * 
025     * You should have received a copy of the GNU Lesser Public License
026     * along with this program.  If not, see http://www.gnu.org/licenses.
027     */
028    package edu.wisc.ssec.mcidasv.data.hydra;
029    
030    import java.rmi.RemoteException;
031    
032    import org.apache.commons.math3.stat.correlation.PearsonsCorrelation;
033    import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
034    
035    import visad.Data;
036    import visad.FlatField;
037    import visad.FunctionType;
038    import visad.MathType;
039    import visad.Real;
040    import visad.RealTuple;
041    import visad.RealTupleType;
042    import visad.RealType;
043    import visad.TupleType;
044    import visad.VisADException;
045    
046    public class Statistics {
047    
048       DescriptiveStatistics[] descriptiveStats = null;
049       double[][] values_x;
050       double[][] rngVals;
051       int rngTupLen;
052       int numPoints;
053       int[] numGoodPoints;
054       MathType statType;
055    
056       PearsonsCorrelation pCorrelation = null;
057    
058    
059       public Statistics(FlatField fltFld) throws VisADException, RemoteException {
060          rngVals = fltFld.getValues(false);
061          rngTupLen = rngVals.length;
062          numPoints = fltFld.getDomainSet().getLength();
063          numGoodPoints = new int[rngTupLen];
064    
065          values_x = new double[rngTupLen][];
066    
067          for (int k=0; k<rngTupLen; k++) {
068            values_x[k] = removeMissing(rngVals[k]);
069            numGoodPoints[k] = values_x[k].length;
070          }
071    
072          descriptiveStats = new DescriptiveStatistics[rngTupLen];
073          for (int k=0; k<rngTupLen; k++) {
074            descriptiveStats[k] = new DescriptiveStatistics(values_x[k]);
075          }
076    
077          MathType rangeType = ((FunctionType)fltFld.getType()).getRange();
078    
079          if (rangeType instanceof RealTupleType) {
080            RealType[] rttypes = ((TupleType)rangeType).getRealComponents();
081            if (rngTupLen > 1) {
082              statType = new RealTupleType(rttypes);
083            }
084            else {
085              statType = (RealType) rttypes[0];
086            }
087          }
088          else if (rangeType instanceof RealType) {
089            statType = (RealType) rangeType;
090          }
091          else {
092             throw new VisADException("incoming type must be RealTupleType or RealType");
093          }
094    
095          pCorrelation = new PearsonsCorrelation();
096       }
097    
098       /** get the number of points in the domain of the FlatField
099       *
100       * @return number of points
101       */
102       public int numPoints() {
103         return numPoints;
104       }
105    
106       /** get the number of non-missing points in each range component
107       *
108       * @return number of non-missing points
109       */
110       public int[] getNumGoodPoints() {
111         return numGoodPoints;
112       }
113    
114       /* get the original range values
115       *
116       *@return the original range values
117       */
118       public double[][] getRngVals() {
119         return rngVals;
120       }
121    
122       /* get the range values actually used (missing removed)
123       *
124       *@return range values used
125       */
126    
127       public double[][] getValues() {
128         return values_x;
129       }
130    
131       private double[] removeMissing(double[] vals) {
132         int num = vals.length;
133         int cnt = 0;
134         int[] good = new int[num];
135         for (int k=0; k<num; k++) {
136            if ( !(Double.isNaN(vals[k])) ) {
137              good[cnt] = k;
138              cnt++;
139            }
140         }
141    
142         if (cnt == num) {
143            return vals;
144         }
145    
146         double[] newVals = new double[cnt];
147         for (int k=0; k<cnt; k++) {
148           newVals[k] = vals[good[k]];
149         }
150    
151         return newVals;
152       }
153    
154       private double[][] removeMissing(double[][] vals) {
155         int tupLen = vals.length;
156         double[][] newVals = new double[tupLen][];
157         for (int k=0; k < tupLen; k++) {
158            newVals[k] = removeMissing(vals[k]);
159         }
160         return newVals;
161       }
162    
163       public Data mean() throws VisADException, RemoteException {
164         double[] stats = new double[rngTupLen];
165         for (int k=0; k<rngTupLen; k++) {
166           stats[k] = descriptiveStats[k].getMean();
167         }
168         return makeStat(stats);
169       }
170    
171       public Data geometricMean() throws VisADException, RemoteException {
172         double[] stats = new double[rngTupLen];
173         for (int k=0; k<rngTupLen; k++) {
174           stats[k] = descriptiveStats[k].getGeometricMean();
175         }
176         return makeStat(stats);
177       }
178    
179    
180       public Data max() throws VisADException, RemoteException {
181         double[] stats = new double[rngTupLen];
182         for (int k=0; k<rngTupLen; k++) {
183           stats[k] = descriptiveStats[k].getMax();
184         }
185         return makeStat(stats);
186       }
187    
188       public Data min() throws VisADException, RemoteException {
189         double[] stats = new double[rngTupLen];
190         for (int k=0; k<rngTupLen; k++) {
191           stats[k] = descriptiveStats[k].getMin();
192         }
193         return makeStat(stats);
194       }
195    
196       public Data median() throws VisADException, RemoteException {
197         double[] stats = new double[rngTupLen];
198         for (int k=0; k<rngTupLen; k++) {
199           stats[k] = descriptiveStats[k].getPercentile(50.0);
200         }
201         return makeStat(stats);
202       }
203    
204       public Data variance() throws VisADException, RemoteException {
205         double[] stats = new double[rngTupLen];
206         for (int k=0; k<rngTupLen; k++) {
207           stats[k] = descriptiveStats[k].getVariance();
208         }
209         return makeStat(stats);
210       }
211    
212       public Data kurtosis() throws VisADException, RemoteException {
213         double[] stats = new double[rngTupLen];
214         for (int k=0; k<rngTupLen; k++) {
215           stats[k] = descriptiveStats[k].getKurtosis();
216         }
217         return makeStat(stats);
218       }
219    
220       public Data standardDeviation() throws VisADException, RemoteException {
221         double[] stats = new double[rngTupLen];
222         for (int k=0; k<rngTupLen; k++) {
223           stats[k] = descriptiveStats[k].getStandardDeviation();
224         }
225         return makeStat(stats);
226       }
227    
228       public Data skewness() throws VisADException, RemoteException {
229         double[] stats = new double[rngTupLen];
230         for (int k=0; k<rngTupLen; k++) {
231           stats[k] = descriptiveStats[k].getSkewness();
232         }
233         return makeStat(stats);
234       }
235    
236       public Data correlation(FlatField fltFld) throws VisADException, RemoteException {
237         double[][] values_x = this.rngVals;
238         double[][] values_y = fltFld.getValues(false);
239    
240         if (values_y.length != rngTupLen) {
241           throw new VisADException("both fields must have same range tuple length");
242         }
243    
244         double[] stats = new double[rngTupLen];
245         
246         for (int k=0; k<rngTupLen; k++) {
247           double[][] newVals = removeMissingAND(values_x[k], values_y[k]);
248           stats[k] = pCorrelation.correlation(newVals[0], newVals[1]);
249         }
250    
251         return makeStat(stats);
252       }
253    
254       private Data makeStat(double[] stats) throws VisADException, RemoteException {
255         if (statType instanceof RealType) {
256           return new Real((RealType)statType, stats[0]);
257         }
258         else if (statType instanceof RealTupleType) {
259           return new RealTuple((RealTupleType)statType, stats);
260         }
261         return null;
262       }
263    
264       private double[][] removeMissingAND(double[] vals_x, double[] vals_y) {
265         int cnt = 0;
266         int[] good = new int[vals_x.length];
267         for (int k=0; k<vals_x.length; k++) {
268           if ( !(Double.isNaN(vals_x[k])) && !(Double.isNaN(vals_y[k]))  ) {
269             good[cnt] = k;
270             cnt++;
271           }
272         }
273    
274         if (cnt == vals_x.length) {
275           return new double[][] {vals_x, vals_y};
276         }
277         else {
278           double[] newVals_x = new double[cnt];
279           double[] newVals_y = new double[cnt];
280           for (int k=0; k<cnt; k++) {
281             newVals_x[k] = vals_x[good[k]];
282             newVals_y[k] = vals_y[good[k]];
283           }
284           return new double[][] {newVals_x, newVals_y};
285         }
286       }
287    
288    }