001/* ***** BEGIN LICENSE BLOCK ***** 002 * Version: MPL 1.1/GPL 2.0/LGPL 2.1 003 * 004 * The contents of this file are subject to the Mozilla Public License Version 005 * 1.1 (the "License"); you may not use this file except in compliance with 006 * the License. You may obtain a copy of the License at 007 * http://www.mozilla.org/MPL/ 008 * 009 * Software distributed under the License is distributed on an "AS IS" basis, 010 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 011 * for the specific language governing rights and limitations under the 012 * License. 013 * 014 * The Original Code is JTransforms. 015 * 016 * The Initial Developer of the Original Code is 017 * Piotr Wendykier, Emory University. 018 * Portions created by the Initial Developer are Copyright (C) 2007-2009 019 * the Initial Developer. All Rights Reserved. 020 * 021 * Alternatively, the contents of this file may be used under the terms of 022 * either the GNU General Public License Version 2 or later (the "GPL"), or 023 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 024 * in which case the provisions of the GPL or the LGPL are applicable instead 025 * of those above. If you wish to allow use of your version of this file only 026 * under the terms of either the GPL or the LGPL, and not to allow others to 027 * use your version of this file under the terms of the MPL, indicate your 028 * decision by deleting the provisions above and replace them with the notice 029 * and other provisions required by the GPL or the LGPL. If you do not delete 030 * the provisions above, a recipient may use your version of this file under 031 * the terms of any one of the MPL, the GPL or the LGPL. 032 * 033 * ***** END LICENSE BLOCK ***** */ 034 035package edu.emory.mathcs.jtransforms.dst; 036 037import edu.emory.mathcs.utils.IOUtils; 038 039/** 040 * Accuracy check of double precision DST's 041 * 042 * @author Piotr Wendykier (piotr.wendykier@gmail.com) 043 * 044 */ 045@SuppressWarnings("javadoc") 046public class AccuracyCheckDoubleDST { 047 048 // private static int[] sizes1D = { 10158, 16384, 32768, 65536, 131072 }; 049 050 private static int[] sizes1D = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 120, 128, 256, 310, 512, 051 1024, 1056, 2048, 8192, 10158, 16384, 32768, 65536, 131072 }; 052 053 private static int[] sizes2D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 120, 128, 256, 310, 511, 054 512, 1024 }; 055 056 private static int[] sizes3D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 128 }; 057 058 private static double eps = Math.pow(2, -52); 059 060 private AccuracyCheckDoubleDST() { 061 062 } 063 064 public static void checkAccuracyDST_1D() { 065 System.out.println("Checking accuracy of 1D DST..."); 066 for (int i = 0; i < sizes1D.length; i++) { 067 DoubleDST_1D dst = new DoubleDST_1D(sizes1D[i]); 068 double err = 0; 069 double[] a = new double[sizes1D[i]]; 070 IOUtils.fillMatrix_1D(sizes1D[i], a); 071 double[] b = new double[sizes1D[i]]; 072 IOUtils.fillMatrix_1D(sizes1D[i], b); 073 dst.forward(a, true); 074 dst.inverse(a, true); 075 err = computeRMSE(a, b); 076 if (err > eps) { 077 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 078 } else { 079 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 080 } 081 a = null; 082 b = null; 083 dst = null; 084 System.gc(); 085 } 086 } 087 088 public static void checkAccuracyDST_2D() { 089 System.out.println("Checking accuracy of 2D DST (double[] input)..."); 090 for (int i = 0; i < sizes2D.length; i++) { 091 DoubleDST_2D dst2 = new DoubleDST_2D(sizes2D[i], sizes2D[i]); 092 double err = 0; 093 double[] a = new double[sizes2D[i] * sizes2D[i]]; 094 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 095 double[] b = new double[sizes2D[i] * sizes2D[i]]; 096 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], b); 097 dst2.forward(a, true); 098 dst2.inverse(a, true); 099 err = computeRMSE(a, b); 100 if (err > eps) { 101 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 102 } else { 103 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 104 } 105 a = null; 106 b = null; 107 dst2 = null; 108 System.gc(); 109 } 110 System.out.println("Checking accuracy of 2D DST (double[][] input)..."); 111 for (int i = 0; i < sizes2D.length; i++) { 112 DoubleDST_2D dst2 = new DoubleDST_2D(sizes2D[i], sizes2D[i]); 113 double err = 0; 114 double[][] a = new double[sizes2D[i]][sizes2D[i]]; 115 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 116 double[][] b = new double[sizes2D[i]][sizes2D[i]]; 117 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], b); 118 dst2.forward(a, true); 119 dst2.inverse(a, true); 120 err = computeRMSE(a, b); 121 if (err > eps) { 122 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 123 } else { 124 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 125 } 126 a = null; 127 b = null; 128 dst2 = null; 129 System.gc(); 130 } 131 132 } 133 134 public static void checkAccuracyDST_3D() { 135 System.out.println("Checking accuracy of 3D DST (double[] input)..."); 136 for (int i = 0; i < sizes3D.length; i++) { 137 DoubleDST_3D dst3 = new DoubleDST_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 138 double err = 0; 139 double[] a = new double[sizes3D[i] * sizes3D[i] * sizes3D[i]]; 140 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 141 double[] b = new double[sizes3D[i] * sizes3D[i] * sizes3D[i]]; 142 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], b); 143 dst3.forward(a, true); 144 dst3.inverse(a, true); 145 err = computeRMSE(a, b); 146 if (err > eps) { 147 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 148 + err); 149 } else { 150 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 151 + err); 152 } 153 a = null; 154 b = null; 155 dst3 = null; 156 System.gc(); 157 } 158 159 System.out.println("Checking accuracy of 3D DST (double[][][] input)..."); 160 for (int i = 0; i < sizes3D.length; i++) { 161 DoubleDST_3D dst3 = new DoubleDST_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 162 double err = 0; 163 double[][][] a = new double[sizes3D[i]][sizes3D[i]][sizes3D[i]]; 164 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 165 double[][][] b = new double[sizes3D[i]][sizes3D[i]][sizes3D[i]]; 166 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], b); 167 dst3.forward(a, true); 168 dst3.inverse(a, true); 169 err = computeRMSE(a, b); 170 if (err > eps) { 171 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 172 + err); 173 } else { 174 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 175 + err); 176 } 177 a = null; 178 b = null; 179 dst3 = null; 180 System.gc(); 181 } 182 } 183 184 private static double computeRMSE(double[] a, double[] b) { 185 if (a.length != b.length) { 186 throw new IllegalArgumentException("Arrays are not the same size."); 187 } 188 double rms = 0; 189 double tmp; 190 for (int i = 0; i < a.length; i++) { 191 tmp = (a[i] - b[i]); 192 rms += tmp * tmp; 193 } 194 return Math.sqrt(rms / a.length); 195 } 196 197 private static double computeRMSE(double[][] a, double[][] b) { 198 if (a.length != b.length || a[0].length != b[0].length) { 199 throw new IllegalArgumentException("Arrays are not the same size."); 200 } 201 double rms = 0; 202 double tmp; 203 for (int r = 0; r < a.length; r++) { 204 for (int c = 0; c < a[0].length; c++) { 205 tmp = (a[r][c] - b[r][c]); 206 rms += tmp * tmp; 207 } 208 } 209 return Math.sqrt(rms / (a.length * a[0].length)); 210 } 211 212 private static double computeRMSE(double[][][] a, double[][][] b) { 213 if (a.length != b.length || a[0].length != b[0].length || a[0][0].length != b[0][0].length) { 214 throw new IllegalArgumentException("Arrays are not the same size."); 215 } 216 double rms = 0; 217 double tmp; 218 for (int s = 0; s < a.length; s++) { 219 for (int r = 0; r < a[0].length; r++) { 220 for (int c = 0; c < a[0][0].length; c++) { 221 tmp = (a[s][r][c] - b[s][r][c]); 222 rms += tmp * tmp; 223 } 224 } 225 } 226 return Math.sqrt(rms / (a.length * a[0].length * a[0][0].length)); 227 } 228 229 public static void main(String[] args) { 230 checkAccuracyDST_1D(); 231 checkAccuracyDST_2D(); 232 checkAccuracyDST_3D(); 233 System.exit(0); 234 } 235}