|
/Users/lyon/j4p/src/math/Mat1.java
|
1 package math;
2
3 import utils.Timer;
4
5 /**
6 * Mat1 is a class that uses 1-D arrays and contains static methods.
7 * @author Douglas Lyon
8 */
9 public abstract class Mat1 {
10
11 public static void main(final String[] args) {
12 //testAdd1();
13
14 //testAdd2();
15 for (int i=0; i < 10; i++)
16 testSum();
17 }
18
19
20 /**
21 * Add a vector times a scalar to a vector.
22 * The first vector contains the sum.
23 *
24 * @param a the first vector
25 * @param b the second vector
26 * @param c the scalar multiplier
27 */
28 public static void add(final double[] a,
29 final double[] b,
30 final double c) {
31 final int aLength = a.length;
32 if (aLength != b.length) {
33 System.out.println(
34 "ERROR: Vectors must be of equal length to add.");
35 return;
36 }
37 for (int i = 0; i < aLength; i++) {
38 a[i] += c * b[i];
39 }
40 }
41
42 /**
43 * Add two vectors. The first vector contains
44 * the sum.
45 * <b>This is not optimized for <i>altivec</i>
46 * </b>
47 *
48 * @param a the first vector <code>a = a + b; </code>
49 * @param b the second vector
50 */
51 public static void add(final double[] a,
52 final double[] b) {
53 final int aLength = a.length;
54 if (aLength != b.length) {
55 System.out.println(
56 "ERROR: Vectors must be of equal length to add.");
57 return;
58 }
59 for (int i = 0; i < aLength; i++) {
60 a[i] += b[i];
61 }
62 }
63
64 /**
65 * Add two vectors. The first vector contains
66 * the sum.
67 *
68 * @param a the first vector a = a + b;
69 * @param b the second vector
70 */
71 public static void add(final float[] a,
72 final float[] b) {
73 final int aLength = a.length;
74 if (aLength != b.length) {
75 System.out.println(
76 "ERROR: Vectors must be of equal length to add.");
77 return;
78 }
79 for (int i = 0; i < aLength; i++) {
80 a[i] += b[i];
81 }
82 }
83
84 public static void testSum() {
85 final int n = 1 << 19;
86 System.out.println("n=" + n);
87 final float[] a = ramp(0, 100, n);
88 final Timer t = new Timer();
89 double ans;
90 ans = sum(a);
91 ans = sum(a);
92 ans = sum(a);
93 t.start();
94 ans = sum(a);
95 t.stop();
96
97 t.print(1, "ans=" + ans);
98 ans = sum4(a);
99 ans = sum4(a);
100 ans = sum4(a);
101 t.start();
102 ans = sum4(a);
103 t.stop();
104 t.print(1, "ans4=" + ans);
105 ans = sum8(a);
106 ans = sum8(a);
107 ans = sum8(a);
108 t.start();
109 ans = sum8(a);
110 t.stop();
111 t.print(1, "ans8=" + ans);
112 ans = sumV(a);
113 ans = sumV(a);
114 ans = sumV(a);
115 t.start();
116 ans = sumV(a);
117 t.stop();
118 t.print(1, "sumV=" + ans);
119 }
120
121 public static void testAdd2() {
122 final int n = 1 << 12;
123 System.out.println("n=" + n);
124 final float[] a = ramp(0, 100, n);
125 final float[] b = ramp(0, 100, n);
126 final Timer t = new Timer();
127 double ans;
128 t.start();
129 add(a, b);
130 ans = add(a);
131 t.stop();
132
133 t.print(1, "ans=" + ans);
134
135 ans = sum8(a);
136 t.stop();
137 t.print(1, "ans=" + ans);
138 }
139
140 /**
141 * Calculate the dot product of two vectors.
142 *
143 * @param a the first vector
144 * @param b the second vector
145 * @return the dot product
146 */
147 static public double dot(final double[] a,
148 final double[] b) {
149 final int aLength = a.length;
150 if (aLength != b.length) {
151 System.out.println(
152 "ERROR: Vectors must be of equal length in dot product.");
153 return 0;
154 }
155 double sum = 0;
156 for (int i = 0; i < aLength; i++) {
157 sum += a[i] * b[i];
158 }
159 return sum;
160 }
161
162
163 /**
164 * Calculate the dot product of a vector with
165 * itself.
166 *
167 * @param a the vector
168 * @return the dot product
169 */
170 static public double dot(final double[] a) {
171 final int aLength = a.length;
172 double sum = 0;
173 for (int i = 0; i < aLength; i++) {
174 sum += a[i] * a[i];
175 }
176 return sum;
177 }
178
179 /**
180 * Calculate the magnitude a vector.
181 *
182 * @param a the vector
183 * @return the magnitude
184 */
185 static public double magnitude(
186 final double[] a) {
187 final int aLength = a.length;
188 double sum = 0;
189 for (int i = 0; i < aLength; i++) {
190 sum += a[i] * a[i];
191 }
192 return Math.sqrt(sum);
193 }
194
195 public static void print(final double[] a) {
196 for (int i = 0; i < a.length; i++) {
197
198 System.out.print(a[i] + " ");
199 System.out.println();
200 }
201 }
202
203 public static void printStats(
204 final String title,
205 final float[] a) {
206 System.out.println(title);
207 printStats(a);
208 }
209
210 public static void printStats(
211 final float[] a) {
212 float min = Float.MAX_VALUE;
213 float max = Float.MIN_VALUE;
214 float aBar = 0;
215
216
217 final double N = a.length;
218 for (int x = 0; x < a.length; x++) {
219
220 aBar += a[x];
221
222 min = Math.min(a[x], min);
223 max = Math.max(a[x], max);
224 }
225
226
227 aBar /= N;
228
229 System.out.println(" aBar=" + aBar +
230 " a min=" +
231 min +
232 " a max=" +
233 max +
234 " a.length=" +
235 a.length +
236 " a[0].length=" +
237 a.length);
238
239 }
240
241 public static void print(final float[] a) {
242 for (int i = 0; i < a.length; i++) {
243 System.out.print(a[i] + " ");
244 System.out.println();
245 }
246 }
247
248 public static final float sum(
249 final float[] a) {
250 float s = 0;
251 for (int i = 0; i < a.length; i++)
252 s += a[i];
253 return s;
254 }
255
256 public static final double sum(
257 final double[] a) {
258 double s = 0.0;
259 for (int i = 0; i < a.length; i++)
260 s += a[i];
261 return s;
262 }
263
264 public static float[] ramp(final float start,
265 final float end,
266 final int n) {
267 final float[] f = new float[n];
268 float t = 0f;
269 for (int i = 0; i < f.length; i++) {
270
271 f[i] = t * end + (1 - t) * start;
272 // t*p1 + (1-t)*p0
273 t = (1f * i) / n;
274 }
275 return f;
276
277 }
278
279 public static void testAdd1() {
280 final int n = 1 << 18;
281 System.out.println("n=" + n);
282 final float[] a = ramp(0, 100, n);
283 //final float[] b = ramp(0, 200, n);
284 //final float[] c = new float[a.length];
285 //final Timer t = new Timer();
286
287
288 testAdd1(a);
289
290
291 }
292
293 /**
294 * Show how loop unwinding can improve speed
295 * for a vector sum, but only to a point.
296 */
297 public static void testAdd1(final float a[]) {
298
299 int numberOfRuns = 100;
300 System.out.println(
301 "each run is executed " +
302 numberOfRuns +
303 " times, averages are reported in ms");
304 Runnable r = new Runnable() {
305 public void run() {
306 add(a);
307 }
308 };
309 System.out.println(" add: " +
310 Timer.benchMark(r,
311 numberOfRuns));
312
313 r = new Runnable() {
314 public void run() {
315 sumV(a);
316 }
317 };
318 System.out.println(" SumV: " +
319 Timer.benchMark(r,
320 numberOfRuns));
321
322
323 r = new Runnable() {
324 public void run() {
325 add16(a);
326 }
327 };
328 System.out.println(" add16:" +
329 Timer.benchMark(r,
330 numberOfRuns));
331
332 r = new Runnable() {
333 public void run() {
334 add8Depend(a);
335 }
336 };
337 System.out.println(" add8Depend:" +
338 Timer.benchMark(r,
339 numberOfRuns));
340
341 r = new Runnable() {
342 public void run() {
343 add8NoDepend(a);
344 }
345 };
346 System.out.println(" add8NoDepend:" +
347 Timer.benchMark(r,
348 numberOfRuns));
349
350
351 //add(a, b, c);
352 // ans = add8(c);
353 // t.print(1, "add c = a + b, ans=" + ans);
354
355
356 // add8(a, b, c);
357 // ans = add8(c);
358 // t.print(1,
359 // "add c = a + b, p; ans=" + ans);
360 }
361
362 public static void add(final float[] a,
363 final float[] b,
364 final float[] c) {
365 for (int i = 0; i < a.length; i++)
366 c[i] = a[i] + b[i];
367 }
368
369 public static void add8(final float[] a,
370 final float[] b,
371 final float[] c) {
372 int i1, i2, i3, i4, i5, i6, i7;
373 for (int i = 0; i < a.length; i += 8) {
374 i1 = i + 1;
375 i2 = i + 2;
376 i3 = i + 3;
377 i4 = i + 4;
378 i5 = i + 5;
379 i6 = i + 6;
380 i7 = i + 7;
381 c[i] = a[i] + b[i];
382 c[i1] = a[i1] + b[i1];
383 c[i2] = a[i2] + b[i2];
384 c[i3] = a[i3] + b[i3];
385 c[i4] = a[i4] + b[i4];
386 c[i5] = a[i5] + b[i5];
387 c[i6] = a[i6] + b[i6];
388 c[i7] = a[i7] + b[i7];
389 }
390 }
391
392 public static final float sum4(final float[] a) {
393 float s = 0;
394 for (int i = 0; i < a.length; i = i + 4) {
395 s = a[i] +
396 a[i + 1] +
397 a[i + 2] +
398 a[i + 3] +
399 s;
400 }
401 return s;
402 }
403
404 public static final float sum8(final float[] a) {
405 float s = 0;
406 for (int i = 0; i < a.length; i = i + 8) {
407 s = a[i] +
408 a[i + 1] +
409 a[i + 2] +
410 a[i + 3] +
411 a[i + 4] +
412 a[i + 5] +
413 a[i + 6] +
414 a[i + 7] +
415 s;
416 }
417 return s;
418 }
419
420 public static double add8Depend(
421 final float[] a) {
422 double s = 0;
423 for (int i = 0; i < a.length; i = i + 8) {
424 s = s + a[i];
425 s = s + a[i + 1];
426 s = s + a[i + 2];
427 s = s + a[i + 3];
428 s = s + a[i + 4];
429 s = s + a[i + 5];
430 s = s + a[i + 6];
431 s = s + a[i + 7];
432 }
433 return s;
434 }
435
436 public static double add8NoDepend(
437 final float[] a) {
438 final float[] s = new float[8];
439 for (int i = 0; i < a.length; i = i + 8) {
440 s[0] = s[0] + a[i];
441 s[1] = s[1] + a[i + 1];
442 s[2] = s[2] + a[i + 2];
443 s[3] = s[3] + a[i + 3];
444 s[4] = s[4] + a[i + 4];
445 s[5] = s[5] + a[i + 5];
446 s[6] = s[6] + a[i + 6];
447 s[7] = s[7] + a[i + 7];
448
449 }
450 // compacting phase.
451 return
452 s[0] +
453 s[1] +
454 s[2] +
455 s[3] +
456 s[4] +
457 s[5] +
458 s[6] +
459 s[7];
460 }
461
462 public static double add16(final float[] a) {
463 double s = 0;
464 for (int i = 0; i < a.length; i = i + 16) {
465 s = a[i] +
466 a[i + 1] +
467 a[i + 2] +
468 a[i + 3] +
469 a[i + 4] +
470 a[i + 5] +
471 a[i + 6] +
472 a[i + 7] +
473 a[i + 8] +
474 a[i + 9] +
475 a[i + 10] +
476 a[i + 11] +
477 a[i + 12] +
478 a[i + 13] +
479 a[i + 14] +
480 a[i + 15] +
481 s;
482 }
483 return s;
484 }
485
486
487 public static double add(final float[] a) {
488 double s = 0.0;
489 for (int i = 0; i < a.length; i++)
490 s = s + a[i];
491 return s;
492 }
493
494 public static double add(final short[] a) {
495 double s = 0.0;
496 if (a == null) return 0;
497 for (int i = 0; i < a.length; i++)
498 s += a[i];
499 return s;
500 }
501
502 public static short clip(final short i) {
503 if (i < 0)
504 return 0;
505
506 if (i > 255)
507 return 255;
508
509 return i;
510 }
511
512 public static final float sumV(
513 final float a[]) {
514 int i = 0;
515 int n = a.length;
516 float s[] = {0, 0, 0, 0};
517 for (i = 0; i < n; i = i + 4) {
518 s[0] = s[0] + a[i];
519 s[1] = s[1] + a[i + 1];
520 s[2] = s[2] + a[i + 2];
521 s[3] = s[3] + a[i + 3];
522 }// compaction phase...
523 return s[0] + s[1] + s[2] + s[3];
524 }
525
526 public static short[][] clip(
527 final short[][] i) {
528 final short[][] i2 = new short[i.length][i[0].length];
529 for (int x = 0; x < i2.length; x++)
530 for (int y = 0; y < i2[0].length; y++)
531 i2[x][y] = clip(i[x][y]);
532 return i2;
533 }
534
535 public static void normalize(
536 final double[] a) {
537 scale(a, 1.0 / sum(a));
538 }
539
540 public static float[] normalize(
541 final short[] a) {
542 return scale(a, (float) (1.0f / add(a)));
543 }
544
545 public static void normalize(final float[] a) {
546 scale(a, 1.0 / add(a));
547 }
548
549 public static double average(
550 final double[] a) {
551 final int n = a.length * a.length;
552 return sum(a) / n;
553 }
554
555 public static double average(final float[] a) {
556 final int n = a.length * a.length;
557 return add(a) / n;
558 }
559
560 public static double average(final short[] a) {
561 final int n = a.length * a.length;
562 return add(a) / n;
563 }
564
565 public static void threshold(final short[] a,
566 final short thresh) {
567 for (int i = 0; i < a.length; i++)
568 if (a[i] < thresh)
569 a[i] = 0;
570 else
571 a[i] = 255;
572 }
573
574 public static void threshold(final short[] a) {
575 threshold(a, (short) average(a));
576 }
577
578
579 public static float[] scale(final short[] a,
580 final float k) {
581 if (a == null) return null;
582 final float[] f = new float[a.length];
583 for (int i = 0; i < a.length; i++) {
584 f[i] = k * a[i];
585 }
586 return f;
587
588 }
589
590 public static void scale(final double[] a,
591 final double k) {
592 System.out.println(
593 "scale(double a[], double k)");
594 for (int i = 0; i < a.length; i++) {
595
596 a[i] *= k;
597 }
598
599 }
600
601 public static void scale(final float[] a,
602 final float k) {
603 if (a == null) return;
604 for (int i = 0; i < a.length; i++) {
605 a[i] = a[i] * k;
606 }
607 }
608
609 public static void scale(final float[] a,
610 final double k) {
611 //System.out.println("scale(float a[][], double k)");
612 // System.out.println("k="+k);
613 scale(a, (float) k);
614 }
615
616 public static float[] shortToFloat(
617 final short[] a) {
618 final int w = a.length;
619
620 final float[] c = new float[w];
621 for (int i = 0; i < w; i++)
622 c[i] = a[i];
623 return c;
624 }
625
626 public static short[] copyArray(
627 final short[] a) {
628 final int w = a.length;
629
630 final short[] c = new short[w];
631 for (int i = 0; i < w; i++)
632 c[i] = a[i];
633 return c;
634 }
635
636 public static double variance(final int[] a) {
637 final double xBar = mean(a);
638 double sum = 0;
639 double dx;
640 for (int i = 0; i < a.length; i++) {
641 dx = a[i] - xBar;
642 sum += dx * dx;
643 }
644 return sum / a.length;
645 }
646
647 public static double mean(final int[] a) {
648 double sum = 0;
649 for (int i = 0; i < a.length; i++)
650 sum += a[i];
651 return sum / a.length;
652 }
653
654 public static double coefficientOfVariation(
655 final int[] a) {
656 final double aBar = mean(a);
657 final double aBar2 = aBar * aBar;
658 return Math.sqrt(variance(a) / aBar2);
659 }
660
661 public static void testVariance() {
662 final int[] a = {1, 2, 3, 5, 4, 3, 2, 5,
663 6, 7};
664 System.out.println(
665 "The variance =" + variance(a));
666 }
667
668 public static void testCoefficientOfVariation() {
669 final int[] a = {0, 85, 87, 90, 100};
670 System.out.println("coefficientOfVariation({0,85,87,90,100}) ="
671 +
672 coefficientOfVariation(
673 a));
674 final int[] b = {95, 85, 87, 90, 100};
675 System.out.println("The coefficientOfVariation({95,85,87,90,100}) ="
676 +
677 coefficientOfVariation(
678 b));
679 }
680
681 public static boolean outlierHere(
682 final int[] a) {
683 return (coefficientOfVariation(a) > .1);
684 }
685
686 public static void testOutlier() {
687 final int[] a = {0, 85, 87, 90, 100};
688 final int[] b = {95, 85, 87, 90, 100};
689 System.out.println("dog ate my homework ={0,85,87,90,100}"
690 + outlierHere(a));
691 System.out.println("dog ate my homework ={95,85,87,90,100}"
692 + outlierHere(b));
693 }
694
695 public static void quickSort(final int[] a,
696 final int lo0,
697 final int hi0) {
698 // Based on the QuickSort method by
699 // James Gosling from Sun's SortDemo applet
700
701 int lo = lo0;
702 int hi = hi0;
703 final int mid;
704 int t;
705
706 if (hi0 > lo0) {
707 mid = a[(lo0 + hi0) / 2];
708 while (lo <= hi) {
709 while ((lo < hi0) &&
710 (a[lo] < mid))
711 ++lo;
712 while ((hi > lo0) &&
713 (a[hi] > mid))
714 --hi;
715 if (lo <= hi) {
716 t = a[lo];
717 a[lo] = a[hi];
718 a[hi] = t;
719 ++lo;
720 --hi;
721 }
722 }
723 if (lo0 < hi)
724 quickSort(a, lo0, hi);
725 if (lo < hi0)
726 quickSort(a, lo, hi0);
727
728 }
729 }
730
731 public static void swap(final int[] x,
732 final int a,
733 final int b) {
734 final int t = x[a];
735 x[a] = x[b];
736 x[b] = t;
737 }
738
739 public static void quickSort(final int[] a) {
740 quickSort(a, 0, a.length - 1);
741 }
742
743 public static double mean(final short[][] a) {
744 double sum = 0;
745 for (int x = 0; x < a.length; x++)
746 for (int y = 0; y < a[0].length; y++) {
747 sum += a[x][y];
748 }
749 return sum / (a.length * a[0].length);
750 }
751
752 public static double variance(
753 final short[][] a) {
754 final double xBar = mean(a);
755 double sum = 0;
756 double dx;
757 for (int x = 0; x < a.length; x++)
758 for (int y = 0; y < a[0].length; y++) {
759 dx = a[x][y] - xBar;
760 sum += dx * dx;
761 }
762 return sum / (a.length * a[0].length);
763 }
764
765 public static double[] getAverage(
766 final double[] a,
767 final double[] b,
768 final double[] c) {
769 final double[] avg = new double[a.length];
770
771 for (int i = 0; i < a.length; i++) {
772 avg[i] = (a[i] + b[i] + c[i]) / 3.0;
773 }
774 return avg;
775 }
776
777 public static short[][] copy(
778 final short[][] r) {
779 final short[][] c = new short[r.length][r[0].length];
780 final int w = r.length;
781 final int h = r[0].length;
782 for (int x = 0; x < w; x++)
783 for (int y = 0; y < h; y++)
784 c[x][y] = r[x][y];
785 return c;
786 }
787
788 public static short getMin(final short[] a) {
789 short min = 255;
790 for (int i = 0; i < a.length; i++)
791 if (a[i] < min)
792 min = a[i];
793 return min;
794 }
795
796 public static short getMax(final short[] a) {
797 short max = -255;
798 for (int i = 0; i < a.length; i++)
799 if (a[i] > max)
800 max = a[i];
801 return max;
802 }
803
804
805 public static void testQuickSort() {
806 final int[] a = {1, 2, 3, 5, 4, 3, 2, 5,
807 6, 7};
808 quickSort(a);
809 for (int i = 0; i < a.length; i++)
810 System.out.println(a[i]);
811 }
812
813 public static int numberOfNonZeros(
814 final short[][] k) {
815 final int umax = k.length;
816 final int vmax = k[0].length;
817 int sum = 0;
818
819 for (int x = 0; x < umax; x++)
820 for (int y = 0; y < vmax; y++)
821 if (k[x][y] != 0) sum++;
822 return sum;
823 }
824
825 public static void printMedian(
826 final short[][] k,
827 final String name) {
828 //printMaple(k);
829 System.out.println("\npublic void " +
830 name +
831 "(){\n"
832 + "\tfloat k[][] = {");
833 final int w = k.length;
834 final int h = k[0].length;
835
836 for (int y = 0; y < h; y++) {
837 System.out.print("\t{");
838 for (int x = 0; x < w - 1; x++)
839 System.out.print(k[x][y] + ", ");
840 String s = k[w - 1][y] + "}";
841 if (y < h - 1)
842 s = s + ",";
843 else
844 s = s + "};";
845 System.out.println(s);
846 }
847
848 final String s = "\n\tmedian(k);\n}";
849 System.out.println(s);
850
851 }
852
853
854 }
855