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