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 3, Thu Apr 30 12:52:16 2015 UTC revision 4, Thu Apr 30 17:32:47 2015 UTC
# Line 10  Line 10 
10  /* by a certain number of length groups, and Wgrowth contains */  /* by a certain number of length groups, and Wgrowth contains */
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   */
14    void AgeBandMatrix::Grow(const Matrix& Lgrowth, const Matrix& 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;
 //    ofstream outfile;  
 //      outfile.open("aaa", ios::out);  
 //    _Lgrowth.Print(outfile);  
 //    int a;  
 //    cout << "1-----------------" << endl;  
 //    cin >> a;  
   
17    
18      maxlgrp = Lgrowth.Nrow();      maxlgrp = Lgrowth.Nrow();
19    
   
   //   double** lGrowth = ToMatrix(Lgrowth, maxlgrp, aux);  
   //   double** wGrowth = ToMatrix(Wgrowth, maxlgrp, aux);  
   
   
20      for (i = 0; i < nrow; i++) {      for (i = 0; i < nrow; i++) {
21          int maxCol = v[i]->maxCol();          int maxCol = v[i]->maxCol();
22          int minCol = v[i]->minCol();          int minCol = v[i]->minCol();
# Line 43  Line 31 
31            wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);            wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);
32          }          }
33        }        }
 //      cout << "wt1:" << wt << endl;  
34    
35        lgrp = maxCol - 1;        lgrp = maxCol - 1;
36        if (isZero(num) || (wt < verysmall)) {        if (isZero(num) || (wt < verysmall)) {
# Line 53  Line 40 
40          (*v[i])[lgrp].N = num;          (*v[i])[lgrp].N = num;
41        }        }
42    
   
43        //the central diagonal part of the length division        //the central diagonal part of the length division
44        for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {        for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
45          num = 0.0;          num = 0.0;
# Line 61  Line 47 
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          }          }
 //        cout << "wt2:" << wt << endl;  
53          if (isZero(num) || (wt < verysmall)) {          if (isZero(num) || (wt < verysmall)) {
54            (*v[i])[lgrp].setToZero();            (*v[i])[lgrp].setToZero();
55          } else {          } else {
# Line 72  Line 58 
58          }          }
59        }        }
60    
   
61        //the lowest part of the length division        //the lowest part of the length division
62        for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {        for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
63          num = 0.0;          num = 0.0;
# Line 80  Line 65 
65          for (grow = 0; grow <= lgrp - 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 90  Line 75 
75            (*v[i])[lgrp].N = num;            (*v[i])[lgrp].N = num;
76          }          }
77        }        }
 //      cout << "wt3:" << wt << endl;  
78      }      }
79    }    }
80    
   
   
   
   
   
   
   
   
81    /********************************************************************************/    /********************************************************************************/
82    
83  /* JMB changed to deal with very small weights a bit better   */  ///* JMB changed to deal with very small weights a bit better   */
84  void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth, const DoubleMatrix& Wgrowth) {  //void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth,
85    int i, lgrp, grow, maxlgrp;  //              const DoubleMatrix& Wgrowth) {
86    double num, wt, tmp;  //      int i, lgrp, grow, maxlgrp;
87    //      double num, wt, tmp;
88    maxlgrp = Lgrowth.Nrow();  //
89    //      maxlgrp = Lgrowth.Nrow();
90  //   double** lGrowth = ToMatrix(Lgrowth, maxlgrp, aux);  //
91  //   double** wGrowth = ToMatrix(Wgrowth, maxlgrp, aux);  //      for (i = 0; i < nrow; i++) {
92    //              int maxCol = v[i]->maxCol();
93    //              int minCol = v[i]->minCol();
94    for (i = 0; i < nrow; i++) {  //              //the part that grows to or above the highest length group
95          int maxCol = v[i]->maxCol();  //              num = 0.0;
96          int minCol = v[i]->minCol();  //              wt = 0.0;
97      //the part that grows to or above the highest length group  //
98      num = 0.0;  //              for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
99      wt = 0.0;  //                      for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
100    //                              tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
101      for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {  //                              num += tmp;
102        for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {  //                              wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);
103          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;  //                      }
104          num += tmp;  //              }
105          wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);  //
106        }  //              lgrp = maxCol - 1;
107      }  //              if (isZero(num) || (wt < verysmall)) {
108    //                      (*v[i])[lgrp].setToZero();
109      lgrp = maxCol - 1;  //              } else {
110      if (isZero(num) || (wt < verysmall)) {  //                      (*v[i])[lgrp].W = wt / num;
111        (*v[i])[lgrp].setToZero();  //                      (*v[i])[lgrp].N = num;
112      } else {  //              }
113        (*v[i])[lgrp].W = wt / num;  //
114        (*v[i])[lgrp].N = num;  //              //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      //the central diagonal part of the length division  //                      for (grow = 0; grow < maxlgrp; grow++) {
119      for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {  //                              tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
120        num = 0.0;  //                              num += tmp;
121        wt = 0.0;  //                              wt += tmp
122        for (grow = 0; grow < maxlgrp; grow++) {  //                                              * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
123          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;  //                      }
124          num += tmp;  //
125          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);  //                      if (isZero(num) || (wt < verysmall)) {
126        }  //                              (*v[i])[lgrp].setToZero();
127    //                      } else {
128        if (isZero(num) || (wt < verysmall)) {  //                              (*v[i])[lgrp].W = wt / num;
129          (*v[i])[lgrp].setToZero();  //                              (*v[i])[lgrp].N = num;
130        } else {  //                      }
131          (*v[i])[lgrp].W = wt / num;  //              }
132          (*v[i])[lgrp].N = num;  //
133        }  //              //the lowest part of the length division
134      }  //              for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
135    //                      num = 0.0;
136      //the lowest part of the length division  //                      wt = 0.0;
137      for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {  //                      for (grow = 0; grow <= lgrp - minCol; grow++) {
138        num = 0.0;  //                              tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
139        wt = 0.0;  //                              num += tmp;
140        for (grow = 0; grow <= lgrp - minCol; grow++) {  //                              wt += tmp
141          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;  //                                              * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
142          num += tmp;  //                      }
143          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);  //
144        }  //                      if (isZero(num) || (wt < verysmall)) {
145    //                              (*v[i])[lgrp].setToZero();
146        if (isZero(num) || (wt < verysmall)) {  //                      } else {
147          (*v[i])[lgrp].setToZero();  //                              (*v[i])[lgrp].W = wt / num;
148        } else {  //                              (*v[i])[lgrp].N = num;
149          (*v[i])[lgrp].W = wt / num;  //                      }
150          (*v[i])[lgrp].N = num;  //              }
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();            int maxCol = v[i]->maxCol();
# Line 222  Line 195 
195        Mat->storeMatureStock(area, age, lgrp, matnum, wt / num);        Mat->storeMatureStock(area, age, lgrp, matnum, wt / num);
196      }      }
197    
   
198      for (lgrp = maxCol - 2; lgrp >= 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 260  Line 234 
234        wt = 0.0;        wt = 0.0;
235        matnum = 0.0;        matnum = 0.0;
236        for (grow = 0; grow <= lgrp - 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 290  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    
# Line 342  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;
# Line 388  Line 365 
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 418  Line 396 
396        num = 0.0;        num = 0.0;
397        matnum = 0.0;        matnum = 0.0;
398        for (grow = 0; grow <= lgrp - 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.3  
changed lines
  Added in v.4
root@forge.cesga.es
 ViewVC Help
Powered by ViewVC 1.0.0  

Powered By FusionForge