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.dht; 036 037import edu.emory.mathcs.utils.IOUtils; 038 039/** 040 * Accuracy check of double precision DHT's 041 * 042 * @author Piotr Wendykier (piotr.wendykier@gmail.com) 043 * 044 */ 045@SuppressWarnings("javadoc") 046public class AccuracyCheckDoubleDHT { 047 048 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, 049 1024, 1056, 2048, 8192, 10158, 16384, 32768, 65536, 131072 }; 050 051 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, 052 512, 1024 }; 053 054 private static int[] sizes3D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 128 }; 055 056 private static double eps = Math.pow(2, -52); 057 058 private AccuracyCheckDoubleDHT() { 059 060 } 061 062 public static void checkAccuracyDHT_1D() { 063 System.out.println("Checking accuracy of 1D DHT..."); 064 for (int i = 0; i < sizes1D.length; i++) { 065 DoubleDHT_1D fht = new DoubleDHT_1D(sizes1D[i]); 066 double err = 0.0; 067 double[] a = new double[sizes1D[i]]; 068 IOUtils.fillMatrix_1D(sizes1D[i], a); 069 double[] b = new double[sizes1D[i]]; 070 IOUtils.fillMatrix_1D(sizes1D[i], b); 071 fht.forward(a); 072 fht.inverse(a, true); 073 err = computeRMSE(a, b); 074 if (err > eps) { 075 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 076 } else { 077 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 078 } 079 a = null; 080 b = null; 081 fht = null; 082 System.gc(); 083 } 084 } 085 086 public static void checkAccuracyDHT_2D() { 087 System.out.println("Checking accuracy of 2D DHT (double[] input)..."); 088 for (int i = 0; i < sizes2D.length; i++) { 089 DoubleDHT_2D fht2 = new DoubleDHT_2D(sizes2D[i], sizes2D[i]); 090 double err = 0.0; 091 double[] a = new double[sizes2D[i] * sizes2D[i]]; 092 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 093 double[] b = new double[sizes2D[i] * sizes2D[i]]; 094 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], b); 095 fht2.forward(a); 096 fht2.inverse(a, true); 097 err = computeRMSE(a, b); 098 if (err > eps) { 099 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 100 } else { 101 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 102 } 103 a = null; 104 b = null; 105 fht2 = null; 106 System.gc(); 107 } 108 System.out.println("Checking accuracy of 2D DHT (double[][] input)..."); 109 for (int i = 0; i < sizes2D.length; i++) { 110 DoubleDHT_2D fht2 = new DoubleDHT_2D(sizes2D[i], sizes2D[i]); 111 double err = 0.0; 112 double[][] a = new double[sizes2D[i]][sizes2D[i]]; 113 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 114 double[][] b = new double[sizes2D[i]][sizes2D[i]]; 115 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], b); 116 fht2.forward(a); 117 fht2.inverse(a, true); 118 err = computeRMSE(a, b); 119 if (err > eps) { 120 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 121 } else { 122 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 123 } 124 a = null; 125 b = null; 126 fht2 = null; 127 System.gc(); 128 } 129 130 } 131 132 public static void checkAccuracyDHT_3D() { 133 System.out.println("Checking accuracy of 3D DHT (double[] input)..."); 134 for (int i = 0; i < sizes3D.length; i++) { 135 DoubleDHT_3D dht3 = new DoubleDHT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 136 double err = 0.0; 137 double[] a = new double[sizes3D[i] * sizes3D[i] * sizes3D[i]]; 138 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 139 double[] b = new double[sizes3D[i] * sizes3D[i] * sizes3D[i]]; 140 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], b); 141 dht3.forward(a); 142 dht3.inverse(a, true); 143 err = computeRMSE(a, b); 144 if (err > eps) { 145 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 146 + err); 147 } else { 148 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 149 + err); 150 } 151 a = null; 152 b = null; 153 dht3 = null; 154 System.gc(); 155 } 156 157 System.out.println("Checking accuracy of 3D DHT (double[][][] input)..."); 158 for (int i = 0; i < sizes3D.length; i++) { 159 DoubleDHT_3D dht3 = new DoubleDHT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 160 double err = 0.0; 161 double[][][] a = new double[sizes3D[i]][sizes3D[i]][sizes3D[i]]; 162 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 163 double[][][] b = new double[sizes3D[i]][sizes3D[i]][sizes3D[i]]; 164 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], b); 165 dht3.forward(a); 166 dht3.inverse(a, true); 167 err = computeRMSE(a, b); 168 if (err > eps) { 169 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 170 + err); 171 } else { 172 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 173 + err); 174 } 175 a = null; 176 b = null; 177 dht3 = null; 178 System.gc(); 179 } 180 } 181 182 private static double computeRMSE(double[] a, double[] b) { 183 if (a.length != b.length) { 184 throw new IllegalArgumentException("Arrays are not the same size"); 185 } 186 double rms = 0; 187 double tmp; 188 for (int i = 0; i < a.length; i++) { 189 tmp = (a[i] - b[i]); 190 rms += tmp * tmp; 191 } 192 return Math.sqrt(rms / a.length); 193 } 194 195 private static double computeRMSE(double[][] a, double[][] b) { 196 if (a.length != b.length || a[0].length != b[0].length) { 197 throw new IllegalArgumentException("Arrays are not the same size"); 198 } 199 double rms = 0; 200 double tmp; 201 for (int r = 0; r < a.length; r++) { 202 for (int c = 0; c < a[0].length; c++) { 203 tmp = (a[r][c] - b[r][c]); 204 rms += tmp * tmp; 205 } 206 } 207 return Math.sqrt(rms / (a.length * a[0].length)); 208 } 209 210 private static double computeRMSE(double[][][] a, double[][][] b) { 211 if (a.length != b.length || a[0].length != b[0].length || a[0][0].length != b[0][0].length) { 212 throw new IllegalArgumentException("Arrays are not the same size"); 213 } 214 double rms = 0; 215 double tmp; 216 for (int s = 0; s < a.length; s++) { 217 for (int r = 0; r < a[0].length; r++) { 218 for (int c = 0; c < a[0][0].length; c++) { 219 tmp = (a[s][r][c] - b[s][r][c]); 220 rms += tmp * tmp; 221 } 222 } 223 } 224 return Math.sqrt(rms / (a.length * a[0].length * a[0][0].length)); 225 } 226 227 public static void main(String[] args) { 228 checkAccuracyDHT_1D(); 229 checkAccuracyDHT_2D(); 230 checkAccuracyDHT_3D(); 231 System.exit(0); 232 } 233}