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}