Log In | Get Help   
Home My Page Projects Code Snippets Project Openings Mareframe
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files
[mareframe] Diff of /trunk/gadget/grow.cc
[mareframe] / trunk / gadget / grow.cc Repository:
ViewVC logotype

Diff of /trunk/gadget/grow.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1, Mon Feb 10 17:09:07 2014 UTC revision 4, Thu Apr 30 17:32:47 2015 UTC
# Line 2  Line 2 
2  #include "grower.h"  #include "grower.h"
3  #include "agebandmatrix.h"  #include "agebandmatrix.h"
4  #include "gadget.h"  #include "gadget.h"
5    #include "doublematrix.h"
6    #include "matrix.h"
7    
8  /* Update the agebandmatrix to reflect the calculated growth  */  /* Update the agebandmatrix to reflect the calculated growth  */
9  /* Lgrowth contains the ratio of each length group that grows */  /* Lgrowth contains the ratio of each length group that grows */
# Line 9  Line 11 
11  /* the weight increase for each entry in Lgrowth              */  /* the weight increase for each entry in Lgrowth              */
12    
13  /* JMB changed to deal with very small weights a bit better   */  /* JMB changed to deal with very small weights a bit better   */
14  void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth, const DoubleMatrix& Wgrowth) {  void AgeBandMatrix::Grow(const Matrix& Lgrowth, const Matrix& Wgrowth) {
15    int i, lgrp, grow, maxlgrp;    int i, lgrp, grow, maxlgrp;
16    double num, wt, tmp;    double num, wt, tmp;
17    
18    maxlgrp = Lgrowth.Nrow();    maxlgrp = Lgrowth.Nrow();
19    
20    for (i = 0; i < nrow; i++) {    for (i = 0; i < nrow; i++) {
21                    int maxCol = v[i]->maxCol();
22                    int minCol = v[i]->minCol();
23      //the part that grows to or above the highest length group      //the part that grows to or above the highest length group
24      num = 0.0;      num = 0.0;
25      wt = 0.0;      wt = 0.0;
26      for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--) {  
27        for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++) {                  for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
28                            for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
29          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
30          num += tmp;          num += tmp;
31          wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);          wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);
32        }        }
33      }      }
34    
35      lgrp = v[i]->maxCol() - 1;                  lgrp = maxCol - 1;
36      if (isZero(num) || (wt < verysmall)) {      if (isZero(num) || (wt < verysmall)) {
37        (*v[i])[lgrp].setToZero();        (*v[i])[lgrp].setToZero();
38      } else {      } else {
# Line 35  Line 41 
41      }      }
42    
43      //the central diagonal part of the length division      //the central diagonal part of the length division
44      for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {                  for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
45        num = 0.0;        num = 0.0;
46        wt = 0.0;        wt = 0.0;
47        for (grow = 0; grow < maxlgrp; grow++) {        for (grow = 0; grow < maxlgrp; grow++) {
48          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
49          num += tmp;          num += tmp;
50          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);                                  wt += tmp
51                                                    * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
52        }        }
   
53        if (isZero(num) || (wt < verysmall)) {        if (isZero(num) || (wt < verysmall)) {
54          (*v[i])[lgrp].setToZero();          (*v[i])[lgrp].setToZero();
55        } else {        } else {
# Line 53  Line 59 
59      }      }
60    
61      //the lowest part of the length division      //the lowest part of the length division
62      for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {                  for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
63        num = 0.0;        num = 0.0;
64        wt = 0.0;        wt = 0.0;
65        for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++) {                          for (grow = 0; grow <= lgrp - minCol; grow++) {
66          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
67          num += tmp;          num += tmp;
68          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);                                  wt += tmp
69                                                    * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
70        }        }
   
71        if (isZero(num) || (wt < verysmall)) {        if (isZero(num) || (wt < verysmall)) {
72          (*v[i])[lgrp].setToZero();          (*v[i])[lgrp].setToZero();
73        } else {        } else {
# Line 72  Line 78 
78    }    }
79  }  }
80    
81    /********************************************************************************/
82    
83    ///* JMB changed to deal with very small weights a bit better   */
84    //void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth,
85    //              const DoubleMatrix& Wgrowth) {
86    //      int i, lgrp, grow, maxlgrp;
87    //      double num, wt, tmp;
88    //
89    //      maxlgrp = Lgrowth.Nrow();
90    //
91    //      for (i = 0; i < nrow; i++) {
92    //              int maxCol = v[i]->maxCol();
93    //              int minCol = v[i]->minCol();
94    //              //the part that grows to or above the highest length group
95    //              num = 0.0;
96    //              wt = 0.0;
97    //
98    //              for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
99    //                      for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
100    //                              tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
101    //                              num += tmp;
102    //                              wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);
103    //                      }
104    //              }
105    //
106    //              lgrp = maxCol - 1;
107    //              if (isZero(num) || (wt < verysmall)) {
108    //                      (*v[i])[lgrp].setToZero();
109    //              } else {
110    //                      (*v[i])[lgrp].W = wt / num;
111    //                      (*v[i])[lgrp].N = num;
112    //              }
113    //
114    //              //the central diagonal part of the length division
115    //              for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
116    //                      num = 0.0;
117    //                      wt = 0.0;
118    //                      for (grow = 0; grow < maxlgrp; grow++) {
119    //                              tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
120    //                              num += tmp;
121    //                              wt += tmp
122    //                                              * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
123    //                      }
124    //
125    //                      if (isZero(num) || (wt < verysmall)) {
126    //                              (*v[i])[lgrp].setToZero();
127    //                      } else {
128    //                              (*v[i])[lgrp].W = wt / num;
129    //                              (*v[i])[lgrp].N = num;
130    //                      }
131    //              }
132    //
133    //              //the lowest part of the length division
134    //              for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
135    //                      num = 0.0;
136    //                      wt = 0.0;
137    //                      for (grow = 0; grow <= lgrp - minCol; grow++) {
138    //                              tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
139    //                              num += tmp;
140    //                              wt += tmp
141    //                                              * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
142    //                      }
143    //
144    //                      if (isZero(num) || (wt < verysmall)) {
145    //                              (*v[i])[lgrp].setToZero();
146    //                      } else {
147    //                              (*v[i])[lgrp].W = wt / num;
148    //                              (*v[i])[lgrp].N = num;
149    //                      }
150    //              }
151    //      }
152    //}
153  //Same program with certain num of fish made mature.  //Same program with certain num of fish made mature.
154  void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth, const DoubleMatrix& Wgrowth, Maturity* const Mat, int area) {  void AgeBandMatrix::Grow(const Matrix& Lgrowth, const Matrix& Wgrowth,
155                    Maturity* const Mat, int area) {
156    
157    int i, lgrp, grow, maxlgrp, age;    int i, lgrp, grow, maxlgrp, age;
158    double num, wt, matnum, tmp, ratio;    double num, wt, matnum, tmp, ratio;
159    
160    maxlgrp = Lgrowth.Nrow();    maxlgrp = Lgrowth.Nrow();
161    for (i = 0; i < nrow; i++) {    for (i = 0; i < nrow; i++) {
162                    int maxCol = v[i]->maxCol();
163                    int minCol = v[i]->minCol();
164      age = i + minage;      age = i + minage;
165      num = 0.0;      num = 0.0;
166      wt = 0.0;      wt = 0.0;
167      matnum = 0.0;      matnum = 0.0;
168      for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--) {                  for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
169        for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++) {                          for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
170          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp].W);          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp].W);
171          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
172          matnum += (tmp * ratio);          matnum += (tmp * ratio);
# Line 94  Line 175 
175        }        }
176      }      }
177    
178      lgrp = v[i]->maxCol() - 1;                  lgrp = maxCol - 1;
179      if (isZero(num) || (wt < verysmall)) {      if (isZero(num) || (wt < verysmall)) {
180        //no fish grow to this length cell        //no fish grow to this length cell
181        (*v[i])[lgrp].setToZero();        (*v[i])[lgrp].setToZero();
# Line 114  Line 195 
195        Mat->storeMatureStock(area, age, lgrp, matnum, wt / num);        Mat->storeMatureStock(area, age, lgrp, matnum, wt / num);
196      }      }
197    
198      for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {                  for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
199        num = 0.0;        num = 0.0;
200        wt = 0.0;        wt = 0.0;
201        matnum = 0.0;        matnum = 0.0;
202        for (grow = 0; grow < maxlgrp; grow++) {        for (grow = 0; grow < maxlgrp; grow++) {
203          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp - grow].W);                                  ratio = Mat->calcMaturation(age, lgrp, grow,
204                                                    (*v[i])[lgrp - grow].W);
205          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
206          matnum += (tmp * ratio);          matnum += (tmp * ratio);
207          num += tmp;          num += tmp;
208          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);                                  wt += tmp
209                                                    * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
210        }        }
211    
212        if (isZero(num) || (wt < verysmall)) {        if (isZero(num) || (wt < verysmall)) {
# Line 146  Line 229 
229        }        }
230      }      }
231    
232      for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {                  for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
233        num = 0.0;        num = 0.0;
234        wt = 0.0;        wt = 0.0;
235        matnum = 0.0;        matnum = 0.0;
236        for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++) {                          for (grow = 0; grow <= lgrp - minCol; grow++) {
237          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp - grow].W);                                  ratio = Mat->calcMaturation(age, lgrp, grow,
238                                                    (*v[i])[lgrp - grow].W);
239          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
240          matnum += (tmp * ratio);          matnum += (tmp * ratio);
241          num += tmp;          num += tmp;
242          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);                                  wt += tmp
243                                                    * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
244        }        }
245    
246        if (isZero(num) || (wt < verysmall)) {        if (isZero(num) || (wt < verysmall)) {
# Line 181  Line 266 
266  }  }
267    
268  //fleksibest formulation - weight read in from file (should be positive)  //fleksibest formulation - weight read in from file (should be positive)
269  void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth, const DoubleVector& Weight) {  void AgeBandMatrix::Grow(const Matrix& Lgrowth, const DoubleVector& Weight) {
270    int i, lgrp, grow, maxlgrp;    int i, lgrp, grow, maxlgrp;
271    double num;    double num;
272    
273    maxlgrp = Lgrowth.Nrow();    maxlgrp = Lgrowth.Nrow();
274    
275    for (i = 0; i < nrow; i++) {    for (i = 0; i < nrow; i++) {
276                    int maxCol = v[i]->maxCol();
277                    int minCol = v[i]->minCol();
278      num = 0.0;      num = 0.0;
279      for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--)                  for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--)
280        for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++)                          for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++)
281          num += (Lgrowth[grow][lgrp] * (*v[i])[lgrp].N);          num += (Lgrowth[grow][lgrp] * (*v[i])[lgrp].N);
282    
283      lgrp = v[i]->maxCol() - 1;                  lgrp = maxCol - 1;
284      if (isZero(num)) {      if (isZero(num)) {
285        (*v[i])[lgrp].setToZero();        (*v[i])[lgrp].setToZero();
286      } else {      } else {
# Line 200  Line 288 
288        (*v[i])[lgrp].W = Weight[lgrp];        (*v[i])[lgrp].W = Weight[lgrp];
289      }      }
290    
291      for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {                  for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
292        num = 0.0;        num = 0.0;
293        for (grow = 0; grow < maxlgrp; grow++)        for (grow = 0; grow < maxlgrp; grow++)
294          num += (Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N);          num += (Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N);
# Line 213  Line 301 
301        }        }
302      }      }
303    
304      for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {                  for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
305        num = 0.0;        num = 0.0;
306        for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++)                          for (grow = 0; grow <= lgrp - minCol; grow++)
307          num += (Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N);          num += (Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N);
308    
309        if (isZero(num)) {        if (isZero(num)) {
# Line 230  Line 318 
318    
319  //fleksibest formulation - weight read in from file (should be positive)  //fleksibest formulation - weight read in from file (should be positive)
320  //Same program with certain num of fish made mature.  //Same program with certain num of fish made mature.
321  void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth, const DoubleVector& Weight, Maturity* const Mat, int area) {  void AgeBandMatrix::Grow(const Matrix& Lgrowth, const DoubleVector& Weight,
322                    Maturity* const Mat, int area) {
323    
324    int i, lgrp, grow, maxlgrp, age;    int i, lgrp, grow, maxlgrp, age;
325    double num, matnum, tmp, ratio;    double num, matnum, tmp, ratio;
326    
327    maxlgrp = Lgrowth.Nrow();    maxlgrp = Lgrowth.Nrow();
328    
329    for (i = 0; i < nrow; i++) {    for (i = 0; i < nrow; i++) {
330                    int maxCol = v[i]->maxCol();
331                    int minCol = v[i]->minCol();
332      age = i + minage;      age = i + minage;
333      num = 0.0;      num = 0.0;
334      matnum = 0.0;      matnum = 0.0;
335      for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--) {                  for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
336        for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++) {                          for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
337          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp].W);          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp].W);
338          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
339          matnum += (tmp * ratio);          matnum += (tmp * ratio);
# Line 249  Line 341 
341        }        }
342      }      }
343    
344      lgrp = v[i]->maxCol() - 1;                  lgrp = maxCol - 1;
345      if (isZero(num)) {      if (isZero(num)) {
346        //no fish grow to this length cell        //no fish grow to this length cell
347        (*v[i])[lgrp].setToZero();        (*v[i])[lgrp].setToZero();
# Line 269  Line 361 
361        Mat->storeMatureStock(area, age, lgrp, matnum, Weight[lgrp]);        Mat->storeMatureStock(area, age, lgrp, matnum, Weight[lgrp]);
362      }      }
363    
364      for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {                  for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
365        num = 0.0;        num = 0.0;
366        matnum = 0.0;        matnum = 0.0;
367        for (grow = 0; grow < maxlgrp; grow++) {        for (grow = 0; grow < maxlgrp; grow++) {
368          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp - grow].W);                                  ratio = Mat->calcMaturation(age, lgrp, grow,
369                                                    (*v[i])[lgrp - grow].W);
370          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
371          matnum += (tmp * ratio);          matnum += (tmp * ratio);
372          num += tmp;          num += tmp;
# Line 299  Line 392 
392        }        }
393      }      }
394    
395      for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {                  for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
396        num = 0.0;        num = 0.0;
397        matnum = 0.0;        matnum = 0.0;
398        for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++) {                          for (grow = 0; grow <= lgrp - minCol; grow++) {
399          ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp - grow].W);                                  ratio = Mat->calcMaturation(age, lgrp, grow,
400                                                    (*v[i])[lgrp - grow].W);
401          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
402          matnum += (tmp * ratio);          matnum += (tmp * ratio);
403          num += tmp;          num += tmp;
Legend:
Removed from v.1  
changed lines
  Added in v.4
root@forge.cesga.es
 ViewVC Help
Powered by ViewVC 1.0.0  

Powered By FusionForge