diff --git a/pykat/testing/test.py b/pykat/testing/test.py
old mode 100644
new mode 100755
index 87c93c60972c4e902417701454d51a4a67a41572..ead6ca71f4a9ccdda89c337ac95ebd552ac015dd
--- a/pykat/testing/test.py
+++ b/pykat/testing/test.py
@@ -47,7 +47,7 @@ def run_kat_file(item):
#try:
start = time.time()
- out,err = utils.runcmd([FINESSE_EXE, "--noheader", kat], cwd=SUITE_PATH)
+ out,err = utils.runcmd(["nice", "--18", FINESSE_EXE, "--noheader", kat], cwd=SUITE_PATH)
runtime = time.time()-start
OUT_FILE = os.path.join(SUITE_PATH,basename + ".out")
diff --git a/pykat/utilities/greedypoints/matched20.txt b/pykat/utilities/greedypoints/matched20.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a0c9fbcb1d14d6f3eb0ac44035c770c2f7fb50f8
--- /dev/null
+++ b/pykat/utilities/greedypoints/matched20.txt
@@ -0,0 +1,99 @@
+
+
+
+
+
+0.001000 0.000000 0 0
+0.001000 10000.000000 1 1
+0.258268 568.561873 1 11
+0.047776 5551.839465 1 15
+0.051117 501.672241 1 17
+0.001000 0.000000 0 20
+0.017706 334.448161 1 19
+0.014365 167.224080 1 19
+0.001000 0.000000 2 18
+0.014365 133.779264 5 19
+0.004341 0.000000 1 20
+0.014365 200.668896 3 20
+0.031070 4247.491639 5 17
+0.051117 3612.040134 7 17
+0.004341 0.000000 3 20
+0.011023 367.892977 7 20
+0.001000 0.000000 4 16
+0.001000 0.000000 4 20
+0.027729 2742.474916 7 16
+0.004341 0.000000 5 20
+0.001000 0.000000 1 19
+0.004341 0.000000 1 19
+0.024388 602.006689 11 17
+0.004341 0.000000 7 20
+0.004341 4481.605351 4 16
+0.014365 334.448161 7 15
+0.004341 0.000000 3 19
+0.001000 0.000000 6 18
+0.011023 535.117057 7 20
+0.004341 0.000000 9 20
+0.004341 0.000000 5 19
+0.017706 267.558528 15 19
+0.051117 301.003344 8 20
+0.001000 0.000000 9 18
+0.004341 0.000000 13 20
+0.004341 200.668896 11 19
+0.014365 568.561873 10 20
+0.001000 0.000000 7 17
+0.054458 2006.688963 18 19
+0.004341 0.000000 10 18
+0.001000 0.000000 10 20
+0.047776 3678.929766 19 20
+0.004341 0.000000 12 19
+0.007682 167.224080 9 20
+0.004341 0.000000 16 20
+0.011023 367.892977 13 20
+0.001000 0.000000 13 20
+0.017706 501.672241 12 20
+0.034411 1404.682274 19 20
+0.004341 0.000000 7 19
+0.004341 0.000000 18 20
+0.001000 0.000000 17 20
+0.007682 367.892977 18 20
+0.007682 66.889632 14 20
+0.007682 33.444816 12 18
+0.007682 100.334448 17 20
+0.001000 0.000000 20 20
+0.037753 903.010033 19 20
+0.004341 0.000000 20 20
+0.007682 66.889632 19 20
+0.014365 468.227425 20 20
+0.001000 0.000000 14 17
+0.011023 167.224080 20 20
+0.004341 66.889632 19 20
+0.024388 903.010033 20 20
+0.004341 33.444816 20 20
+0.021047 635.451505 19 20
+0.007682 200.668896 20 20
+0.044435 3712.374582 10 20
+0.004341 66.889632 20 20
+0.014365 66.889632 20 20
+0.004341 0.000000 19 20
+0.044435 802.675585 19 20
+0.001000 0.000000 19 20
+0.004341 33.444816 18 20
+0.014365 1036.789298 20 20
+0.007682 100.334448 20 20
+0.017706 702.341137 20 20
+0.001000 0.000000 17 18
+0.014365 167.224080 20 20
+0.007682 33.444816 20 20
+0.031070 1103.678930 20 20
+0.004341 0.000000 15 20
+0.007682 133.779264 20 20
+0.034411 0.000000 20 20
+0.004341 33.444816 17 19
+0.017706 434.782609 20 20
+0.001000 0.000000 11 18
+0.011023 234.113712 20 20
+0.007682 0.000000 20 20
+0.021047 869.565217 20 20
+0.064482 2240.802676 19 20
+0.004341 0.000000 11 19
+0.014365 367.892977 20 20
diff --git a/pykat/utilities/greedypoints/mismatched10.txt b/pykat/utilities/greedypoints/mismatched10.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c73211b349cccaaaf774dd0ae2f389736dd10a93
--- /dev/null
+++ b/pykat/utilities/greedypoints/mismatched10.txt
@@ -0,0 +1,237 @@
+
+
+
+
+
+0.001000 0.001000 -10000.000000 -10000.000000 0 0
+0.001000 0.911853 -10000.000000 -7058.823529 0 0
+0.001000 0.471118 -10000.000000 -2352.941176 0 0
+0.559265 0.471118 3529.411765 7058.823529 0 10
+0.001000 0.001000 10000.000000 -5882.352941 0 6
+1.000000 0.001000 -1176.470588 10000.000000 0 0
+0.001000 0.471118 10000.000000 -2352.941176 0 6
+0.001000 0.911853 -10000.000000 2352.941176 0 8
+0.001000 0.911853 -8823.529412 -1176.470588 0 2
+0.001000 0.001000 -7058.823529 -4117.647059 0 10
+0.001000 0.471118 10000.000000 -2352.941176 0 10
+0.001000 1.000000 -5294.117647 2352.941176 1 9
+0.001000 0.823706 -7647.058824 -1176.470588 0 4
+0.001000 0.500500 4705.882353 -2352.941176 0 4
+0.001000 0.001000 6470.588235 -4705.882353 8 10
+0.970618 0.001000 2352.941176 -9411.764706 8 10
+0.001000 0.764941 10000.000000 -1176.470588 0 6
+0.353588 0.382971 -8823.529412 -10000.000000 1 9
+0.001000 0.001000 7058.823529 -4705.882353 5 9
+0.001000 0.706176 10000.000000 -1176.470588 1 10
+0.001000 0.441735 -8823.529412 -2352.941176 6 10
+0.382971 0.382971 1176.470588 -4705.882353 1 9
+0.588647 0.441735 -7058.823529 1764.705882 1 9
+0.001000 0.001000 7058.823529 -4117.647059 2 10
+0.001000 1.000000 -10000.000000 -4705.882353 1 9
+0.001000 0.529882 10000.000000 -2352.941176 0 2
+0.001000 1.000000 -10000.000000 4705.882353 1 9
+0.001000 0.882471 10000.000000 -1176.470588 1 10
+0.001000 1.000000 -10000.000000 -4705.882353 1 2
+0.001000 1.000000 -10000.000000 4705.882353 1 2
+0.001000 0.941235 -10000.000000 -4705.882353 1 7
+0.001000 0.941235 10000.000000 4705.882353 1 7
+0.001000 0.823706 8235.294118 2352.941176 2 8
+0.001000 0.588647 10000.000000 -2352.941176 1 10
+0.001000 1.000000 7058.823529 2352.941176 1 2
+0.001000 0.001000 8235.294118 8823.529412 9 10
+0.001000 0.001000 7647.058824 8235.294118 8 9
+0.001000 0.001000 8235.294118 -2941.176471 3 10
+0.353588 0.382971 4705.882353 -1176.470588 3 9
+0.882471 0.941235 -3529.411765 -7058.823529 9 10
+1.000000 0.970618 -2352.941176 -4705.882353 2 3
+0.001000 0.823706 9411.764706 -1176.470588 5 10
+0.441735 0.471118 -2352.941176 -9411.764706 3 6
+0.823706 0.882471 -4705.882353 -4705.882353 8 9
+1.000000 0.970618 2352.941176 4705.882353 5 6
+1.000000 0.853088 -1176.470588 -4705.882353 3 4
+0.001000 1.000000 -7647.058824 2352.941176 10 10
+0.588647 0.706176 -9411.764706 2352.941176 8 9
+1.000000 0.823706 -2352.941176 -4705.882353 7 10
+0.001000 0.001000 5882.352941 -5294.117647 6 10
+0.001000 0.735559 8823.529412 -9411.764706 5 9
+0.500500 0.588647 -9411.764706 -2352.941176 5 8
+0.001000 0.911853 -4705.882353 4705.882353 1 5
+0.588647 0.706176 -9411.764706 1176.470588 8 9
+0.001000 1.000000 -3529.411765 -9411.764706 9 10
+0.001000 0.001000 7058.823529 -7647.058824 7 8
+0.970618 1.000000 -8823.529412 8235.294118 8 9
+0.001000 0.911853 -4705.882353 -4705.882353 1 5
+0.559265 0.676794 -9411.764706 -2352.941176 7 10
+0.001000 0.001000 9411.764706 10000.000000 9 10
+0.177294 0.147912 4117.647059 2941.176471 6 8
+1.000000 0.911853 -1176.470588 -4705.882353 8 10
+1.000000 1.000000 4705.882353 1176.470588 4 5
+0.676794 0.764941 9411.764706 4705.882353 3 9
+1.000000 0.764941 -9411.764706 1176.470588 1 9
+0.618029 0.647412 -9411.764706 -2352.941176 9 10
+0.001000 1.000000 -3529.411765 9411.764706 9 10
+0.823706 0.647412 9411.764706 -4705.882353 5 10
+0.001000 0.001000 8823.529412 -4117.647059 8 9
+0.177294 0.147912 -6470.588235 -8823.529412 6 8
+1.000000 1.000000 588.235294 6470.588235 8 9
+0.882471 0.882471 -3529.411765 -7647.058824 6 7
+1.000000 1.000000 1764.705882 -7647.058824 9 10
+0.706176 0.764941 7058.823529 -7647.058824 8 9
+0.970618 1.000000 -6470.588235 3529.411765 9 10
+0.265441 0.412353 -2352.941176 4705.882353 8 8
+0.823706 0.823706 -10000.000000 -8235.294118 9 10
+0.001000 0.970618 -7058.823529 2352.941176 7 9
+0.294824 0.294824 8823.529412 3529.411765 8 9
+0.001000 1.000000 -4117.647059 1176.470588 6 10
+0.823706 0.823706 8235.294118 10000.000000 9 10
+0.412353 0.412353 7647.058824 -8823.529412 9 10
+1.000000 1.000000 8823.529412 -9411.764706 9 10
+0.353588 0.382971 4705.882353 588.235294 3 9
+0.794324 0.794324 7647.058824 -8235.294118 9 10
+0.853088 0.647412 -9411.764706 1176.470588 10 10
+0.294824 0.294824 -8823.529412 -3529.411765 8 9
+0.001000 1.000000 3529.411765 -10000.000000 9 10
+0.676794 0.618029 4705.882353 9411.764706 9 10
+0.001000 1.000000 -4117.647059 -5294.117647 7 10
+0.500500 0.294824 -7058.823529 -1176.470588 0 10
+0.882471 0.882471 5294.117647 8823.529412 9 10
+0.001000 1.000000 8235.294118 7058.823529 2 10
+0.118529 0.030382 -9411.764706 7647.058824 3 6
+1.000000 0.970618 2352.941176 4705.882353 1 5
+1.000000 1.000000 3529.411765 -9411.764706 9 10
+1.000000 1.000000 3529.411765 7647.058824 9 10
+0.559265 0.764941 -10000.000000 8823.529412 10 10
+0.001000 1.000000 3529.411765 7647.058824 9 10
+0.001000 1.000000 -3529.411765 4705.882353 10 10
+0.882471 1.000000 -10000.000000 3529.411765 9 10
+0.970618 0.911853 -2352.941176 9411.764706 9 10
+0.618029 0.559265 4705.882353 -9411.764706 5 8
+0.001000 0.001000 5882.352941 1764.705882 3 9
+0.559265 0.588647 8235.294118 7647.058824 9 10
+0.001000 1.000000 4705.882353 -7058.823529 6 10
+0.764941 0.794324 6470.588235 5882.352941 9 10
+0.970618 1.000000 -3529.411765 4705.882353 9 10
+0.823706 0.823706 1764.705882 2941.176471 9 10
+0.001000 1.000000 4117.647059 -3529.411765 7 10
+0.294824 0.324206 -2352.941176 9411.764706 9 10
+1.000000 1.000000 6470.588235 7647.058824 9 10
+0.324206 0.265441 -4705.882353 -2352.941176 6 8
+1.000000 0.001000 -4705.882353 3529.411765 10 10
+0.001000 0.001000 -7058.823529 8823.529412 6 9
+1.000000 1.000000 5882.352941 9411.764706 9 10
+0.853088 0.647412 -4705.882353 9411.764706 10 10
+0.001000 1.000000 -4705.882353 -7647.058824 6 10
+0.001000 1.000000 -6470.588235 588.235294 4 10
+0.500500 0.412353 6470.588235 -7647.058824 6 10
+0.001000 1.000000 -3529.411765 -2352.941176 10 10
+0.001000 1.000000 -3529.411765 5294.117647 10 10
+1.000000 0.001000 -1764.705882 3529.411765 10 10
+0.529882 0.559265 -1176.470588 9411.764706 7 8
+0.823706 0.853088 -4117.647059 5882.352941 9 10
+0.001000 1.000000 -3529.411765 2352.941176 9 9
+0.001000 1.000000 5294.117647 8235.294118 5 10
+0.001000 1.000000 -3529.411765 -1176.470588 10 10
+0.588647 0.559265 -6470.588235 1764.705882 9 10
+0.001000 1.000000 4705.882353 -8823.529412 6 10
+0.001000 1.000000 -3529.411765 4117.647059 10 10
+0.941235 0.823706 -2352.941176 5882.352941 1 7
+0.177294 0.412353 -8823.529412 -8823.529412 8 9
+0.001000 1.000000 -5882.352941 3529.411765 5 10
+0.911853 0.794324 10000.000000 -8823.529412 8 10
+0.911853 0.970618 -10000.000000 -10000.000000 9 10
+0.706176 0.823706 -3529.411765 -6470.588235 6 9
+0.941235 0.970618 6470.588235 4117.647059 9 10
+0.001000 0.001000 5882.352941 7058.823529 10 10
+1.000000 0.001000 -4117.647059 3529.411765 10 10
+0.911853 0.706176 5294.117647 8823.529412 9 9
+0.001000 1.000000 6470.588235 8823.529412 4 10
+1.000000 0.001000 -5882.352941 -3529.411765 10 10
+0.001000 0.001000 -10000.000000 -10000.000000 10 10
+0.001000 1.000000 -5882.352941 6470.588235 4 10
+0.001000 1.000000 4117.647059 -5294.117647 9 10
+0.970618 0.764941 3529.411765 2941.176471 6 10
+1.000000 0.001000 -8235.294118 -3529.411765 10 10
+0.001000 1.000000 4117.647059 -7058.823529 4 10
+1.000000 1.000000 -8235.294118 9411.764706 10 10
+0.618029 0.706176 4117.647059 10000.000000 9 10
+0.001000 1.000000 -10000.000000 -588.235294 2 10
+0.001000 1.000000 6470.588235 -6470.588235 4 10
+0.001000 1.000000 -5294.117647 7058.823529 6 10
+1.000000 0.001000 1764.705882 3529.411765 10 10
+0.001000 1.000000 4117.647059 10000.000000 8 10
+0.001000 0.001000 4705.882353 5882.352941 10 10
+0.001000 1.000000 -5294.117647 2941.176471 6 10
+0.001000 1.000000 4117.647059 7647.058824 9 10
+0.001000 1.000000 -3529.411765 -2941.176471 10 10
+1.000000 1.000000 -5294.117647 5882.352941 10 10
+0.001000 1.000000 4117.647059 1176.470588 9 10
+0.911853 0.970618 -10000.000000 8823.529412 8 9
+0.001000 0.001000 -9411.764706 -8235.294118 10 10
+0.001000 1.000000 -4117.647059 8235.294118 10 10
+1.000000 1.000000 1176.470588 8823.529412 10 10
+0.001000 1.000000 -4117.647059 9411.764706 10 10
+0.001000 1.000000 -4117.647059 6470.588235 10 10
+0.001000 1.000000 5882.352941 5882.352941 5 10
+0.001000 1.000000 -4705.882353 -7647.058824 8 10
+0.001000 1.000000 4117.647059 8823.529412 10 10
+1.000000 0.001000 -9411.764706 4117.647059 10 10
+0.001000 1.000000 -6470.588235 -10000.000000 5 10
+1.000000 0.001000 -4705.882353 4117.647059 10 10
+1.000000 0.001000 3529.411765 4705.882353 10 10
+0.001000 0.001000 -8235.294118 -7058.823529 10 10
+0.001000 1.000000 4705.882353 8235.294118 10 10
+0.001000 1.000000 -4117.647059 -2352.941176 10 10
+0.001000 1.000000 -4705.882353 10000.000000 10 10
+0.001000 1.000000 4705.882353 588.235294 9 10
+0.001000 1.000000 -5882.352941 -1176.470588 6 10
+1.000000 0.001000 -6470.588235 -4705.882353 10 10
+1.000000 0.001000 2941.176471 4705.882353 10 10
+1.000000 0.001000 -7058.823529 -4705.882353 10 10
+0.001000 1.000000 -4705.882353 -3529.411765 10 10
+0.001000 1.000000 -4705.882353 7058.823529 10 10
+0.001000 1.000000 -4705.882353 -8823.529412 9 10
+1.000000 0.001000 2352.941176 4705.882353 10 10
+0.001000 1.000000 -4705.882353 4705.882353 10 10
+1.000000 0.001000 -8235.294118 -4705.882353 10 10
+0.001000 1.000000 -4705.882353 -2941.176471 10 10
+0.001000 0.001000 4117.647059 4705.882353 10 10
+0.001000 1.000000 -4705.882353 4117.647059 10 10
+0.001000 0.001000 8823.529412 9411.764706 10 10
+0.001000 1.000000 -4705.882353 -4117.647059 10 10
+1.000000 0.001000 -5294.117647 -4705.882353 10 10
+0.001000 1.000000 -4705.882353 -588.235294 10 10
+0.001000 1.000000 -5294.117647 8823.529412 9 10
+0.001000 1.000000 -4705.882353 5294.117647 10 10
+0.001000 0.001000 -7058.823529 -5294.117647 10 10
+0.001000 1.000000 4705.882353 9411.764706 10 10
+0.001000 1.000000 -4705.882353 -1764.705882 10 10
+1.000000 0.001000 -5882.352941 -4705.882353 10 10
+1.000000 0.001000 -8823.529412 -5294.117647 10 10
+0.001000 1.000000 5294.117647 6470.588235 9 10
+1.000000 0.001000 1176.470588 4705.882353 10 10
+1.000000 0.001000 -9411.764706 -5294.117647 10 10
+0.001000 1.000000 4705.882353 7647.058824 10 10
+0.001000 1.000000 -9411.764706 1764.705882 3 10
+0.001000 1.000000 -5294.117647 588.235294 10 10
+0.001000 0.001000 -10000.000000 -9411.764706 10 10
+0.001000 1.000000 5294.117647 10000.000000 10 10
+1.000000 0.001000 -10000.000000 -5294.117647 10 10
+0.001000 1.000000 -5294.117647 -588.235294 10 10
+0.001000 1.000000 -5882.352941 -1176.470588 9 10
+0.001000 1.000000 -5882.352941 -7058.823529 9 10
+0.001000 1.000000 5882.352941 -6470.588235 9 10
+0.001000 1.000000 5882.352941 4117.647059 9 10
+0.001000 0.001000 -8235.294118 -6470.588235 10 10
+0.001000 1.000000 -5882.352941 4705.882353 9 10
+0.001000 1.000000 5294.117647 5882.352941 10 10
+0.001000 1.000000 5882.352941 -5882.352941 9 10
+0.001000 1.000000 -5882.352941 -4705.882353 9 10
+1.000000 0.001000 -2352.941176 5294.117647 10 10
+0.001000 1.000000 -5882.352941 -2941.176471 9 10
+1.000000 0.001000 -7647.058824 -5294.117647 10 10
+0.001000 0.001000 4117.647059 5294.117647 10 10
+0.001000 1.000000 -5882.352941 -4117.647059 9 10
+0.001000 1.000000 -6470.588235 -5294.117647 8 10
+0.001000 0.001000 -8235.294118 7647.058824 10 10
+0.001000 1.000000 5882.352941 7058.823529 9 10
+0.001000 1.000000 -6470.588235 7647.058824 8 10
diff --git a/pykat/utilities/greedypoints/mismatched20.txt b/pykat/utilities/greedypoints/mismatched20.txt
new file mode 100644
index 0000000000000000000000000000000000000000..205e3761be1348d886a2297e4559ccfbfae3df52
--- /dev/null
+++ b/pykat/utilities/greedypoints/mismatched20.txt
@@ -0,0 +1,249 @@
+
+
+
+
+
+0.001000 0.001000 0.000000 0.000000 0 0
+0.001000 0.001000 4000.000000 4000.000000 1 1
+0.013000 0.001000 689.655172 0.000000 14 20
+0.001000 0.001000 275.862069 1103.448276 1 15
+0.001000 0.001000 3724.137931 137.931034 1 3
+0.006333 0.001000 3172.413793 551.724138 3 3
+0.001000 0.001000 137.931034 4000.000000 3 3
+0.003667 0.001000 3724.137931 137.931034 11 15
+0.001000 0.003667 965.517241 3172.413793 5 9
+0.001000 0.001000 137.931034 3310.344828 15 19
+0.010333 0.001000 413.793103 0.000000 12 12
+0.006333 0.001000 3172.413793 551.724138 3 15
+0.001000 0.001000 3724.137931 137.931034 1 7
+0.007667 0.001000 3448.275862 275.862069 1 1
+0.001000 0.010333 137.931034 3862.068966 1 7
+0.001000 0.003667 551.724138 3862.068966 1 13
+0.007667 0.001000 4000.000000 2068.965517 3 13
+0.006333 0.001000 3448.275862 275.862069 1 9
+0.001000 0.001000 275.862069 1103.448276 8 11
+0.001000 0.001000 275.862069 3862.068966 2 3
+0.001000 0.001000 689.655172 0.000000 0 8
+0.001000 0.001000 827.586207 275.862069 1 1
+0.001000 0.001000 2896.551724 965.517241 1 11
+0.001000 0.002333 1655.172414 2758.620690 1 17
+0.001000 0.001000 413.793103 137.931034 1 4
+0.001000 0.002333 413.793103 2482.758621 1 5
+0.002333 0.001000 3172.413793 413.793103 1 1
+0.001000 0.001000 689.655172 3172.413793 5 5
+0.001000 0.001000 689.655172 0.000000 0 6
+0.001000 0.001000 137.931034 413.793103 1 1
+0.003667 0.001000 3724.137931 137.931034 15 19
+0.009000 0.001000 137.931034 0.000000 20 20
+0.001000 0.001000 3586.206897 827.586207 1 5
+0.001000 0.001000 1241.379310 275.862069 1 1
+0.001000 0.001000 413.793103 3586.206897 1 1
+0.001000 0.006333 965.517241 3172.413793 3 18
+0.001000 0.001000 413.793103 137.931034 1 12
+0.001000 0.001000 275.862069 2482.758621 2 3
+0.001000 0.001000 827.586207 275.862069 1 10
+0.002333 0.002333 551.724138 3862.068966 4 19
+0.006333 0.001000 2482.758621 689.655172 1 1
+0.001000 0.001000 137.931034 413.793103 6 9
+0.001000 0.001000 413.793103 137.931034 1 10
+0.001000 0.001000 3448.275862 275.862069 1 13
+0.001000 0.001000 689.655172 2206.896552 1 7
+0.003667 0.001000 689.655172 137.931034 2 3
+0.001000 0.001000 2896.551724 1241.379310 1 15
+0.001000 0.001000 689.655172 0.000000 0 4
+0.001000 0.001000 2206.896552 413.793103 1 4
+0.002333 0.001000 689.655172 137.931034 2 7
+0.001000 0.003667 3172.413793 551.724138 3 8
+0.001000 0.003667 137.931034 3862.068966 19 19
+0.001000 0.001000 1241.379310 275.862069 1 10
+0.001000 0.011667 1931.034483 3448.275862 1 17
+0.001000 0.001000 413.793103 1655.172414 1 7
+0.006333 0.001000 3724.137931 137.931034 13 18
+0.001000 0.001000 275.862069 1793.103448 2 3
+0.001000 0.001000 3586.206897 689.655172 1 3
+0.001000 0.003667 137.931034 3862.068966 6 11
+0.001000 0.001000 1655.172414 413.793103 1 1
+0.010333 0.005000 1793.103448 0.000000 19 20
+0.001000 0.001000 1241.379310 2896.551724 3 9
+0.001000 0.003667 275.862069 1103.448276 3 9
+0.001000 0.001000 551.724138 2344.827586 1 1
+0.001000 0.001000 137.931034 689.655172 2 3
+0.001000 0.001000 413.793103 137.931034 1 20
+0.002333 0.001000 413.793103 137.931034 1 2
+0.001000 0.001000 1793.103448 275.862069 1 3
+0.001000 0.001000 3034.482759 1103.448276 0 1
+0.001000 0.001000 137.931034 413.793103 11 15
+0.001000 0.001000 275.862069 827.586207 2 3
+0.002333 0.010333 2068.965517 4000.000000 5 19
+0.002333 0.001000 827.586207 275.862069 1 6
+0.001000 0.001000 137.931034 689.655172 5 5
+0.001000 0.001000 689.655172 137.931034 1 10
+0.013000 0.001000 3448.275862 137.931034 16 17
+0.001000 0.002333 551.724138 1793.103448 1 3
+0.011667 0.006333 0.000000 0.000000 1 15
+0.001000 0.003667 137.931034 3586.206897 17 18
+0.001000 0.005000 3034.482759 1793.103448 1 20
+0.005000 0.001000 2482.758621 827.586207 1 1
+0.005000 0.001000 965.517241 3310.344828 14 15
+0.010333 0.005000 0.000000 0.000000 1 13
+0.005000 0.001000 827.586207 137.931034 7 9
+0.001000 0.003667 0.000000 0.000000 17 17
+0.001000 0.002333 137.931034 275.862069 2 3
+0.001000 0.006333 137.931034 827.586207 10 11
+0.001000 0.001000 413.793103 137.931034 1 7
+0.001000 0.001000 1931.034483 551.724138 1 1
+0.001000 0.006333 275.862069 1103.448276 1 9
+0.005000 0.001000 1241.379310 137.931034 3 8
+0.001000 0.001000 275.862069 137.931034 3 10
+0.001000 0.001000 1379.310345 413.793103 1 1
+0.001000 0.001000 137.931034 413.793103 13 13
+0.005000 0.001000 1241.379310 137.931034 3 11
+0.001000 0.006333 3448.275862 1931.034483 1 20
+0.001000 0.001000 965.517241 0.000000 0 2
+0.003667 0.001000 1241.379310 275.862069 1 6
+0.001000 0.003667 137.931034 3724.137931 8 9
+0.013000 0.003667 1241.379310 0.000000 15 17
+0.002333 0.010333 689.655172 3862.068966 3 3
+0.001000 0.001000 413.793103 1379.310345 1 1
+0.001000 0.005000 137.931034 551.724138 16 16
+0.001000 0.001000 413.793103 1931.034483 0 1
+0.003667 0.001000 551.724138 137.931034 8 9
+0.001000 0.001000 137.931034 413.793103 9 16
+0.013000 0.001000 3724.137931 137.931034 17 19
+0.001000 0.001000 2206.896552 413.793103 0 3
+0.005000 0.001000 137.931034 0.000000 19 19
+0.005000 0.001000 1241.379310 965.517241 1 10
+0.001000 0.007667 137.931034 965.517241 7 8
+0.001000 0.007667 827.586207 2068.965517 3 20
+0.001000 0.001000 137.931034 275.862069 8 9
+0.002333 0.001000 0.000000 0.000000 19 20
+0.001000 0.001000 275.862069 137.931034 3 7
+0.001000 0.003667 275.862069 1103.448276 8 9
+0.001000 0.010333 137.931034 1241.379310 10 12
+0.001000 0.001000 137.931034 689.655172 11 12
+0.001000 0.001000 275.862069 137.931034 1 2
+0.001000 0.005000 275.862069 1103.448276 3 14
+0.002333 0.001000 2206.896552 689.655172 8 9
+0.003667 0.001000 0.000000 0.000000 14 19
+0.005000 0.005000 0.000000 0.000000 8 18
+0.005000 0.001000 1793.103448 275.862069 2 5
+0.001000 0.003667 137.931034 3448.275862 15 19
+0.001000 0.001000 689.655172 137.931034 1 16
+0.001000 0.001000 137.931034 413.793103 18 19
+0.002333 0.001000 1241.379310 137.931034 9 13
+0.001000 0.001000 413.793103 3724.137931 3 3
+0.011667 0.002333 965.517241 0.000000 19 20
+0.001000 0.001000 965.517241 275.862069 1 2
+0.001000 0.001000 3724.137931 137.931034 1 20
+0.001000 0.002333 137.931034 275.862069 7 13
+0.001000 0.001000 137.931034 413.793103 8 8
+0.006333 0.005000 0.000000 0.000000 8 20
+0.001000 0.001000 965.517241 137.931034 1 4
+0.011667 0.001000 3586.206897 137.931034 7 17
+0.002333 0.001000 413.793103 137.931034 2 19
+0.003667 0.001000 827.586207 2068.965517 3 20
+0.001000 0.001000 275.862069 1793.103448 7 9
+0.001000 0.001000 137.931034 689.655172 8 8
+0.001000 0.001000 137.931034 1241.379310 4 5
+0.001000 0.002333 137.931034 275.862069 9 20
+0.002333 0.001000 0.000000 0.000000 13 19
+0.002333 0.001000 0.000000 0.000000 2 18
+0.011667 0.001000 3724.137931 137.931034 18 18
+0.005000 0.006333 0.000000 0.000000 18 19
+0.003667 0.001000 551.724138 137.931034 8 12
+0.001000 0.003667 137.931034 551.724138 10 19
+0.010333 0.001000 3448.275862 413.793103 2 13
+0.005000 0.005000 0.000000 0.000000 11 18
+0.001000 0.003667 137.931034 551.724138 9 14
+0.002333 0.001000 275.862069 137.931034 7 18
+0.001000 0.007667 137.931034 3586.206897 17 17
+0.003667 0.001000 413.793103 137.931034 16 18
+0.009000 0.005000 0.000000 0.000000 20 20
+0.001000 0.003667 137.931034 275.862069 5 15
+0.001000 0.002333 275.862069 827.586207 2 16
+0.007667 0.002333 4000.000000 137.931034 11 19
+0.001000 0.007667 137.931034 965.517241 13 19
+0.003667 0.007667 137.931034 3586.206897 19 19
+0.003667 0.001000 413.793103 137.931034 13 18
+0.002333 0.007667 0.000000 0.000000 19 20
+0.010333 0.007667 3586.206897 1793.103448 3 18
+0.001000 0.002333 137.931034 413.793103 13 20
+0.001000 0.003667 137.931034 413.793103 5 19
+0.001000 0.001000 275.862069 137.931034 9 16
+0.001000 0.011667 137.931034 0.000000 1 20
+0.002333 0.003667 0.000000 0.000000 14 18
+0.009000 0.003667 0.000000 137.931034 9 20
+0.001000 0.010333 137.931034 3724.137931 18 19
+0.002333 0.001000 0.000000 0.000000 13 18
+0.003667 0.002333 0.000000 0.000000 19 20
+0.001000 0.001000 137.931034 413.793103 14 16
+0.011667 0.003667 0.000000 137.931034 19 20
+0.010333 0.001000 3172.413793 551.724138 3 10
+0.002333 0.013000 137.931034 0.000000 16 20
+0.009000 0.001000 0.000000 137.931034 19 20
+0.001000 0.001000 137.931034 275.862069 13 16
+0.013000 0.002333 0.000000 137.931034 18 18
+0.003667 0.002333 0.000000 0.000000 18 19
+0.002333 0.013000 137.931034 0.000000 11 17
+0.003667 0.001000 551.724138 137.931034 5 8
+0.001000 0.001000 275.862069 137.931034 1 16
+0.001000 0.003667 137.931034 413.793103 10 19
+0.013000 0.002333 137.931034 137.931034 15 20
+0.001000 0.001000 0.000000 0.000000 17 20
+0.005000 0.006333 0.000000 0.000000 20 20
+0.002333 0.001000 275.862069 137.931034 15 18
+0.001000 0.002333 137.931034 275.862069 11 18
+0.001000 0.002333 137.931034 413.793103 3 18
+0.001000 0.013000 137.931034 4000.000000 19 20
+0.003667 0.005000 0.000000 0.000000 19 20
+0.013000 0.002333 0.000000 137.931034 20 20
+0.005000 0.006333 137.931034 0.000000 19 20
+0.005000 0.001000 413.793103 137.931034 18 18
+0.011667 0.002333 4000.000000 137.931034 19 19
+0.002333 0.002333 0.000000 0.000000 19 20
+0.003667 0.010333 137.931034 0.000000 20 20
+0.001000 0.002333 137.931034 275.862069 19 20
+0.007667 0.003667 0.000000 137.931034 19 20
+0.002333 0.007667 137.931034 4000.000000 20 20
+0.002333 0.002333 0.000000 0.000000 16 19
+0.002333 0.001000 275.862069 137.931034 18 20
+0.003667 0.003667 0.000000 0.000000 20 20
+0.003667 0.001000 413.793103 137.931034 20 20
+0.001000 0.003667 137.931034 413.793103 18 20
+0.001000 0.001000 0.000000 0.000000 20 20
+0.005000 0.005000 0.000000 137.931034 20 20
+0.002333 0.013000 137.931034 0.000000 20 20
+0.013000 0.001000 4000.000000 137.931034 20 20
+0.001000 0.011667 137.931034 1517.241379 13 14
+0.002333 0.002333 0.000000 0.000000 20 20
+0.001000 0.001000 413.793103 1241.379310 11 12
+0.010333 0.003667 0.000000 137.931034 20 20
+0.005000 0.005000 0.000000 0.000000 20 20
+0.006333 0.001000 965.517241 137.931034 9 12
+0.002333 0.001000 0.000000 0.000000 20 20
+0.003667 0.007667 137.931034 0.000000 20 20
+0.002333 0.013000 137.931034 4000.000000 20 20
+0.003667 0.003667 0.000000 0.000000 18 20
+0.013000 0.002333 4000.000000 137.931034 20 20
+0.002333 0.003667 0.000000 0.000000 20 20
+0.006333 0.003667 0.000000 137.931034 20 20
+0.001000 0.013000 137.931034 2896.551724 20 20
+0.011667 0.001000 3724.137931 137.931034 3 14
+0.005000 0.005000 137.931034 0.000000 20 20
+0.002333 0.003667 0.000000 0.000000 9 14
+0.002333 0.002333 0.000000 0.000000 18 19
+0.011667 0.002333 0.000000 137.931034 20 20
+0.002333 0.011667 137.931034 0.000000 20 20
+0.003667 0.005000 0.000000 0.000000 20 20
+0.001000 0.001000 137.931034 689.655172 18 20
+0.001000 0.001000 0.000000 0.000000 19 20
+0.009000 0.003667 0.000000 137.931034 20 20
+0.002333 0.013000 137.931034 137.931034 20 20
+0.013000 0.001000 0.000000 137.931034 20 20
+0.003667 0.006333 137.931034 0.000000 20 20
+0.001000 0.013000 137.931034 3586.206897 20 20
+0.013000 0.001000 3586.206897 137.931034 20 20
+0.005000 0.002333 0.000000 0.000000 20 20
+0.003667 0.009000 137.931034 0.000000 20 20
+0.002333 0.001000 0.000000 0.000000 14 20
+0.002333 0.001000 551.724138 137.931034 5 18
diff --git a/pykat/utilities/hermite.py b/pykat/utilities/hermite.py
new file mode 100644
index 0000000000000000000000000000000000000000..c9a5ea836f6f6b78abdeb53d8aca1eaa416c7acc
--- /dev/null
+++ b/pykat/utilities/hermite.py
@@ -0,0 +1,69 @@
+import numpy as np
+
+def hermite(n, x):
+
+ if n == 0:
+ return 1.0
+
+ elif n == 1:
+
+ return (2.0 * x)
+
+ elif n == 2:
+
+ return (4.0 * x * x - 2.0)
+
+ elif n == 3:
+
+ return (8.0 * x * x * x - 12.0 * x)
+
+
+ elif n == 4:
+ x_sq = x * x
+ return (16.0 * x_sq*x_sq - 48.0 * x_sq + 12.0)
+
+
+ elif n == 5:
+ x_sq = x * x
+ x_cb = x * x_sq
+ return (32.0 * x_cb * x_sq - 160.0 * x_cb + 120.0 * x)
+
+
+ elif n == 6:
+ x_sq = x * x
+ x_four = x_sq * x_sq
+ return (64.0 * x_four*x_sq - 480.0 * x_four + 720.0 * x_sq - 120.0)
+
+
+ elif n == 7:
+
+ x_sq = x*x
+ x_cb = x_sq*x
+ x_four = x_cb*x
+ return (128.0 * x_cb*x_four - 1344.0 * x_cb*x_sq + 3360.0 * x_cb - 1680.0 * x)
+
+
+ elif n == 8:
+
+ x_sq = x*x
+ x_four = x_sq*x_sq
+ x_six = x_four*x_sq
+ return (256.0 * x_four*x_four - 3584.0 * x_six + 13440.0 * x_four - 13440.0 * x_sq + 1680.0)
+
+ elif n == 9:
+ x_sq = x*x
+ x_cb = x_sq*x
+ x_four = x_cb*x
+ return (512.0 * x_cb*x_cb*x_cb - 9216.0 * x_four*x_cb + 48384.0 * x_cb*x_sq - 80640.0 * x_cb + 30240.0 * x)
+
+ elif n == 10:
+
+ x_sq = x*x
+ x_cb = x_sq*x
+ x_four = x_cb*x
+ x_five = x_four * x
+ return (1024.0 * x_five*x_five - 23040.0 * x_four*x_four + 161280.0 * x_cb*x_cb - 403200.0 * x_four + 302400.0 * x_sq - 30240.0)
+
+ else :
+
+ return (2 * x * hermite(n - 1, x) - 2 * (n - 1) * hermite(n - 2, x))
diff --git a/pykat/utilities/knm.py b/pykat/utilities/knm.py
new file mode 100644
index 0000000000000000000000000000000000000000..53ca9988938eae101f4d8e3f1f3638d1c4ee7805
--- /dev/null
+++ b/pykat/utilities/knm.py
@@ -0,0 +1,83 @@
+from itertools import combinations_with_replacement as combinations
+from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
+from pykat.exceptions import BasePyKatException
+
+import maps
+import os.path
+import numpy as np
+import pykat
+import collections
+import math
+
+def makeCouplingMatrix(max_order):
+ pass
+
+
+def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None):
+
+ if Axy == None:
+ Axy == np.ones((len(x), len(y)))
+
+ if q1y == None:
+ q1y = q1
+
+ if q2y == None:
+ q2y = q2
+
+ if len(mode_in) != 2 or len(mode_out) != 2:
+ raise BasePyKatException("Both mode in and out should be a container with modes [n m]")
+
+ Hg_in = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
+ Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
+
+ dx = abs(x[1] - x[0])
+ dy = abs(y[1] - y[0])
+
+ U1 = Hg_in.Unm(x,y)
+ U2 = Hg_out.Unm(x,y).conjugate()
+
+ return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
+
+
+def ROM_HG_knm(m, mode_in, mode_out, q1, q2, q1y=None, q2y=None):
+ pass
+
+
+def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
+
+ couplings = np.array(couplings)
+
+ a = couplings.size / 4.0
+
+ if int(a) - a != 0:
+ raise BasePyKatException("Iterator should be product of 4, each element of coupling array should be [n,m,n',m']")
+
+ if "phase" in surface_map.type:
+ Axy = np.exp(2j * 2 * math.pi / q1.wavelength * surface_map.data * surface_map.scaling)
+
+ x = surface_map.x
+ y = surface_map.y
+
+ K = np.zeros((couplings.size/4,), dtype=np.complex128)
+
+ it = np.nditer(couplings, flags=['refs_ok','f_index'])
+
+ i = 0
+
+ while not it.finished:
+ try:
+
+ mode_in = [int(it.next()), int(it.next())]
+ mode_out = [int(it.next()), int(it.next())]
+
+ if method == "riemann":
+ K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy)
+ else:
+ raise BasePyKatException("method value not accepted")
+
+ i +=1
+
+ except StopIteration:
+ break
+
+ return K.reshape(couplings.shape[:-1])
\ No newline at end of file
diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py
index 02bd5128006f3ba31e2bf7020138564b8428c9a1..e9f1ab23d19591cf67cfce8ee18f89af990fd370 100644
--- a/pykat/utilities/maps.py
+++ b/pykat/utilities/maps.py
@@ -1,6 +1,6 @@
+from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights
import numpy
-
class surfacemap:
def __init__(self, name, maptype, size, center, step_size, scaling, data=None):
@@ -11,6 +11,7 @@ class surfacemap:
self.step_size = step_size
self.scaling = scaling
self.data = data
+ self._rom_weights = None
def write_map(self, filename):
with open(filename,'w') as mapfile:
@@ -31,16 +32,32 @@ class surfacemap:
@property
def x(self):
- return self.step_size[0] * (numpy.array(range(0, self.data.shape[0]))- self.center[0])
+ return self.step_size[0] * (numpy.array(range(1, self.data.shape[0]+1)) - self.center[0])
@property
def y(self):
- return self.step_size[1] * (numpy.array(range(0, self.data.shape[1]))- self.center[1])
-
+ return self.step_size[1] * (numpy.array(range(1, self.data.shape[1]+1))- self.center[1])
+
+ @property
+ def offset(self):
+ return numpy.array(self.step_size)*(self.center - numpy.array(self.size)/2)
+
+ @property
+ def ROMWeights(self):
+ return self._rom_weights
+
+ def generateROMWeights(self):
+ b = makeReducedBasis(self.x[0:(len(self.x)/2)], offset=self.offset)
+ EI = makeEmpiricalInterpolant(b)
+ self._rom_weights = makeWeights(self, EI)
+
+ return self.ROMWeights
+
def plot(self, show=True, clabel=None):
import pylab
+ # 100 factor for scaling to cm
xrange = 100*self.x
yrange = 100*self.y
@@ -60,7 +77,9 @@ class surfacemap:
pylab.show()
return fig
-
+
+
+
def read_map(filename):
with open(filename, 'r') as f:
diff --git a/pykat/utilities/optics/gaussian_beams.py b/pykat/utilities/optics/gaussian_beams.py
index 79b67aa7b3477ea5ce31149c3bc344a3578c81dc..02ecfbfcfd28a7fbab6133b5f0230e39ba9d02eb 100644
--- a/pykat/utilities/optics/gaussian_beams.py
+++ b/pykat/utilities/optics/gaussian_beams.py
@@ -81,8 +81,58 @@ class gauss_param(object):
@property
def w(self):
return abs(self.__q)*math.sqrt(self.__lambda / (self.__nr * math.pi * self.__q.imag))
- #return self.w0 * math.sqrt(1 + (self.__q.real/self.__q.imag)**2)
+ def w(self, z=None, wavelength=None, nr=None, w0=None):
+
+ if z == None:
+ z = self.z
+ else:
+ z = np.array(z)
+
+ if wavelength == None:
+ wavelength = self.wavelength
+ else:
+ wavelength = np.array(wavelength)
+
+ if nr == None:
+ nr = self.nr
+ else:
+ nr = np.array(nr)
+
+ if w0 == None:
+ w0 = self.w0
+ else:
+ w0 = np.array(w0)
+
+ q = z + 1j * math.pi * w0 **2 / wavelength
+
+ return np.abs(q)*np.sqrt(wavelength / (nr * math.pi * q.imag))
+
+ def gouy(self, z=None, wavelength=None, nr=None, w0=None):
+ if z == None:
+ z = self.z
+ else:
+ z = np.array(z)
+
+ if wavelength == None:
+ wavelength = self.wavelength
+ else:
+ wavelength = np.array(wavelength)
+
+ if nr == None:
+ nr = self.nr
+ else:
+ nr = np.array(nr)
+
+ if w0 == None:
+ w0 = self.w0
+ else:
+ w0 = np.array(w0)
+
+ q = z + 1j * math.pi * w0 **2 / wavelength
+
+ return np.arctan2(q.real, q.imag)
+
@property
def w0(self):
return math.sqrt(self.__q.imag * self.__lambda / (self.__nr * math.pi))
@@ -94,6 +144,31 @@ class gauss_param(object):
else:
return float("inf")
+ def Rc(self, z=None, wavelength=None, nr=None, w0=None):
+ if z == None:
+ z = self.z
+ else:
+ z = np.array(z)
+
+ if wavelength == None:
+ wavelength = self.wavelength
+ else:
+ wavelength = np.array(wavelength)
+
+ if nr == None:
+ nr = self.nr
+ else:
+ nr = np.array(nr)
+
+ if w0 == None:
+ w0 = self.w0
+ else:
+ w0 = np.array(w0)
+
+ q = z + 1j * math.pi * w0 **2 / wavelength
+
+ return q.real * (1+ (q.imag/q.real)**2)
+
def conjugate(self):
return beam_param(self.__lambda, self.__nr, self.__q.conjugate())
@@ -145,6 +220,9 @@ class gauss_param(object):
return beam_param(self.__lambda, self.__nr, -self.__q)
def __eq__(self, q):
+ if q == None:
+ return False
+
return complex(q) == self.__q
@property
@@ -161,10 +239,8 @@ class gauss_param(object):
def reverse(self):
self.__q = -1.0 * self.__q.real + 1j * self.__q.imag
-
class beam_param(gauss_param):
pass
-
class HG_beam(object):
@@ -172,15 +248,15 @@ class HG_beam(object):
self._qx = copy.deepcopy(qx)
self._2pi_qrt = math.pow(2.0/math.pi, 0.25)
- if qy.__class__ == beam_param:
+ if qy == None:
self._qy = copy.deepcopy(qx)
else:
self._qy = copy.deepcopy(qy)
- self._n = n
- self._m = m
- self._hn = hermite(n)
- self._hm = hermite(m)
+ self._n = int(n)
+ self._m = int(m)
+ self._hn = hermite(self._n)
+ self._hm = hermite(self._m)
self._calc_constants()
@property
@@ -238,12 +314,12 @@ class HG_beam(object):
def _calc_constants(self):
self.__xpre_const = math.pow(2.0/math.pi, 0.25)
- self.__xpre_const *= math.sqrt(1.0/(2**self._n * math.factorial(self._n)))
+ self.__xpre_const *= math.sqrt(1.0/(self._qx.w*2**self._n * math.factorial(self._n)))
self.__xpre_const *= cmath.sqrt(1j*self._qx.imag / self._qx.q)
self.__xpre_const *= ((1j*self._qx.imag * self._qx.q.conjugate())/(-1j*self._qx.imag * self._qx.q)) ** ( self._n/2.0)
self.__ypre_const = math.pow(2.0/math.pi, 0.25)
- self.__ypre_const *= math.sqrt(1.0/(2**self._m * math.factorial(self._m)))
+ self.__ypre_const *= math.sqrt(1.0/(self._qy.w*2**self._m * math.factorial(self._m)))
self.__ypre_const *= cmath.sqrt(1j*self._qy.imag / self._qy.q)
self.__ypre_const *= ((1j*self._qy.imag * self._qy.q.conjugate())/(-1j*self._qy.imag * self._qy.q)) **(self._m/2.0)
@@ -262,8 +338,10 @@ class HG_beam(object):
def Um(self, y):
return self.__ypre_const * self._hm(self.__sqrt2_wyz * y) * np.exp(-0.5j * self.__ky * y*y * self.__invqy)
- def Unm(self, x, y):
- return self.Un(x) * self.Um(y)
+ def Unm(self, x, y):
+ _un = self.Un(x)
+ _um = self.Um(y)
+ return np.outer(_un, _um)
def plot(self, ndx=100, ndy=100, xscale=4, yscale=4):
import pylab
@@ -274,9 +352,7 @@ class HG_beam(object):
dx = xrange[1]-xrange[0]
dy = yrange[1]-yrange[0]
- xx, yy = np.meshgrid(xrange,yrange)
-
- data = self.Unm(xx, yy)
+ data = self.Unm(xrange,yrange)
fig = pylab.figure()
axes = pylab.imshow(np.abs(data), aspect=dx/dy, extent=[min(xrange),max(xrange),min(yrange),max(yrange)])
diff --git a/pykat/utilities/romhom.py b/pykat/utilities/romhom.py
new file mode 100644
index 0000000000000000000000000000000000000000..840e7aee915dd4ff155b76f6f72939e76aa9a44c
--- /dev/null
+++ b/pykat/utilities/romhom.py
@@ -0,0 +1,371 @@
+import maps
+import os.path
+import numpy as np
+import pykat
+import collections
+from itertools import combinations_with_replacement as combinations
+from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
+from scipy.linalg import inv
+from math import factorial
+from hermite import *
+import math
+
+EmpiricalInterpolant = collections.namedtuple('EmpiricalInterpolant', 'B nodes node_indices limits x')
+ReducedBasis = collections.namedtuple('ReducedBasis', 'RB limits x')
+ROMLimits = collections.namedtuple('ROMLimits', 'zmin zmax w0min w0max max_order')
+
+class ROMWeights:
+
+ def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, limits):
+ self.w_ij_Q1 = w_ij_Q1
+ self.w_ij_Q2 = w_ij_Q2
+ self.w_ij_Q3 = w_ij_Q3
+ self.w_ij_Q4 = w_ij_Q4
+ self.limits = limits
+
+ def writeToFile(self, filename):
+ f = open(filename, 'w+')
+
+ f.write(repr(self.limits) + "\n")
+
+ f.write(repr(self.w_ij_Q1.shape) + "\n")
+
+ for v in self.w_ij_Q1.flatten():
+ f.write("%s " % repr(complex(v))[1:-1])
+
+ f.write("\n")
+
+ f.write(repr(self.w_ij_Q2.shape) + "\n")
+
+ for v in self.w_ij_Q2.flatten():
+ f.write("%s " % repr(complex(v))[1:-1])
+
+ f.write("\n")
+
+ f.write(repr(self.w_ij_Q3.shape) + "\n")
+
+ for v in self.w_ij_Q3.flatten():
+ f.write("%s " % repr(complex(v))[1:-1])
+
+ f.write("\n")
+
+ f.write(repr(self.w_ij_Q4.shape) + "\n")
+
+ for v in self.w_ij_Q4.flatten():
+ f.write("%s " % repr(complex(v))[1:-1])
+
+ f.write("\n")
+
+ f.close()
+
+def combs(a, r):
+ """
+ Return successive r-length combinations of elements in the array a.
+ Should produce the same output as array(list(combinations(a, r))), but
+ faster.
+ """
+ a = np.asarray(a)
+ dt = np.dtype([('', a.dtype)]*r)
+ b = np.fromiter(combinations(a, r), dt)
+ return b.view(a.dtype).reshape(-1, r)
+
+def project_onto_basis(integration_weights, e, h, projections, proj_coefficients, idx):
+
+ for j in range(len(h)):
+ proj_coefficients[idx][j] = np.vdot(integration_weights* e[idx], h[j])
+ projections[j] += proj_coefficients[idx][j]*e[idx]
+
+ return projections
+
+def B_matrix(invV, e):
+ return np.inner(invV.T, e[0:(invV.shape[0])].T)
+
+def emp_interp(B_matrix, func, indices):
+ # B : RB matrix
+ assert B_matrix.shape[0] == len(indices)
+ interpolant = np.inner(func[indices].T, B_matrix.T)
+
+ return interpolant
+
+
+def w(w0, im_q, re_q):
+ return w0 * np.sqrt( 1 + (re_q / im_q)**2. )
+
+
+def u(re_q1, w0_1, n1, x):
+
+ im_q1 = np.pi*w0_1**2 / 1064e-9
+ q_z1 = re_q1 + 1j*im_q1
+
+ A_n1 = (2./np.pi)**(1./4.) * (1./((2.**n1)*factorial(n1)*w0_1))**(1./2.) * (im_q1 / q_z1)**(1./2.) * ( im_q1*q_z1.conjugate() / (-im_q1*q_z1) )**(n1/2.)
+
+ wz1 = w(w0_1, im_q1, re_q1)
+
+ return A_n1 * hermite(n1, np.sqrt(2.)*x / wz1) * np.exp(-1j*(2*math.pi/(1064e-9))* x**2 /(2.*q_z1))
+
+
+def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x):
+ return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x).conjugate()
+
+
+def makeReducedBasis(x, offset=np.array([0,0]), isModeMatched=True, tolerance = 1e-12, sigma = 1):
+
+ if isModeMatched:
+ greedypts = 'matched20.txt'
+ else:
+ greedypts = 'mismatched20.txt'
+
+ greedyfile = os.path.join(pykat.__path__[0],'utilities','greedypoints',greedypts)
+
+ limits = np.loadtxt(greedyfile, usecols=(1,))[:5]
+
+ romlimits = ROMLimits(w0min=limits[0], w0max=limits[1], zmin=limits[2], zmax=limits[3], max_order=limits[4])
+
+ w_min = limits[0]
+ w_max = limits[1]
+
+ re_q_min = limits[2]
+ re_q_max = limits[3]
+
+ max_order = int(limits[4])
+
+ indices = range(0, max_order+1)
+ nm_pairs = combs(indices, 2)
+
+ dx = x[1] - x[0]
+
+ params = np.loadtxt(greedyfile, skiprows=5)
+
+ TS_size = len(params) # training space of TS_size number of waveforms
+
+ #### allocate memory for training space ####
+
+ TS = np.zeros(TS_size*len(x), dtype = complex).reshape(TS_size, len(x)) # store training space in TS_size X len(x) array
+
+ IDx = 0 #TS index
+
+ for i in range(len(params)):
+
+ q1 = beam_param(w0=params[i][0], z=params[i][2])
+
+ if isModeMatched:
+ q2 = q1
+ n = int(params[i][2])
+ m = int(params[i][3])
+ w0_1 = params[i][0]
+ w0_2 = w0_1
+ re_q1 = params[i][1]
+ re_q2 = re_q1
+ else:
+ q2 = beam_param(w0=params[i][1], z=params[i][3])
+ n = int(params[i][4])
+ m = int(params[i][5])
+
+ TS[IDx] = u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, x)
+
+ # normalize
+ norm = np.sqrt(abs(np.vdot(dx * TS[IDx], TS[IDx])))
+
+ if norm != 0:
+ TS[IDx] /= norm
+ IDx += 1
+ else:
+ np.delete(TS, IDx)
+
+ #### Begin greedy: see Field et al. arXiv:1308.3565v2 ####
+
+ tolerance = 1e-12 # set maximum RB projection error
+
+ sigma = 1 # projection error at 0th iteration
+
+ RB_matrix = [TS[0]] # seed greedy algorithm (arbitrary)
+
+ proj_coefficients = np.zeros(TS_size*TS_size, dtype = complex).reshape(TS_size, TS_size)
+ projections = np.zeros(TS_size*len(x), dtype = complex).reshape(TS_size, len(x))
+
+ iter = 0
+
+ while sigma > tolerance:
+ #for k in range(TS_size-1):
+ # go through training set and find worst projection error
+ projections = project_onto_basis(dx, RB_matrix, TS, projections, proj_coefficients, iter)
+ residual = TS - projections
+ projection_errors = [np.vdot(dx* residual[i], residual[i]) for i in range(len(residual))]
+ sigma = abs(max(projection_errors))
+ index = np.argmax(projection_errors)
+
+ #Gram-Schmidt to get the next basis and normalize
+ next_basis = TS[index] - projections[index]
+ next_basis /= np.sqrt(abs(np.vdot(dx* next_basis, next_basis)))
+
+ RB_matrix.append(next_basis)
+
+ iter += 1
+
+ return ReducedBasis(RB=np.array(RB_matrix), limits=romlimits, x=x)
+
+def makeEmpiricalInterpolant(RB):
+
+ e_x = RB.RB
+ x = RB.x
+
+ node_indices = []
+ x_nodes = []
+ V = np.zeros((len(e_x), len(e_x)), dtype = complex)
+
+ # seed EIM algorithm
+ node_indices.append( int(np.argmax( np.abs(e_x[0]) )))
+ x_nodes.append(x[node_indices])
+
+ for i in range(1, len(e_x)):
+ #build empirical interpolant for e_iter
+ for j in range(len(node_indices)):
+ for k in range(len(node_indices)):
+ V[k][j] = e_x[j][node_indices[k]]
+
+ invV = inv(V[0:len(node_indices), 0:len(node_indices)])
+ B = B_matrix(invV, e_x)
+
+ interpolant = emp_interp(B, e_x[i], node_indices)
+ res = interpolant - e_x[i]
+
+ index = int(np.argmax(np.abs(res)))
+ node_indices.append(index)
+ x_nodes.append( x[index] )
+
+ # make B matrix with all the indices
+ for j in range(len(node_indices)):
+ for k in range(len(node_indices)):
+ V[k][j] = e_x[j][node_indices[k]]
+
+ invV = inv(V[0:len(node_indices), 0:len(node_indices)])
+ B = B_matrix(invV, e_x)
+
+ return EmpiricalInterpolant(B=np.array(B), nodes=np.array(x_nodes), node_indices=np.array(node_indices), limits=RB.limits, x=RB.x)
+
+
+def makeWeights(smap, EI):
+
+ # get full A_xy
+ A_xy = smap.data
+
+ hx = len(smap.x)/2
+ hy = len(smap.y)/2
+ fx = hx*2
+ fy = hy*2
+
+ # get A_xy in the four quadrants of the x-y plane
+ A_xy_Q1 = A_xy[0:hx][0:hy, hy:fy]
+ A_xy_Q2 = A_xy[hx:fx][0:hy, hy:fy]
+ A_xy_Q3 = A_xy[hx:fx][0:hy, 0:hy]
+ A_xy_Q4 = A_xy[0:hx][0:hy, 0:hy]
+
+ full_x = smap.x
+ full_y = smap.y
+
+ dx = full_x[1] - full_x[0]
+ dy = full_y[1] - full_y[0]
+
+ B = EI.B
+ nodes_nv = EI.nodes
+ nodes_pv = - EI.nodes
+
+ # make integration weights
+ w_ij_Q1 = np.zeros((len(B), len(B)), dtype = complex)
+ w_ij_Q2 = np.zeros((len(B), len(B)), dtype = complex)
+ w_ij_Q3 = np.zeros((len(B), len(B)), dtype = complex)
+ w_ij_Q4 = np.zeros((len(B), len(B)), dtype = complex)
+
+ for i in range(len(B)):
+ for j in range(len(B)):
+ B_ij_Q1 = np.outer(B[i], B[j][::-1])
+ w_ij_Q1[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q1, A_xy_Q1)
+
+ B_ij_Q2 = np.outer(B[i][::-1], B[j][::-1])
+ w_ij_Q2[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q2, A_xy_Q2)
+
+ B_ij_Q3 = np.outer(B[i][::-1], B[j])
+ w_ij_Q3[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q3, A_xy_Q3)
+
+ B_ij_Q4 = np.outer(B[i], B[j])
+ w_ij_Q4[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q4, A_xy_Q4)
+
+ return ROMWeights(w_ij_Q1=w_ij_Q1, w_ij_Q2=w_ij_Q2, w_ij_Q3=w_ij_Q3, w_ij_Q4=w_ij_Q4, limits=EI.limits)
+
+
+def ROMKnm(W, max_order, q1, q2, q1y=None, q2y=None):
+
+ if q1y == None:
+ q1y = q1
+
+ if q2y == None:
+ q2y = q2
+
+ # get symmetric and anti-symmetric w_ij's
+ w_ij_Q1Q3 = W.w_ij_Q1 + W.w_ij_Q3
+ w_ij_Q2Q4 = W.w_ij_Q2 + W.w_ij_Q4
+ w_ij_Q1Q2 = W.w_ij_Q1 + W.w_ij_Q2
+ w_ij_Q1Q4 = W.w_ij_Q1 + W.w_ij_Q4
+ w_ij_Q2Q3 = W.w_ij_Q2 + W.w_ij_Q3
+ w_ij_Q3Q4 = W.w_ij_Q3 + W.w_ij_Q4
+
+ w_ij_Q1Q2Q3Q4 = W.w_ij_Q1 + W.w_ij_Q3 + W.w_ij_Q2 + W.w_ij_Q4
+
+ num_fields = int((max_order + 1) * (max_order + 2) / 2);
+
+ K_ROQ = np.array((num_fields, num_fields))
+
+ for i in range(len(nm_pairs)):
+ for j in range(len(nm_pairs)):
+
+ # full quadrature
+ n = nm_pairs[i][0]
+ m = nm_pairs[i][1]
+ npr = nm_pairs[j][0]
+ mpr = nm_pairs[j][1]
+ u_xy = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, full_x), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, full_x))
+ k = dx*dy*np.einsum('ij,ij', u_xy, A_xy)
+
+ # ROQ
+ if nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 and nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0:
+ u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
+
+ elif nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1 and nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
+ u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
+
+ elif nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 and nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
+ u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
+
+ elif nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1 and nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0:
+ u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
+
+ elif nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 1 or nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 0:
+ u_x_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, n,m, nodes_nv)
+ u_y_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv)
+
+ if nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0 or nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
+ u_xy_nodes_Q1Q4 = np.outer(u_x_nodes, u_y_nodes)
+ u_xy_nodes_Q2Q3 = - u_xy_nodes_Q1Q4
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q1Q4, w_ij_Q1Q4) + np.einsum('ij,ij', u_xy_nodes_Q2Q3, w_ij_Q2Q3)
+ else:
+ u_xy_nodes_Q2Q4 = np.outer(u_x_nodes, u_y_nodes)
+ u_xy_nodes_Q1Q3 = - u_xy_nodes_Q2Q4
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q2Q4, w_ij_Q2Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q3, w_ij_Q1Q3)
+
+ elif nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 1 or nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 0:
+ u_x_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv)
+ u_y_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv)
+
+ if nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 or nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1:
+ u_xy_nodes_Q3Q4 = np.outer(u_x_nodes, u_y_nodes)
+ u_xy_nodes_Q1Q2 = - u_xy_nodes_Q3Q4
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q3Q4, w_ij_Q3Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q2, w_ij_Q1Q2)
+ else:
+ u_xy_nodes_Q2Q4 = np.outer(u_x_nodes, u_y_nodes)
+ u_xy_nodes_Q1Q3 = - u_xy_nodes_Q2Q4
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q2Q4, w_ij_Q2Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q3, w_ij_Q1Q3)
+