Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 1939:2730d17df021

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added the implementation of the 'Bezugslinienverfahren'. flys-artifacts/trunk@3320 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 25 Nov 2011 18:52:11 +0000
parents 1d991c91285b
children 2898b1ff6013
comparison
equal deleted inserted replaced
1938:1d991c91285b 1939:2730d17df021
87 this.weight = weight; 87 this.weight = weight;
88 return this; 88 return this;
89 } 89 }
90 90
91 } // class Position 91 } // class Position
92
93 public static final class SplineFunction {
94
95 public PolynomialSplineFunction spline;
96 public double [] splineQs;
97 public double [] splineWs;
98
99 public SplineFunction(
100 PolynomialSplineFunction spline,
101 double [] splineQs,
102 double [] splineWs
103 ) {
104 this.spline = spline;
105 this.splineQs = splineQs;
106 this.splineWs = splineWs;
107 }
108
109 public double [][] sample(
110 int numSamples,
111 double km,
112 Calculation errors
113 ) {
114 double minQ = getQMin();
115 double maxQ = getQMax();
116
117 double [] outWs = new double[numSamples];
118 double [] outQs = new double[numSamples];
119
120 Arrays.fill(outWs, Double.NaN);
121 Arrays.fill(outQs, Double.NaN);
122
123 double stepWidth = (maxQ - minQ)/numSamples;
124
125 try {
126 double q = minQ;
127 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
128 outWs[i] = spline.value(outQs[i] = q);
129 }
130 }
131 catch (ArgumentOutsideDomainException aode) {
132 if (errors != null) {
133 // TODO: I18N
134 errors.addProblem(km, "spline interpolation failed");
135 }
136 log.error("spline interpolation failed.", aode);
137 }
138
139 return new double [][] { outWs, outQs };
140 }
141
142 public double getQMin() {
143 return Math.min(splineQs[0], splineQs[splineQs.length-1]);
144 }
145
146 public double getQMax() {
147 return Math.max(splineQs[0], splineQs[splineQs.length-1]);
148 }
149
150 /** Constructs a continues index between the columns to Qs. */
151 public PolynomialSplineFunction createIndexQRelation() {
152
153 double [] indices = new double[splineQs.length];
154 for (int i = 0; i < indices.length; ++i) {
155 indices[i] = i;
156 }
157
158 try {
159 SplineInterpolator interpolator = new SplineInterpolator();
160 return interpolator.interpolate(indices, splineQs);
161 }
162 catch (MathIllegalArgumentException miae) {
163 // Ignore me!
164 }
165 return null;
166 }
167 } // class SplineFunction
92 168
93 /** 169 /**
94 * A row, typically a position where measurements were taken. 170 * A row, typically a position where measurements were taken.
95 */ 171 */
96 public static final class Row 172 public static final class Row
164 ows[i] = Double.NaN; 240 ows[i] = Double.NaN;
165 } 241 }
166 } 242 }
167 } 243 }
168 244
169 public static final class SplineFunction {
170
171 public PolynomialSplineFunction spline;
172 public double [] splineQs;
173 public double [] splineWs;
174
175 public SplineFunction(
176 PolynomialSplineFunction spline,
177 double [] splineQs,
178 double [] splineWs
179 ) {
180 this.spline = spline;
181 this.splineQs = splineQs;
182 this.splineWs = splineWs;
183 }
184
185 public double [][] sample(
186 int numSamples,
187 double km,
188 Calculation errors
189 ) {
190 double q1 = splineQs[0], qN = splineQs[splineQs.length-1];
191 double minQ = Math.min(q1, qN);
192 double maxQ = Math.max(q1, qN);
193
194 double [] outWs = new double[numSamples];
195 double [] outQs = new double[numSamples];
196
197 Arrays.fill(outWs, Double.NaN);
198 Arrays.fill(outQs, Double.NaN);
199
200 double stepWidth = (maxQ - minQ)/numSamples;
201
202 try {
203 double q = minQ;
204 for (int i = 0; i < outWs.length; ++i, q += stepWidth) {
205 outWs[i] = spline.value(outQs[i] = q);
206 }
207 }
208 catch (ArgumentOutsideDomainException aode) {
209 if (errors != null) {
210 // TODO: I18N
211 errors.addProblem(km, "spline interpolation failed");
212 }
213 log.error("spline interpolation failed.", aode);
214 }
215
216 return new double [][] { outWs, outQs };
217 }
218 } // class SplineFunction
219 245
220 public SplineFunction createSpline( 246 public SplineFunction createSpline(
221 WstValueTable table, 247 WstValueTable table,
222 Calculation errors 248 Calculation errors
223 ) { 249 ) {
819 return new double [][] {qs, ws}; 845 return new double [][] {qs, ws};
820 } 846 }
821 847
822 public double [] findQsForW(double km, double w) { 848 public double [] findQsForW(double km, double w) {
823 849
824
825 int rowIndex = Collections.binarySearch(rows, new Row(km)); 850 int rowIndex = Collections.binarySearch(rows, new Row(km));
826 851
827 if (rowIndex >= 0) { 852 if (rowIndex >= 0) {
828 return rows.get(rowIndex).findQsForW(w, this); 853 return rows.get(rowIndex).findQsForW(w, this);
829 } 854 }
830 855
831 rowIndex = -rowIndex - 1; 856 rowIndex = -rowIndex - 1;
832 857
833 if (rowIndex < 1 || rowIndex >= rows.size()) { 858 if (rowIndex < 1 || rowIndex >= rows.size()) {
834 // do not extrapolate 859 // Do not extrapolate.
835 return new double[0]; 860 return new double[0];
836 } 861 }
837 862
838 863 // Needs bilinear interpolation.
839 // Needs bilinear interpolation
840 Row r1 = rows.get(rowIndex-1); 864 Row r1 = rows.get(rowIndex-1);
841 Row r2 = rows.get(rowIndex); 865 Row r2 = rows.get(rowIndex);
842 866
843 return r1.findQsForW(r2, w, km, this); 867 return r1.findQsForW(r2, w, km, this);
868 }
869
870 protected SplineFunction createSpline(double km, Calculation errors) {
871
872 int rowIndex = Collections.binarySearch(rows, new Row(km));
873
874 if (rowIndex >= 0) {
875 SplineFunction sf = rows.get(rowIndex).createSpline(this, errors);
876 if (sf == null && errors != null) {
877 // TODO: I18N
878 errors.addProblem(km, "cannot create w/q relation");
879 }
880 return sf;
881 }
882
883 rowIndex = -rowIndex - 1;
884
885 if (rowIndex < 1 || rowIndex >= rows.size()) {
886 // Do not extrapolate.
887 if (errors != null) {
888 // TODO: I18N
889 errors.addProblem(km, "km not found");
890 }
891 return null;
892 }
893
894 // Needs bilinear interpolation.
895 Row r1 = rows.get(rowIndex-1);
896 Row r2 = rows.get(rowIndex);
897
898 SplineFunction sf = r1.createSpline(r2, km, this, errors);
899 if (sf == null && errors != null) {
900 // TODO: I18N
901 errors.addProblem(km, "cannot create w/q relation");
902 }
903
904 return sf;
905 }
906
907 /** 'Bezugslinienverfahren' */
908 public double [][] relateWs(
909 double km1,
910 double km2,
911 int numSamples,
912 Calculation errors
913 ) {
914 SplineFunction sf1 = createSpline(km1, errors);
915 if (sf1 == null) {
916 return new double[2][0];
917 }
918
919 SplineFunction sf2 = createSpline(km1, errors);
920 if (sf2 == null) {
921 return new double[2][0];
922 }
923
924 PolynomialSplineFunction iQ1 = sf1.createIndexQRelation();
925 if (iQ1 == null) {
926 if (errors != null) {
927 // TODO: I18N
928 errors.addProblem(km1, "cannot create index/q relation");
929 }
930 return new double[2][0];
931 }
932
933 PolynomialSplineFunction iQ2 = sf1.createIndexQRelation();
934 if (iQ2 == null) {
935 if (errors != null) {
936 // TODO: I18N
937 errors.addProblem(km2, "cannot create index/q relation");
938 }
939 return new double[2][0];
940 }
941
942 int N = sf1.splineQs.length;
943 double stepWidth = N/(double)numSamples;
944
945 PolynomialSplineFunction qW1 = sf1.spline;
946 PolynomialSplineFunction qW2 = sf2.spline;
947
948 double [] ws1 = new double[numSamples];
949 double [] ws2 = new double[numSamples];
950
951 Arrays.fill(ws1, Double.NaN);
952 Arrays.fill(ws2, Double.NaN);
953
954 boolean hadErrors = false;
955
956 double p = 0d;
957 for (int i = 0; i < numSamples; ++i, p += stepWidth) {
958 try {
959 double q1 = iQ1.value(p);
960 double w1 = qW1.value(q1);
961 double q2 = iQ2.value(p);
962 double w2 = qW2.value(q2);
963 ws1[i] = w1;
964 ws2[i] = w2;
965 }
966 catch (ArgumentOutsideDomainException aode) {
967 if (!hadErrors) {
968 // XXX: I'm not sure if this really can happen
969 // and if we should report this more than once.
970 hadErrors = true;
971 if (errors != null) {
972 // TODO: I18N
973 log.debug("W~W failed", aode);
974 errors.addProblem("W~W failed");
975 }
976 }
977 }
978 }
979
980 return new double [][] { ws1, ws2 };
844 } 981 }
845 982
846 public QPosition getQPosition(double km, double q) { 983 public QPosition getQPosition(double km, double q) {
847 return getQPosition(km, q, new QPosition()); 984 return getQPosition(km, q, new QPosition());
848 } 985 }
http://dive4elements.wald.intevation.org