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 }