Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/exports/fixings/FixATWriter.java @ 3641:1dc904370a64

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed issue687: Rewrote AT export in FixA flys-artifacts/trunk@5361 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 04 Sep 2012 17:32:37 +0000
parents 0c8a6145098b
children 83c0735092a3
comparison
equal deleted inserted replaced
3640:1d2856de489d 3641:1dc904370a64
12 12
13 import java.io.IOException; 13 import java.io.IOException;
14 import java.io.PrintWriter; 14 import java.io.PrintWriter;
15 import java.io.Writer; 15 import java.io.Writer;
16 16
17 import java.util.ArrayDeque;
18 import java.util.Collection;
19 import java.util.Locale; 17 import java.util.Locale;
20 18
21 import org.apache.log4j.Logger; 19 import org.apache.log4j.Logger;
22 20
23 public class FixATWriter 21 public class FixATWriter
30 public static final String I18N_HEADER_DEFAULT = 28 public static final String I18N_HEADER_DEFAULT =
31 "Exported fixings discharge curve for {0} {0}-km: {1}"; 29 "Exported fixings discharge curve for {0} {0}-km: {1}";
32 30
33 public static final String [] Q_MAX_COLUMN = new String [] { "max_q" }; 31 public static final String [] Q_MAX_COLUMN = new String [] { "max_q" };
34 32
35 public static class WQ { 33 private static final int MAX_ITERATIONS = 10000;
36 34 private static final double EPSILON = 1e-8;
37 protected double w;
38 protected double q;
39
40 public WQ() {
41 }
42
43 public WQ(double w, double q) {
44 this.w = w;
45 this.q = q;
46 }
47 } // class WQ
48 35
49 protected Function function; 36 protected Function function;
50 protected Parameters parameters; 37 protected Parameters parameters;
51 38
52 public FixATWriter() { 39 public FixATWriter() {
92 Function inverse = function.getInverse(); 79 Function inverse = function.getInverse();
93 80
94 de.intevation.flys.artifacts.math.Function invInst = 81 de.intevation.flys.artifacts.math.Function invInst =
95 inverse.instantiate(coeffs); 82 inverse.instantiate(coeffs);
96 83
97 /* TODO: Do some fancy root finding stuff with the 84 double wMin = minW(invInst, wMax, qMax[0]);
98 * first derivative of the inverse function to find
99 * the minimum of function and directly export
100 * the longest strict monotonous ascending run of
101 * ws of the inverse in range [0, wMax].
102 *
103 * To simplify the situation we iterate in steps
104 * of 10cm over the range and export the longest
105 * run.
106 */
107 85
108 ArrayDeque<WQ> longest = new ArrayDeque<WQ>(); 86 double wMinCM = wMin * 100d;
109 ArrayDeque<WQ> current = new ArrayDeque<WQ>(); 87 double wMaxCM = wMax * 100d;
110 88
111 for (double w = 0d; w <= wMax; w += 0.1d) { 89 int wRow = ((int)wMinCM / 10) * 10;
112 double q = invInst.value(w);
113 90
114 if (Double.isNaN(q)) { 91 if ((wMinCM - (int)wMinCM) > 0d) {
115 log.debug("Eval of inverse for " + w + " failed."); 92 wMinCM = (int)wMinCM + 1d;
116 continue; 93 }
117 }
118 94
119 WQ wq = new WQ(w, q); 95 double w = wMinCM / 100.0;
120 96
121 if (current.isEmpty() || current.getLast().q < q) { 97 int wcm = ((int)wMinCM) % 10;
122 current.add(wq); 98
123 } 99 if (log.isDebugEnabled()) {
124 else { 100 log.debug("wMinCM: " + wMinCM);
125 if (current.size() >= longest.size()) { 101 log.debug("wMaxCM: " + wMaxCM);
126 longest = current; 102 log.debug("wcm: " + wcm);
127 } 103 }
128 current = new ArrayDeque<WQ>(); 104
129 current.add(wq); 105 out.printf(Locale.US, "%8d", wRow);
106
107 if (wcm > 0) {
108 int rest = 10 - wcm;
109 while (rest-- > 0) {
110 out.print(ATWriter.EMPTY);
130 } 111 }
131 } 112 }
132 113
133 if (current.size() >= longest.size()) { 114 for (;;) {
134 longest = current; 115 while (wcm++ < 10) {
116 if (w > wMax) {
117 break;
118 }
119 double q = invInst.value(w);
120 if (Double.isNaN(w)) {
121 out.print(ATWriter.EMPTY);
122 }
123 else {
124 ATWriter.printQ(out, q);
125 }
126 w += 0.01d;
127 }
128 out.println();
129 if (w > wMax) {
130 break;
131 }
132 out.printf(Locale.US, "%8d", wRow += 10);
133 wcm = 0;
135 } 134 }
136 135
137 printWQs(out, longest);
138 out.flush(); 136 out.flush();
139 }
140
141 protected void printWQs(PrintWriter out, Collection<WQ> wqs) {
142 int lastColumn = 10;
143 for (WQ wq: wqs) {
144 int column = (int)(wq.w * 10d) % 10;
145
146 if (lastColumn > column) {
147 for (; lastColumn < 10; ++lastColumn) {
148 out.print(ATWriter.EMPTY);
149 }
150 out.println();
151 out.printf(Locale.US, "%8d", (int)Math.round(wq.w*100.0));
152 lastColumn = 0;
153 }
154 for (;lastColumn < column-1; ++lastColumn) {
155 out.print(ATWriter.EMPTY);
156 }
157 ATWriter.printQ(out, wq.q);
158
159 lastColumn = column;
160 }
161 out.println();
162 } 137 }
163 138
164 protected void printHeader( 139 protected void printHeader(
165 PrintWriter out, 140 PrintWriter out,
166 CallMeta meta, 141 CallMeta meta,
171 meta, 146 meta,
172 I18N_HEADER_KEY, 147 I18N_HEADER_KEY,
173 I18N_HEADER_DEFAULT, 148 I18N_HEADER_DEFAULT,
174 river, km)); 149 river, km));
175 } 150 }
151
152 private static double minW(
153 de.intevation.flys.artifacts.math.Function function,
154 double maxW,
155 double maxQ
156 ) {
157 double stepWidth = 10d;
158
159 double lastW = maxW;
160 double lastQ = maxQ;
161
162 for (int i = 0; i < MAX_ITERATIONS; ++i) {
163 double w = lastW - stepWidth;
164 double q = function.value(w);
165
166 if (Double.isNaN(q) || q > lastQ || q < 0d) {
167 if (stepWidth < EPSILON) {
168 break;
169 }
170 stepWidth *= 0.5d;
171 continue;
172 }
173
174 lastW = w;
175 lastQ = q;
176 }
177
178 return lastW;
179 }
176 } 180 }
177 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 : 181 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org