search

Benthic-Pelagic-Size-Spectrum-Model / lme_scale_calibration_ISMIP3a

Apache License 2.0
0 stars 0 forks source link

Check translation of `rungridbyLME` into Python #10

Open lidefi87 opened 1 month ago

lidefi87 commented 1 month ago

The rungridbyLME function from the dbpm_model_functions.R script has been translated into Python using copilot.

Currently, the translated code uses numpy and uses matrices as does the original R script. The translated code needs to be tested to ensure it runs correctly. Then we will need to adapt to it to be used with xarray so they model can be run directly onto the netcdf files.

The aim is to avoid having to extract data for a region, store it as a csv file, and then process data (done in script 01_getinputs_ISIMIP3a.R) before the model is run.

lidefi87 commented 3 weeks ago

The loop defined in line 212 is not really looping through anything. I am not sure what this section is expected to do. If it is not looping through anything, then the for loop should be removed from here. 

   #---------------------------------------------------------------------------------------
    # Numerical integration
    #---------------------------------------------------------------------------------------

    # set up with the initial values from param - same for all grid cells
    for (j in 1:Ngrid){
    U[1:Ngrid,1:(ref-1),1]<-U.init[1:(ref-1)]    #(phyto+zoo)plankton size spectrum  
    U[1:Ngrid,ref:120,1]<-U.init[ref:120]           # set initial consumer size spectrum 
    V[1:Ngrid,ref.det:120,1]<-V.init[ref.det:120]  # set initial detritivore spectrum  
    W[1:Ngrid,1]<-W.init             # set initial detritus biomass density (g.m^-3) 
    }
lidefi87 commented 3 weeks ago

From line 251 when fishing mortality is being estimated, a decimal number is used to subset data

Fvec.u[,Fref.u:Nx,1] = Fmort.u*effort[,1]
Fvec.v[,Fref.v:Nx,1] = Fmort.v*effort[,1]

The variable Fref.u comes from the gridded_params variable calculated in script 03_rungridbyLME.R. The Nx variable is an integer, but because a range is being created using Fref.u, the range misses the end of the array. 

> Fref.u:Nx
 [1] 137.3278 138.3278 139.3278 140.3278 141.3278 142.3278 143.3278 144.3278 145.3278 146.3278 147.3278
[12] 148.3278 149.3278 150.3278 151.3278 152.3278 153.3278 154.3278 155.3278 156.3278 157.3278 158.3278
[23] 159.3278 160.3278 161.3278 162.3278 163.3278 164.3278 165.3278 166.3278 167.3278 168.3278 169.3278
[34] 170.3278 171.3278 172.3278 173.3278 174.3278 175.3278 176.3278 177.3278 178.3278 179.3278 180.3278

This results in results in rows 137 to 180 being changed in lines 251 and 252 shown at the top, but the Fvec.u and Fvec.v arrays have 181 rows.

Julia to confirm if this is what was intended from the code.

lidefi87 commented 3 weeks ago

Potential issues were detected in the dbpm_model_functions.R script where matrices were changed by column instead of by row.

The lines quoted here are from the last version before I made any changes to the code (ID commit: a0fbc790f39c044cdecb93030b2b198ca9153f9b)

The loop from line 213:

for (j in 1:Ngrid){
    U[1:Ngrid,1:(ref-1),1]<-U.init[1:(ref-1)]    #(phyto+zoo)plankton size spectrum  
    U[1:Ngrid,ref:120,1]<-U.init[ref:120]           # set initial consumer size spectrum 
    V[1:Ngrid,ref.det:120,1]<-V.init[ref.det:120]  # set initial detritivore spectrum  
    W[1:Ngrid,1]<-W.init             # set initial detritus biomass density (g.m^-3) 
    }

Did not apply the vectors in the same row, but rather along the first dimension. For example, U.init[1:(ref-1)] is a vector containing the following:

[1] **7.551273e+12** 6.034953e+12 4.823118e+12 3.854624e+12 3.080608e+12 2.462017e+12 1.967641e+12
 [8] 1.572537e+12 1.256771e+12 1.004412e+12 8.027263e+11 6.415395e+11 5.127192e+11 4.097660e+11
[15] 3.274858e+11 2.617275e+11 **2.091733e+11** 1.671720e+11 1.336045e+11 1.067772e+11 8.533685e+10
[22] 6.820163e+10 5.450710e+10 4.356238e+10 3.481532e+10 2.782463e+10 2.223763e+10 1.777248e+10
[29] 1.420390e+10 1.135187e+10 9.072512e+09 7.250832e+09 **5.794933e+09** 4.631366e+09 3.701435e+09

However, U[1,1:(ref-1),1] did not contain the same information, instead it contained:

[1] **7.551273e+12** **2.091733e+11** **5.794933e+09** 1.605633e+08 4.449386e+06 1.233133e+05 1.967641e+12
  [8] 5.450710e+10 1.510136e+09 4.184417e+07 1.159604e+06 3.213959e+04 5.127192e+11 1.420390e+10
 [15] 3.935425e+08 1.090514e+07 3.022228e+05 4.823118e+12 1.336045e+11 3.701435e+09 1.025593e+08
 [22] 2.842077e+06 7.876848e+04 1.256771e+12 3.481532e+10 9.645847e+08 2.672799e+07 7.407096e+05
 [29] 2.052984e+04 3.274858e+11 9.072512e+09 2.513729e+08 6.965709e+06 1.930492e+05 3.080608e+12

I highlighted the matching numbers with double stars around them. These matching numbers are separated by 16 elements, which is the size of the first dimension of U.

The code was changed to the following:

for (j in 1:Ngrid){
      #(phyto+zoo)plankton size spectrum (from 1:(ref-1))
      # set initial consumer size spectrum  (from ref:120)
      U[j, 1:120, 1] <- U_init[1:120]  
      # set initial detritivore spectrum  
      V[j, ref_det:120, 1] <- V_init[ref_det:120]  
      # set initial detritus biomass density (g.m^-3) 
      W[j, 1] <- W_init
    }

This produces the desired output of U[1,1:(ref-1),1], which matches vector U.init[1:(ref-1)] printed above.

  [1] 7.551273e+12 6.034953e+12 4.823118e+12 3.854624e+12 3.080608e+12 2.462017e+12 1.967641e+12 1.572537e+12 1.256771e+12
 [10] 1.004412e+12 8.027263e+11 6.415395e+11 5.127192e+11 4.097660e+11 3.274858e+11 2.617275e+11 2.091733e+11 1.671720e+11
 [19] 1.336045e+11 1.067772e+11 8.533685e+10 6.820163e+10 5.450710e+10 4.356238e+10 3.481532e+10 2.782463e+10 2.223763e+10
 [28] 1.777248e+10 1.420390e+10 1.135187e+10 9.072512e+09 7.250832e+09 5.794933e+09 4.631366e+09 3.701435e+09 2.958225e+09
juliablanchard commented 3 weeks ago

Thanks Denisse, Agree with this change. For the other places where this was spotted were these inside or outside of the main time (i) and space (j) loops? cheers, Julia

Professor Julia Blanchard (she/her) ARC Future Fellow Institute for Marine and Antarctic Studies & Centre for Marine Socioecology College of Science and Engineering University of Tasmania

Please note I work Tuesdays-Fridays

E: @.**@.> T: +61 (0)3 62.26.69.32 W: http://www.utas.edu.au/profiles/staff/imas/julia-blanchardhttp://www.utas.edu.au/profiles/staff/imas/julia-blanchard ; http://www.sizeecology.com/ www.sizeecology.comhttp://www.sizeecology.com/

I acknowledge the Traditional Owners of the lands and waters on which I work and live, and pay my respects to the Elders, past and present.

[Logo Description automatically generated with low confidence]

From: Lidefi87 @.> Sent: Wednesday, 21 August 2024 10:04 AM To: Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a @.> Cc: Subscribed @.***> Subject: Re: [Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a] Check translation of rungridbyLME into Python (Issue #10)

Potential issues were detected in the dbpm_model_functions.Rhttps://github.com/Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a/blob/main/dbpm_model_functions.R script where matrices were changed by column instead of by row.

The lines quoted here are from the last version before I made any changes to the code (ID commit: https://github.com/Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a/commit/a0fbc790f39c044cdecb93030b2b198ca9153f9b

a0fbc79https://github.com/Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a/commit/a0fbc790f39c044cdecb93030b2b198ca9153f9b )

The loop from line 213:

for (j in 1:Ngrid){ U[1:Ngrid,1:(ref-1),1]<-U.init[1:(ref-1)] #(phyto+zoo)plankton size spectrum U[1:Ngrid,ref:120,1]<-U.init[ref:120] # set initial consumer size spectrum V[1:Ngrid,ref.det:120,1]<-V.init[ref.det:120] # set initial detritivore spectrum W[1:Ngrid,1]<-W.init # set initial detritus biomass density (g.m^-3) }

Did not apply the vectors in the same row, but rather along the first dimension. For example, U.init[1:(ref-1)] is a vector containing the following:

[1] 7.551273e+12 6.034953e+12 4.823118e+12 3.854624e+12 3.080608e+12 2.462017e+12 1.967641e+12 [8] 1.572537e+12 1.256771e+12 1.004412e+12 8.027263e+11 6.415395e+11 5.127192e+11 4.097660e+11 [15] 3.274858e+11 2.617275e+11 2.091733e+11 1.671720e+11 1.336045e+11 1.067772e+11 8.533685e+10 [22] 6.820163e+10 5.450710e+10 4.356238e+10 3.481532e+10 2.782463e+10 2.223763e+10 1.777248e+10 [29] 1.420390e+10 1.135187e+10 9.072512e+09 7.250832e+09 5.794933e+09 4.631366e+09 3.701435e+09

However, U[1,1:(ref-1),1] did not contain the same information, instead it contained:

[1] 7.551273e+12 2.091733e+11 5.794933e+09 1.605633e+08 4.449386e+06 1.233133e+05 1.967641e+12 [8] 5.450710e+10 1.510136e+09 4.184417e+07 1.159604e+06 3.213959e+04 5.127192e+11 1.420390e+10 [15] 3.935425e+08 1.090514e+07 3.022228e+05 4.823118e+12 1.336045e+11 3.701435e+09 1.025593e+08 [22] 2.842077e+06 7.876848e+04 1.256771e+12 3.481532e+10 9.645847e+08 2.672799e+07 7.407096e+05 [29] 2.052984e+04 3.274858e+11 9.072512e+09 2.513729e+08 6.965709e+06 1.930492e+05 3.080608e+12

I bold I highlighted the matching numbers, which are separated by 16 elements, which is the size of the first dimension of U.

The code was changed to the following:

for (j in 1:Ngrid){

(phyto+zoo)plankton size spectrum (from 1:(ref-1))

  # set initial consumer size spectrum  (from ref:120)
  U[j, 1:120, 1] <- U_init[1:120]
  # set initial detritivore spectrum
  V[j, ref_det:120, 1] <- V_init[ref_det:120]
  # set initial detritus biomass density (g.m^-3)
  W[j, 1] <- W_init
}

This produces the desired output of U[1,1:(ref-1),1], which matches vector U.init[1:(ref-1)] printed above.

[1] 7.551273e+12 6.034953e+12 4.823118e+12 3.854624e+12 3.080608e+12 2.462017e+12 1.967641e+12 1.572537e+12 1.256771e+12 [10] 1.004412e+12 8.027263e+11 6.415395e+11 5.127192e+11 4.097660e+11 3.274858e+11 2.617275e+11 2.091733e+11 1.671720e+11 [19] 1.336045e+11 1.067772e+11 8.533685e+10 6.820163e+10 5.450710e+10 4.356238e+10 3.481532e+10 2.782463e+10 2.223763e+10 [28] 1.777248e+10 1.420390e+10 1.135187e+10 9.072512e+09 7.250832e+09 5.794933e+09 4.631366e+09 3.701435e+09 2.958225e+09

— Reply to this email directly, view it on GitHubhttps://github.com/Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a/issues/10#issuecomment-2299963777, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABGD7OWB5RVWOUGUKGQKZP3ZSPKR7AVCNFSM6AAAAABMPOTJBCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEOJZHE3DGNZXG4. You are receiving this because you are subscribed to this thread.

This email is confidential, and is for the intended recipient only. Access, disclosure, copying, distribution, or reliance on any of it by anyone outside the intended recipient organisation is prohibited and may be a criminal offence. Please delete if obtained in error and email confirmation to the sender. The views expressed in this email are not necessarily the views of the University of Tasmania, unless clearly intended otherwise.

lidefi87 commented 3 weeks ago

Hi @juliablanchard, yes, there were a couple other occasions where I spotted this. I will share them here so we can have a record of why we changed the code for future reference

lidefi87 commented 3 weeks ago

Lines 270-272 in the original script were:

Y.u[,Fref.u:Nx,1]<-Fvec.u[,Fref.u:Nx,1]*U[,Fref.u:Nx,1]*10^x[Fref.u:Nx] 
 #output fisheries catches per yr at size
 Y.v[,Fref.v:Nx,1]<-Fvec.v[,Fref.v:Nx,1]*V[,Fref.v:Nx,1]*10^x[Fref.v:Nx]

I am focusing on the first line, but the same applies to both. The dimensions of each elements on the right hand side of this equation were as follows:

> dim(Fvec.u[,Fref.u:Nx,1])
[1] 16 44
> dim(U[,Fref.u:Nx,1])
[1] 16 44
> dim(10^x[Fref.u:Nx])
NULL
> length(10^x[Fref.u:Nx])
[1] 44

We were expecting the vector (10^x[Fref.u:Nx]) will be multiplied along each row of the matrices. However, this did not happen. Multiplying the matrices gave us this result (showing the first few columns and rows below):

> Fvec.u[,Fref.u:Nx,1]*U[,Fref.u:Nx,1]
              [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]
 [1,] 7.217646e-08 1.135655e-08 1.654789e-09 2.221844e-10 2.733968e-11 3.065023e-12 3.111112e-13
 [2,] 7.217646e-08 1.135655e-08 1.654789e-09 2.221844e-10 2.733968e-11 3.065023e-12 3.111112e-13
 [3,] 7.217646e-08 1.135655e-08 1.654789e-09 2.221844e-10 2.733968e-11 3.065023e-12 3.111112e-13
 [4,] 7.217646e-08 1.135655e-08 1.654789e-09 2.221844e-10 2.733968e-11 3.065023e-12 3.111112e-13
 [5,] 7.217646e-08 1.135655e-08 1.654789e-09 2.221844e-10 2.733968e-11 3.065023e-12 3.111112e-13
 [6,] 1.443529e-08 2.271310e-09 3.309578e-10 4.443688e-11 5.467936e-12 6.130046e-13 6.222225e-14

The vector included in this multiplication looks as follows:

> 10^x[Fref.u:Nx]
 [1]     39.81072     50.11872     63.09573     79.43282    100.00000    125.89254    158.48932
 [8]    199.52623    251.18864    316.22777    398.10717    501.18723    630.95734    794.32823
[15]   1000.00000   1258.92541   1584.89319   1995.26231   2511.88643   3162.27766   3981.07171
[22]   5011.87234   6309.57344   7943.28235  10000.00000  12589.25412  15848.93192  19952.62315
[29]  25118.86432  31622.77660  39810.71706  50118.72336  63095.73445  79432.82347 100000.00000
[36] 125892.54118 158489.31925 199526.23150 251188.64315 316227.76602 398107.17055 501187.23363
[43] 630957.34448 794328.23472

The result we got was:

> (Fvec.u[,Fref.u:Nx,1]*U[,Fref.u:Nx,1]*10^x[Fref.u:Nx])[1:10, 1:5]
              [,1]         [,2]         [,3]         [,4]         [,5]
 [1,] 2.873397e-06 1.799892e-05 1.044101e-04 2.221844e-08 1.088412e-07
 [2,] 3.617392e-06 2.265929e-05 1.314446e-04 2.797136e-08 1.370230e-07
 [3,] 4.554027e-06 2.852636e-05 1.654789e-04 3.521385e-08 1.725017e-07
 [4,] 5.733180e-06 3.591256e-05 2.083256e-04 4.433162e-08 2.171668e-07
 [5,] 7.217646e-06 4.521124e-05 2.622664e-04 5.581020e-08 2.733968e-07
 [6,] 1.817296e-06 1.138351e-05 6.603476e-05 1.405218e-08 6.883724e-08
 [7,] 2.287840e-06 1.433100e-05 8.313283e-05 1.769064e-08 8.666095e-08
 [8,] 2.880219e-06 1.804165e-05 1.046580e-04 2.227120e-08 1.090997e-07
 [9,] 3.625982e-06 2.271310e-05 1.317567e-04 2.803778e-08 1.373483e-07
[10,] 4.564840e-06 2.859410e-05 1.658718e-04 3.529747e-08 1.729113e-07

To obtain the first number [1,1] 2.873397e-06, the first element of the vector and matrix were multiplied, that is 39.81072 x 7.217646e-08. However, to obtain [1,2] 1.799892e-05, matrix element [1,2] was multiplied by element 17 of the vector, that is 1.135655e-8 x 1584.89319. As discussed above, the multiplication was applied along the first axis of the matrix which has a length of 16.

The code was changed to

catch_u <- Fvec_u[, Fref_u:Nx, 1]*U[, Fref_u:Nx, 1]
Y_u[, Fref_u:Nx, 1] <- t(apply(catch_u, 1, function(i) i*10^x[Fref_u:Nx]))

#output fisheries catches per yr at size
catch_v <- Fvec_v[, Fref_v:Nx, 1]*V[, Fref_v:Nx, 1]
 Y_v[, Fref_v:Nx, 1] <- t(apply(catch_v, 1, function(i) i*10^x[Fref_v:Nx]))

This produced the following result:

> t(apply(catch_u, 1, function(i) i*10^x[Fref_u:Nx]))[1:10, 1:5]
              [,1]         [,2]         [,3]         [,4]         [,5]
 [1,] 2.873397e-06 5.691757e-07 1.044101e-07 1.764873e-08 2.733968e-09
 [2,] 2.873397e-06 5.691757e-07 1.044101e-07 1.764873e-08 2.733968e-09
 [3,] 2.873397e-06 5.691757e-07 1.044101e-07 1.764873e-08 2.733968e-09
 [4,] 2.873397e-06 5.691757e-07 1.044101e-07 1.764873e-08 2.733968e-09
 [5,] 2.873397e-06 5.691757e-07 1.044101e-07 1.764873e-08 2.733968e-09
 [6,] 5.746793e-07 1.138351e-07 2.088202e-08 3.529747e-09 5.467936e-10
 [7,] 5.746793e-07 1.138351e-07 2.088202e-08 3.529747e-09 5.467936e-10
 [8,] 5.746793e-07 1.138351e-07 2.088202e-08 3.529747e-09 5.467936e-10
 [9,] 5.746793e-07 1.138351e-07 2.088202e-08 3.529747e-09 5.467936e-10
[10,] 5.746793e-07 1.138351e-07 2.088202e-08 3.529747e-09 5.467936e-10

Multiplying vector (10^x[Fref_u:Nx]) along each row of the matrix.

lidefi87 commented 3 weeks ago

From line 285:

for (i in 1:(Neq)) { # time

      if (i < Neq){

        # get proportion of total fishable biomass for each grid cell
        #output rates of fisheries catches per yr at size
        prop.b<-(apply(U[,Fref.u:Nx,i]*10^x[Fref.u:Nx]*dx,1,sum) + apply(V[,Fref.v:Nx,i]*10^x[Fref.v:Nx]*dx,1,sum))/sum(apply(U[,Fref.u:Nx,i]*10^x[Fref.u:Nx]*dx,1,sum) + apply(V[,Fref.v:Nx,i]*10^x[Fref.v:Nx]*dx,1,sum))
        #check sums to 1
        #sum(prop.b)

        # redistribute total effort across gridcells according to proportion of biomass in that gridcell
        effort[,i+1] = gravitymodel(effort[,i+1], prop.b,depth = depth,iter=1)
        # for option 2 iter = 1

      } # end gravity model

The matrix and vector multiplications have the same issues reported above. We will take U[,Fref.u:Nx,i]*10^x[Fref.u:Nx]*dx as an example. The dimensions of these variables are as follows:

> dim(U[,Fref.u:Nx,i])
[1] 16 44
> dim(10^x[Fref.u:Nx]*dx)
NULL
> length(10^x[Fref.u:Nx]*dx)
[1] 44

The content of each variable are as follows:

> U[,Fref.u:Nx,i][1:10,1:7]
              [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]
 [1,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [2,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [3,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [4,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [5,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [6,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [7,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [8,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
 [9,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
[10,] 7.356206e-09 1.157456e-09 1.686556e-10 2.264497e-11 2.786453e-12 3.123863e-13 3.170837e-14
> 10^x[Fref.u:Nx]*dx
 [1]     3.981072     5.011872     6.309573     7.943282    10.000000    12.589254    15.848932
 [8]    19.952623    25.118864    31.622777    39.810717    50.118723    63.095734    79.432823
[15]   100.000000   125.892541   158.489319   199.526231   251.188643   316.227766   398.107171
[22]   501.187234   630.957344   794.328235  1000.000000  1258.925412  1584.893192  1995.262315
[29]  2511.886432  3162.277660  3981.071706  5011.872336  6309.573445  7943.282347 10000.000000
[36] 12589.254118 15848.931925 19952.623150 25118.864315 31622.776602 39810.717055 50118.723363
[43] 63095.734448 79432.823472

We expect that the vector (10^x[Fref.u:Nx]*dx) will be multiplied by each row in the matrix, but instead we get this result:

(U[,Fref.u:Nx,i]*10^x[Fref.u:Nx]*dx)[1:10, 1:7]
              [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]
 [1,] 2.928558e-08 1.834445e-07 1.064145e-06 2.264497e-10 1.109307e-09 4.950989e-09 7.964783e-13
 [2,] 3.686836e-08 2.309429e-07 1.339679e-06 2.850833e-10 1.396535e-09 6.232926e-09 1.002707e-12
 [3,] 4.641452e-08 2.907399e-07 1.686556e-06 3.588987e-10 1.758133e-09 7.846789e-09 1.262333e-12
 [4,] 5.843242e-08 3.660199e-07 2.123249e-06 4.518266e-10 2.213358e-09 9.878523e-09 1.589183e-12
 [5,] 7.356206e-08 4.607917e-07 2.673012e-06 5.688160e-10 2.786453e-09 1.243632e-08 2.000663e-12
 [6,] 9.260914e-08 5.801024e-07 3.365122e-06 7.160970e-10 3.507936e-09 1.565640e-08 2.518686e-12
 [7,] 1.165880e-07 7.303056e-07 4.236438e-06 9.015127e-10 4.416230e-09 1.971024e-08 3.170837e-12
 [8,] 1.467756e-07 9.194003e-07 5.333360e-06 1.134937e-09 5.559704e-09 2.481373e-08 3.991848e-12
 [9,] 1.847795e-07 1.157456e-06 6.714302e-06 1.428801e-09 6.999253e-09 1.243632e-12 5.025439e-12
[10,] 2.326236e-07 1.457151e-06 8.452805e-06 1.798754e-09 8.811538e-09 1.565640e-12 6.326652e-12

Once again, the vector was multiplied along the first dimension of the matrix. That is to obtain the result in [1,2] in the matrix above 1.834445e-07, element 17 in the vector was multiplied by matrix element [1,2] 158.489319 x 1.157456e-9.

This was changed to:

x_u <- 10^x[Fref_u:Nx]*dx
x_v <- 10^x[Fref_v:Nx]*dx

if (i < Neq){
        # get proportion of total fishable biomass for each grid cell
        #output rates of fisheries catches per yr at size
        prop_bu <- t(apply(U[, Fref_u:Nx, i], 1, function(j) j*x_u))
        prop_bv <- t(apply(V[, Fref_v:Nx, i], 1, function(j) j*x_v))
}

which resulted in:

> (t(apply(U[, Fref_u:Nx, i], 1, function(j) j*x_u)))[1:10,1:7]
              [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]
 [1,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [2,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [3,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [4,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [5,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [6,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [7,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [8,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
 [9,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13
[10,] 2.928558e-08 5.801024e-09 1.064145e-09 1.798754e-10 2.786453e-11 3.932711e-12 5.025439e-13

The vector (10^x[Fref.u:Nx]) was moved out of the loop as it remained the same regardless of the value of i in the loop, so it does not need to be recalculated at every time step.

To improve readability the original code (at the beginning of this comment) was changed to:

for(i in 1:(Neq)){ 
     #These variables are not recalculated at every time step because they remain the same
      x_u <- 10^x[Fref_u:Nx]*dx
      x_v <- 10^x[Fref_v:Nx]*dx

      if (i < Neq){
        # get proportion of total fishable biomass for each grid cell
        #output rates of fisheries catches per yr at size
       #each row in the matrix is multiplied by the vector and all elements are summed
        prop_bu <- rowSums(t(apply(U[, Fref_u:Nx, i], 1, function(j) j*x_u)))
        prop_bv <- rowSums(t(apply(V[, Fref_v:Nx, i], 1, function(j) j*x_v)))

       #Proportions are are calculated 
        prop_b <- (prop_bu+prop_bv)/sum(prop_bu+prop_bv)
        #Proportion should add up to 1
        #sum(prop_b)

        # redistribute total effort across gridcells according to proportion of 
        # biomass in that gridcell
        effort[, i+1] <- gravitymodel(effort[, i+1], prop_b, depth = depth,
                                      iter = 1)
        # for option 2 iter = 1

      } # end gravity model
lidefi87 commented 3 weeks ago

Line 296 uses a gravitymodel function, which is defined as follows:

gravitymodel<-function(effort=effort[,3000],prop.b, depth, iter){

  # redistribute total effort across gridcells according to proportion of biomass in that gridcell
  # using graivity model, Walters & Bonfil, 1999, Gelchu & Pauly 2007
  # ideal free distribution - Blanchard et al 2008
  # new_effort= prop.b*effort

  eff<-as.vector(effort)
  d<-unlist(as.vector(depth))

  for(j in 1:iter) {

    a<-eff
    # suit = prop.b*(1-d/max(d))*(1-a/0.001)
    # # CN option 2 from below 
    suit = prop.b*(1-d/max(d))
    # # CN option 3 from below 
    # suit = prop.b
    # rescale:
    rel.suit = suit/sum(suit)
    neweffort <- eff + rel.suit*eff
    mult <- sum(eff)/sum(neweffort)
    eff<-neweffort*mult     #gradient drives individuals to best locations at equilibrium (assumed to reached after 10000 iterations)

  }

  return(eff)
}

The variable a is never used (a<-eff) so it has been removed.

lidefi87 commented 3 weeks ago

Total biomass eaten by predators in line 396 is calculated as follows: eatenbypred <- 10^x*f_pel*U[j, ,i]+10^x*f_ben*U[j, ,i]

Need to check if this is correct, or if it should be: eatenbypred <- 10^(x*f_pel)*U[j, ,i]+10^(x*f_ben)*U[j, ,i]

I do not know the source of these equations, but before this point, x is multiplied by a constant before being used as a power.

This will not be changed, until confirmation is received.

Same question for detritus outputs in line 402: output_w <- sum(10^x*f_det*V[j, ,i]*dx)

Should this change to: output_w <- sum(10^(x*f_det)*V[j, ,i]*dx)?

juliablanchard commented 3 weeks ago

x is log10 ( body mass) so this raises it back to body mass

———— Sent from my phone

From: Lidefi87 @.> Sent: Wednesday, August 21, 2024 4:16:09 PM To: Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a @.> Cc: Julia Blanchard @.>; Mention @.> Subject: Re: [Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a] Check translation of rungridbyLME into Python (Issue #10)

Total biomass eaten by predators in line 396 is calculated as follows: eatenbypred <- 10^xf_pelU[j, ,i]+10^xf_benU[j, ,i]

Need to check if this is correct, or if it should be: eatenbypred <- 10^(xf_pel)U[j, ,i]+10^(xf_ben)U[j, ,i]

I do not know the source of these equations, but before this point, x is multiplied by a constant before being used as a power.

This will not be changed, until confirmation is received.

— Reply to this email directly, view it on GitHubhttps://github.com/Benthic-Pelagic-Size-Spectrum-Model/lme_scale_calibration_ISMIP3a/issues/10#issuecomment-2301215039, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABGD7OS36JO66O7TL2I4GOTZSQWCTAVCNFSM6AAAAABMPOTJBCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBRGIYTKMBTHE. You are receiving this because you were mentioned.Message ID: @.***>

This email is confidential, and is for the intended recipient only. Access, disclosure, copying, distribution, or reliance on any of it by anyone outside the intended recipient organisation is prohibited and may be a criminal offence. Please delete if obtained in error and email confirmation to the sender. The views expressed in this email are not necessarily the views of the University of Tasmania, unless clearly intended otherwise.

lidefi87 commented 3 weeks ago

Lines 491 to 493 have an if statement, but regardless of the condition met, the result is the same:

if (use.init==T)U[j,1:ref,i+1]<-ui0[j,i]*10^(r.plank[j,i]*x)[1:(ref)] 

if (use.init==F)U[j,1:ref,i+1]<-ui0[j,i]*10^(r.plank[j,i]*x)[1:(ref)]

This has been replaced by: U[j, 1:ref, i+1] <- ui0[j, i]*10^(r_plank[j, i]*x)[1:(ref)]

lidefi87 commented 3 weeks ago

The burnin.len parameter from the gridded_sizemodel function was removed because it was not used anywhere within the function, and it was never declared in the 03_rungridbyLME.R script