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