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 }