From e55eb63276a22d445940230fd975c732f2a27906 Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Fri, 12 Apr 2024 16:15:21 +1100 Subject: [PATCH 01/10] Initial commit --- src/Map.cpp | 43 ++++++++++++++++++++++++-- src/Map.h | 11 ++++--- src/Param.h | 3 ++ src/SeapodymCoupled_EditRunCoupled.cpp | 7 +++++ src/SeapodymCoupled_OnRunCoupled.cpp | 16 +++++++--- src/SeapodymCoupled_ReadTags.cpp | 17 +++++++--- src/SeapodymDocConsole.h | 1 + src/VarParamCoupled.cpp | 15 ++++++++- 8 files changed, 96 insertions(+), 17 deletions(-) diff --git a/src/Map.cpp b/src/Map.cpp index 948dc0b..19ffbdf 100755 --- a/src/Map.cpp +++ b/src/Map.cpp @@ -103,6 +103,45 @@ void PMap::lit_map(CParam ¶m) definit_lim_infsup(nti,ntj); // appel de la fonction } +void PMap::lit_tagmap(CParam ¶m, int tagpop) +{ + const int nti=param.get_nbi(); + const int ntj=param.get_nbj(); + nbl=param.nb_layer; + +// creation des tableaux + + bord_cell.allocate(0, nti - 1, 0, ntj - 1); + nbl_bord_cell.allocate(0, nti - 1, 0, ntj - 1); + carte.allocate(0, nti - 1, 0, ntj - 1);// carte + carte.initialize(); + +// lecture du fichier + + ifstream litcarte(param.file_tag_masks[tagpop].c_str(), ios::in); + + if (!litcarte) {cout<> carte[i][j]; + if (carte[i][j]>nbl) carte[i][j] = nbl; + + } + } + litcarte.close(); + + //initialize + global = 0;//Attn: for optimization no global until fixed! + deltaX = param.deltaX; + double reso_x = param.deltaX/60.0; + domain_type((nti-2)*reso_x); +//global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed. + definit_cell_bords(nti,ntj); // appel de la fonction + definit_lim_infsup(nti,ntj); // appel de la fonction +} + void PMap::reg_indices(CParam& param) { //Regional indices @@ -138,7 +177,7 @@ char PMap::get_bord_layer_y(const int i, const int j){ return pos; } */ -void PMap::definit_cell_bords(const int nti, const int ntj) // définition des bords des cellules de la zone +void PMap::definit_cell_bords(const int nti, const int ntj) // définition des bords des cellules de la zone { CBord bordures; //bordures, objet de la classe Cbords @@ -160,7 +199,7 @@ char PMap::get_bord_layer_y(const int i, const int j){ } //*********************** // Cadre de la grille - // Definition des cotés des cellules situées aux bords de la grille + // Definition des cotés des cellules situées aux bords de la grille //*********************** // bord horizontal haut de la grille, donc ferme a gauche i = 1; diff --git a/src/Map.h b/src/Map.h index 4de7218..3d9bab1 100755 --- a/src/Map.h +++ b/src/Map.h @@ -25,6 +25,7 @@ class PMap public: void lit_map(CParam ¶m); + void lit_tagmap(CParam ¶m, int tagpop); void delete_map(const CParam ¶m) {/*DoesNothing*/} void reg_indices(CParam ¶m); @@ -35,8 +36,8 @@ class PMap IMATRIX bord_cell; // bord des cellules (voir Classe CBord plus loin) IMATRIX nbl_bord_cell; // border for the cells with all layers being exist - IMATRIX carte; // carte (mettre des zéros quand terre ou ile) - DMATRIX itopo; // carte (mettre des zéros quand terre ou ile) + IMATRIX carte; // carte (mettre des z�ros quand terre ou ile) + DMATRIX itopo; // carte (mettre des z�ros quand terre ou ile) IMATRIX maskEEZ; // carte des EEZ avec un code par EEZ IMATRIX maskMPA; // carte des MPA avec un code par MPA @@ -62,9 +63,9 @@ class PMap ivector regjmax; private: - void definit_cell_bords(const int nti, const int ntj); // définition des bords des cellules - void definit_lim_infsup(const int nti, const int ntj); // définition des limites inférieures - // et supérieures de la matrice zone pour optimiser le calcul + void definit_cell_bords(const int nti, const int ntj); // d�finition des bords des cellules + void definit_lim_infsup(const int nti, const int ntj); // d�finition des limites inf�rieures + // et sup�rieures de la matrice zone pour optimiser le calcul void domain_type(const int nlon); }; diff --git a/src/Param.h b/src/Param.h index 68f08ad..cb75f0a 100755 --- a/src/Param.h +++ b/src/Param.h @@ -140,6 +140,7 @@ class CParam string str_dir_init; string str_dir_fisheries; string str_dir_tags; + string str_dir_tagmasks; string strfile_pp; string strfile_sst; string strfile_vld; @@ -157,6 +158,7 @@ class CParam vector strfile_umc,strfile_vmc,strfile_tmc, strfile_oxymc; // array of string [nb_layer] string strdir_output; string strout_tags; + int use_tag_masks; // whether to use individual masks for tag cohorts int write_all_cohorts_dym; int write_all_fisheries_dym; @@ -279,6 +281,7 @@ class CParam vector file_catch_data; vector file_frq_data; vector file_tag_data; + vector file_tag_masks; int nb_catch_files, nb_frq_files, nb_tag_files; int tag_gauss_kernel_on; float dx_tags, dy_tags; // setup of the grid to aggregate tagging data diff --git a/src/SeapodymCoupled_EditRunCoupled.cpp b/src/SeapodymCoupled_EditRunCoupled.cpp index 719de17..c7f7de0 100755 --- a/src/SeapodymCoupled_EditRunCoupled.cpp +++ b/src/SeapodymCoupled_EditRunCoupled.cpp @@ -39,6 +39,13 @@ int SeapodymCoupled::EditRunCoupled(const char* parfile) .5*param->sp_unit_cohort[sp][a-1]+.5*param->sp_unit_cohort[sp][a]; } map.lit_map(*param); + if (param->tag_like[0] && param->use_tag_masks){ + for (int tagpop=0; tagpopstrfile_pp, param->nlong, param->nlat, param->nlevel); diff --git a/src/SeapodymCoupled_OnRunCoupled.cpp b/src/SeapodymCoupled_OnRunCoupled.cpp index 530ce82..a7fa7c3 100755 --- a/src/SeapodymCoupled_OnRunCoupled.cpp +++ b/src/SeapodymCoupled_OnRunCoupled.cpp @@ -432,10 +432,18 @@ double SeapodymCoupled::OnRunCoupled(dvar_vector x, const bool writeoutputfiles) if (!tagpop_age_solve(spop-1,t_count-1,age)) continue; } - pop.Precalrec_Calrec_adult(map,mat,*param,rw, - mat.dvarDensity[spop][age],Mortality, - tcur,fishing,age,sp,year,month,jday, - step_fishery_count,mortality_off);//checked 20150210 + + PMap* map_ptr=nullptr; + if (spop>0 && param->use_tag_masks){ + if ((spop-1) < static_cast(tagmaps.size())) + map_ptr = &tagmaps[spop-1]; + }else{ + map_ptr = ↦ + } + pop.Precalrec_Calrec_adult(*map_ptr,mat,*param,rw, + mat.dvarDensity[spop][age],Mortality, + tcur,fishing,age,sp,year,month,jday, + step_fishery_count,mortality_off);//checked 20150210 } if (sum(param->mask_fishery_sp) && !tags_only) fishing = true; } diff --git a/src/SeapodymCoupled_ReadTags.cpp b/src/SeapodymCoupled_ReadTags.cpp index 4ca3960..a96f88a 100755 --- a/src/SeapodymCoupled_ReadTags.cpp +++ b/src/SeapodymCoupled_ReadTags.cpp @@ -176,6 +176,13 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) if (rtxt){ int n = 0; for (; nuse_tag_masks){ + if (n < static_cast(tagmaps.size())) + map_ptr = &tagmaps[n]; + }else{ + map_ptr = ↦ + } rtxt >> id >> tag_no >> yy >> mm >> dd >> lat >> lon >> yy_rec >> mm_rec >> dd_rec >> lat_rec >> lon_rec >> len_rel >> len_rec; int jday_rel, jday_rec; //with default date_mode=3 use 360-days calendar: @@ -208,7 +215,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) j_mod = param->lattoj(lat_rec); if (lon_rec>xlon[0] && lon_rec<=xlon[nx_obs] && lat_rec<=ylat[0] && lat_rec>ylat[ny_obs]){ - if (map.carte(i_mod,j_mod)) + if (map_ptr->carte(i_mod,j_mod)) rec_added = 1; } @@ -217,7 +224,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) i_mod = param->lontoi(lon); j_mod = param->lattoj(lat); if (i_mod>=map.imin && i_mod<=map.imax && j_mod>=map.jmin && j_mod<=map.jmax){ - if (map.carte(i_mod,j_mod)) + if (map_ptr->carte(i_mod,j_mod)) rel_added = 1; //attn: generic way is to search over 8 surrounding cells //if (!map.carte(i_mod,j_mod) && map.carte(i_mod,j_mod-1)) { @@ -244,7 +251,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) cout << "WARNING: tag recapture is on land: " << p << " " << id << " " << yy_rec << " "<< mm_rec << " " << lat_rec << " " << lon_rec << "; mask: " - << map.carte(i,j)<< endl; + << map_ptr->carte(i,j)<< endl; } else cout<< "WARNING: tag recapture out of observational space: " << p << " " << id << " " << yy_rec << " "<< mm_rec @@ -258,7 +265,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) cout<< "WARNING: tag release out of model domain: " << p << " " << id << " " << yy << " " << mm << " " << lat << " " << lon << "; mask: " - << map.carte(i_mod,j_mod) << endl; + << map_ptr->carte(i_mod,j_mod) << endl; continue; } //One more check might be necessary: discrepancy between length(age) at recapture @@ -325,7 +332,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) int jto = jfrom + yr_tags; for (int i=ifrom; icarte[i][j]){ rec_obs(p,ii,jj) += gkernel(i,j); tlib_obs(p,ii,jj) += gkernel(i,j)*dtlib; } diff --git a/src/SeapodymDocConsole.h b/src/SeapodymDocConsole.h index 58bc39e..68c11f8 100755 --- a/src/SeapodymDocConsole.h +++ b/src/SeapodymDocConsole.h @@ -26,6 +26,7 @@ class SeapodymDocConsole VarParamCoupled* param; // objet param derive de la classe CParam VarMatrices mat; // objet mat derive de la classe CMatrices PMap map; // objet map derive de la classe CMap + vector tagmaps; // tableau de maps for tag populations VarSimtunaFunc func; // objet func derive de la classe CSimtunafunc CNumfunc nfunc; // objet nfunc derive de la classe CNumfunc CSaveTimeArea save; // objet save derive de la classe CSaveTimeArea diff --git a/src/VarParamCoupled.cpp b/src/VarParamCoupled.cpp index d722388..3f3684d 100755 --- a/src/VarParamCoupled.cpp +++ b/src/VarParamCoupled.cpp @@ -115,7 +115,8 @@ bool VarParamCoupled::read(const string& parfile) str_dir_tags = str_dir; if (!doc.get("/strdir_tags", "value").empty()) str_dir_tags = doc.get("/strdir_tags", "value"); - + if (!doc.get("/strdir_tagmasks", "value").empty()) + str_dir_tagmasks = doc.get("/strdir_tagmasks", "value"); strfile_pp = str_dir + doc.get("/strfile_pp", "value"); use_sst = 0; @@ -990,6 +991,11 @@ bool VarParamCoupled::read(const string& parfile) lonmin_tags = 180; lonmax_tags = 290;latmin_tags = -50;latmax_tags = 3; } nb_tag_files = 0; + if (!doc.get("/use_tag_masks", "flag").empty()){ + use_tag_masks=doc.getInteger("/use_tag_masks", "flag"); + }else{ + use_tag_masks=0; + } if (tag_like(sp)){ if (!doc.get("/tags_grid").empty()){ dx_tags = doc.getDouble(string("/tags_grid/reso"),"dx"); @@ -1016,6 +1022,13 @@ bool VarParamCoupled::read(const string& parfile) std::ostringstream ostr; ostr << "file" << nf+1; file_tag_data.push_back(str_dir_tags+doc.get("/file_tag_data/"+sp_name[sp],ostr.str())); + if (use_tag_masks){ + if (!str_dir_tagmasks.empty()){ + file_tag_masks.push_back(str_dir_tagmasks+file_tag_data[nf]); + }else{ + cout << "WARNING: no field in parameter file" << endl; + } + } } } else { cout << "WARNING: TAG DATA ARE ABSENT" << endl; From e621b28160218f332b8e499ab7f57ca66e5ebfdc Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Tue, 16 Apr 2024 10:54:55 +1100 Subject: [PATCH 02/10] Fixes --- src/SeapodymCoupled.h | 5 +++-- src/SeapodymCoupled_EditRunCoupled.cpp | 2 ++ src/SeapodymCoupled_Funcs.cpp | 12 ++++++++---- src/SeapodymCoupled_OnRunCoupled.cpp | 7 +++++-- src/SeapodymCoupled_OnRunFirstStep.cpp | 7 +++++-- src/SeapodymCoupled_ReadTags.cpp | 14 +++++++------- src/VarMatrices.h | 23 ++++++++++++++++++++++- src/VarParamCoupled.cpp | 5 +++-- src/dv_survival.cpp | 22 +++++++++++++++------- src/fd_survival.cpp | 17 +++++++++++++---- src/like.cpp | 23 +++++++++++++++-------- 11 files changed, 98 insertions(+), 39 deletions(-) diff --git a/src/SeapodymCoupled.h b/src/SeapodymCoupled.h index a006652..632a250 100755 --- a/src/SeapodymCoupled.h +++ b/src/SeapodymCoupled.h @@ -162,8 +162,9 @@ friend class tag_release; //void Age_class_survivals(dvar_matrix& N_a, const dmatrix N_a_1, const int nt_a, const int nt_a_1); ///void Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, dvar_matrix& Mort_a, dvar_matrix& Mort_a_1, const int a, const int sp, const bool adult); void Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp); - void Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1); - void AgePlus(dvar_matrix& N_a, dvar_matrix& N_a_1); + void Survival_tagpop(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp, const int tagpop); + void Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1, PMap& map); + void AgePlus(dvar_matrix& N_a, dvar_matrix& N_a_1, PMap& map); void Spawning(dvar_matrix& J, dvar_matrix& Hs, dvar_matrix& N_a, const int jday, const int sp, const int pop_built, const int t_count); void spawning_spinup_comp(dmatrix& J, dmatrix& Hs, dmatrix& SST, double nb_recruitment, double Tmin); diff --git a/src/SeapodymCoupled_EditRunCoupled.cpp b/src/SeapodymCoupled_EditRunCoupled.cpp index c7f7de0..b443c33 100755 --- a/src/SeapodymCoupled_EditRunCoupled.cpp +++ b/src/SeapodymCoupled_EditRunCoupled.cpp @@ -39,6 +39,8 @@ int SeapodymCoupled::EditRunCoupled(const char* parfile) .5*param->sp_unit_cohort[sp][a-1]+.5*param->sp_unit_cohort[sp][a]; } map.lit_map(*param); + nb_tagpops = param->nb_tag_files; + if (param->tag_like[0] && param->use_tag_masks){ for (int tagpop=0; tagpoptag_like[sp]){ for (int aa=a0_adult[sp]; aanb_tag_files;pop++) - total_tags += value(mat.dvarDensity(pop,aa,i,j)); + for (int pop=1; pop<=param->nb_tag_files;pop++){ + + if (!param->use_tag_masks | (i > tagmaps[pop].imin1 && i < tagmaps[pop].imax1 && j > tagmaps[pop].jinf[i] && j < tagmaps[pop].jsup[i])){ + total_tags += value(mat.dvarDensity(pop,aa,i,j)); + + if (pop==1 && value(mat.dvarDensity(pop,aa,i,j))<0) cerr << "NEGATIVE BIOMASS for " << aa << " " << value(mat.dvarDensity(pop,aa,i,j)) << endl; + } + } mat.recruit[0][i][j] += total_tags; - if (value(mat.dvarDensity(1,aa,i,j))<0) cerr << "NEGATIVE BIOMASS for " << aa << " " - << value(mat.dvarDensity(1,aa,i,j)) << endl; } } //mat.sum_B_recruit[sp] += value(mat.dvarDensity[sp][age_recruits][i][j])*W_mt(age_recruits)*lat_corrected_area; diff --git a/src/SeapodymCoupled_OnRunCoupled.cpp b/src/SeapodymCoupled_OnRunCoupled.cpp index a7fa7c3..14fbc50 100755 --- a/src/SeapodymCoupled_OnRunCoupled.cpp +++ b/src/SeapodymCoupled_OnRunCoupled.cpp @@ -486,7 +486,11 @@ double SeapodymCoupled::OnRunCoupled(dvar_vector x, const bool writeoutputfiles) } } //6.2. Ageing of population density - Survival(mat.dvarDensity[spop][a], mat.dvarDensity[spop][a-1] , a, sp); + if (spop>0 && param->use_tag_masks){ + Survival_tagpop(mat.dvarDensity[spop][a], mat.dvarDensity[spop][a-1] , a, sp, spop-1); + }else{ + Survival(mat.dvarDensity[spop][a], mat.dvarDensity[spop][a-1] , a, sp); + } } } nt_dtau=0; @@ -559,7 +563,6 @@ double SeapodymCoupled::OnRunCoupled(dvar_vector x, const bool writeoutputfiles) } // end of simulation loop - if (writeoutputfiles) {SaveDistributions(year, month); cout << "total catch in simulation (optimization): " << SUM_CATCH << endl; } diff --git a/src/SeapodymCoupled_OnRunFirstStep.cpp b/src/SeapodymCoupled_OnRunFirstStep.cpp index d9cacba..ac5a19b 100755 --- a/src/SeapodymCoupled_OnRunFirstStep.cpp +++ b/src/SeapodymCoupled_OnRunFirstStep.cpp @@ -24,7 +24,11 @@ void SeapodymCoupled::OnRunFirstStep() mat.createMatForage(map, nb_forage, t0, nbt, nbi, nbj); int nb_pops = nb_species*(1+param->nb_tag_files); - mat.CreateMatSpecies(map,t0, nbt, nbi, nbj, nb_pops, a0_adult, param->sp_nb_cohorts); + if (param->use_tag_masks){ + mat.CreateMatSpecies(map, tagmaps, t0, nbt, nbi, nbj, nb_pops, a0_adult, param->sp_nb_cohorts); + }else{ + mat.CreateMatSpecies(map,t0, nbt, nbi, nbj, nb_pops, a0_adult, param->sp_nb_cohorts); + } mat.CreateMatHabitat(map,nb_species,nb_forage,nb_layer,max(param->sp_nb_cohorts),t0, nbt,nbi,nbj,a0_adult,param->sp_nb_cohorts,param->age_compute_habitat); past_month=0; past_qtr=0; @@ -78,7 +82,6 @@ void SeapodymCoupled::OnRunFirstStep() } func.allocate_dvmatr(map.imin,map.imax,map.jinf,map.jsup); //TAG data reading and allocation section - nb_tagpops = param->nb_tag_files; if (param->tag_like[0]){ nb_rel.allocate(0,nb_tagpops-1); diff --git a/src/SeapodymCoupled_ReadTags.cpp b/src/SeapodymCoupled_ReadTags.cpp index a96f88a..f99880e 100755 --- a/src/SeapodymCoupled_ReadTags.cpp +++ b/src/SeapodymCoupled_ReadTags.cpp @@ -141,6 +141,13 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) int nbtot_tags_files = 0; for (int p=0; puse_tag_masks){ + if (p < static_cast(tagmaps.size())) + map_ptr = &tagmaps[p]; + }else{ + map_ptr = ↦ + } string file_in = param->file_tag_data[p]; //date of all recaptures in the cohort will be read from the name of the file @@ -176,13 +183,6 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec) if (rtxt){ int n = 0; for (; nuse_tag_masks){ - if (n < static_cast(tagmaps.size())) - map_ptr = &tagmaps[n]; - }else{ - map_ptr = ↦ - } rtxt >> id >> tag_no >> yy >> mm >> dd >> lat >> lon >> yy_rec >> mm_rec >> dd_rec >> lat_rec >> lon_rec >> len_rel >> len_rec; int jday_rel, jday_rec; //with default date_mode=3 use 360-days calendar: diff --git a/src/VarMatrices.h b/src/VarMatrices.h index d69b29e..eb77c7a 100755 --- a/src/VarMatrices.h +++ b/src/VarMatrices.h @@ -75,9 +75,9 @@ class VarMatrices : public CMatrices for (int age = 0; age < agemax_sp; age++) { dvarDensity(sp, age).allocate(map.imin1, map.imax1, map.jinf1, map.jsup1); - dvarDensity(sp, age).initialize(); } } + /* const int age_max_1 = 48-1; dvarDensity_age.allocate(0,age_max-1); @@ -89,6 +89,27 @@ class VarMatrices : public CMatrices } + void CreateMatSpecies(PMap& map, vector tagmaps, int t0, int nbt, int nbi, int nbj, int nb_species, const ivector a0_adult, const ivector& sp_nb_cohorts) { + //CMatrices::createMatSpecies(map, t0, nbt, nbi, nbj, nb_species, a0_adult, sp_nb_cohorts); + CMatrices::createMatSpecies(map, t0, nbt, nbi, nbj, 1, a0_adult, sp_nb_cohorts); + dvarDensity.allocate(0, nb_species - 1); + for (int sp = 0; sp < nb_species; sp++) { + //const int agemax_sp = sp_nb_cohorts(sp); + const int agemax_sp = sp_nb_cohorts(0); + dvarDensity(sp).allocate(0, agemax_sp - 1); + dvarDensity(sp).initialize(); + + for (int age = 0; age < agemax_sp; age++) { + if (sp==0){ + dvarDensity(sp, age).allocate(map.imin1, map.imax1, map.jinf1, map.jsup1); + }else{ + dvarDensity(sp, age).allocate(tagmaps[sp-1].imin1, tagmaps[sp-1].imax1, tagmaps[sp-1].jinf1, tagmaps[sp-1].jsup1); + } + dvarDensity(sp, age).initialize(); + } + } + } + void CreateMatCatch(PMap& map, int nbi, int nbj, int nb_species, const IVECTOR& nb_fleet, const ivector a0_adult, const IVECTOR& nb_cohorts, const IVECTOR& nb_region) { CMatrices::createMatCatch(map, nbi, nbj, nb_species, nb_fleet, a0_adult, nb_cohorts, nb_region); diff --git a/src/VarParamCoupled.cpp b/src/VarParamCoupled.cpp index 3f3684d..0287dbd 100755 --- a/src/VarParamCoupled.cpp +++ b/src/VarParamCoupled.cpp @@ -1021,10 +1021,11 @@ bool VarParamCoupled::read(const string& parfile) for (int nf=0; nf field in parameter file" << endl; } diff --git a/src/dv_survival.cpp b/src/dv_survival.cpp index b9e672f..5c1c69e 100755 --- a/src/dv_survival.cpp +++ b/src/dv_survival.cpp @@ -16,19 +16,27 @@ unsigned long int restore_long_int_value(void); void SeapodymCoupled::Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp) { int nb_cohorts = param->sp_nb_cohorts[sp]; + PMap* map_ptr=↦ if (asp_nb_cohorts[sp]; + PMap* map_ptr = &tagmaps[tagpop]; + if (asp_nb_cohorts[sp]; + PMap* map_ptr=↦ if (asp_nb_cohorts[sp]; + PMap* map_ptr = &tagmaps[tagpop]; + if (ause_tag_masks){ + if (p < static_cast(tagmaps.size())) + map_ptr = &tagmaps[p]; + }else{ + map_ptr = ↦ + } + const int imin = map_ptr->imin; + const int imax = map_ptr->imax; for (int i = imin; i <= imax; i++){ - const int jmin = map.jinf[i]; - const int jmax = map.jsup[i]; + const int jmin = map_ptr->jinf[i]; + const int jmax = map_ptr->jsup[i]; for (int j = jmin ; j <= jmax; j++){ - if (map.carte[i][j]){ + if (map_ptr->carte[i][j]){ double xx = param->itolon(i); double yy = param->jtolat(j); @@ -166,10 +173,10 @@ double SeapodymCoupled::get_tag_like(dvariable& likelihood, bool writeoutputs) } //temporally write recaptures to catch matrix (comment in CalcSums) for (int i = imin; i <= imax; i++){ - const int jmin = map.jinf[i]; - const int jmax = map.jsup[i]; + const int jmin = map_ptr->jinf[i]; + const int jmax = map_ptr->jsup[i]; for (int j = jmin ; j <= jmax; j++){ - if (map.carte[i][j]){ + if (map_ptr->carte[i][j]){ double xx = param->itolon(i); double yy = param->jtolat(j); if (writeoutputs){ From 2e4d9bdd351dc273dcb4ff7f615e3652dd73ab78 Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Tue, 16 Apr 2024 11:01:57 +1100 Subject: [PATCH 03/10] Fix indentation Map.cpp --- src/Map.cpp | 1301 +++++++++++++++++++++++++++------------------------ 1 file changed, 692 insertions(+), 609 deletions(-) diff --git a/src/Map.cpp b/src/Map.cpp index 19ffbdf..81344c8 100755 --- a/src/Map.cpp +++ b/src/Map.cpp @@ -1,165 +1,187 @@ #include "Map.h" //////////////////////////////////////////////////////////// -// fonction lit_map de la classe PMap +// fonction lit_map de la classe PMap // lit le fichier source et charge les donnees dans les tableaux int nbl; double deltaX; -//CBord bord_layer; - +// CBord bord_layer; void PMap::lit_map(CParam ¶m) { - const int nti=param.get_nbi(); - const int ntj=param.get_nbj(); - nbl=param.nb_layer; - -// creation des tableaux + const int nti = param.get_nbi(); + const int ntj = param.get_nbj(); + nbl = param.nb_layer; + // creation des tableaux bord_cell.allocate(0, nti - 1, 0, ntj - 1); nbl_bord_cell.allocate(0, nti - 1, 0, ntj - 1); - carte.allocate(0, nti - 1, 0, ntj - 1);// carte - itopo.allocate(0, nti - 1, 0, ntj - 1);// carte + carte.allocate(0, nti - 1, 0, ntj - 1); // carte + itopo.allocate(0, nti - 1, 0, ntj - 1); // carte carte.initialize(); - itopo.initialize(); itopo = 1; //default value + itopo.initialize(); + itopo = 1; // default value - //indices for regions + // indices for regions reg_indices(param); - -// lecture du fichier - + + // lecture du fichier + ifstream litcarte(param.str_file_mask.c_str(), ios::in); - if (!litcarte) {cout<> carte[i][j]; - if (carte[i][j]>nbl) carte[i][j] = nbl; - + if (carte[i][j] > nbl) + carte[i][j] = nbl; } } litcarte.close(); ifstream littopo(param.str_file_topo.c_str(), ios::in); - if (!littopo) {cout<> itopo[i][j]; - //necessary to avoid Ha=0 - if (carte[i][j] && itopo(i,j)==0) itopo(i,j)=0.0001; + // necessary to avoid Ha=0 + if (carte[i][j] && itopo(i, j) == 0) + itopo(i, j) = 0.0001; } } littopo.close(); - if(param.nb_EEZ) + if (param.nb_EEZ) { // mask contenant les codes des Zones exclusives - maskEEZ.allocate(0, nti-1, 0, ntj-1); + maskEEZ.allocate(0, nti - 1, 0, ntj - 1); maskEEZ.initialize(); ifstream rtxt(param.str_file_maskEEZ.c_str(), ios::in); - if (!rtxt) {cout<> maskEEZ[i][j]; } - rtxt.close(); + rtxt.close(); } - if(param.mpa_simulation) + if (param.mpa_simulation) { // mask contenant les codes des Zones exclusives - maskMPA.allocate(0, nti-1, 0, ntj-1); + maskMPA.allocate(0, nti - 1, 0, ntj - 1); maskMPA.initialize(); ifstream rtxt(param.str_file_maskMPA.c_str(), ios::in); - if (!rtxt) {cout<> maskMPA[i][j]; } - rtxt.close(); + rtxt.close(); } - - //initialize - global = 0;//Attn: for optimization no global until fixed! + // initialize + global = 0; // Attn: for optimization no global until fixed! deltaX = param.deltaX; - double reso_x = param.deltaX/60.0; - domain_type((nti-2)*reso_x); -//global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed. - definit_cell_bords(nti,ntj); // appel de la fonction - definit_lim_infsup(nti,ntj); // appel de la fonction + double reso_x = param.deltaX / 60.0; + domain_type((nti - 2) * reso_x); + // global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed. + definit_cell_bords(nti, ntj); // appel de la fonction + definit_lim_infsup(nti, ntj); // appel de la fonction } void PMap::lit_tagmap(CParam ¶m, int tagpop) { - const int nti=param.get_nbi(); - const int ntj=param.get_nbj(); - nbl=param.nb_layer; + const int nti = param.get_nbi(); + const int ntj = param.get_nbj(); + nbl = param.nb_layer; -// creation des tableaux + // creation des tableaux bord_cell.allocate(0, nti - 1, 0, ntj - 1); nbl_bord_cell.allocate(0, nti - 1, 0, ntj - 1); - carte.allocate(0, nti - 1, 0, ntj - 1);// carte + carte.allocate(0, nti - 1, 0, ntj - 1); // carte carte.initialize(); - -// lecture du fichier - + + // lecture du fichier + ifstream litcarte(param.file_tag_masks[tagpop].c_str(), ios::in); - if (!litcarte) {cout<> carte[i][j]; - if (carte[i][j]>nbl) carte[i][j] = nbl; - + if (carte[i][j] > nbl) + carte[i][j] = nbl; } } litcarte.close(); - //initialize - global = 0;//Attn: for optimization no global until fixed! + // initialize + global = 0; // Attn: for optimization no global until fixed! deltaX = param.deltaX; - double reso_x = param.deltaX/60.0; - domain_type((nti-2)*reso_x); -//global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed. - definit_cell_bords(nti,ntj); // appel de la fonction - definit_lim_infsup(nti,ntj); // appel de la fonction + double reso_x = param.deltaX / 60.0; + domain_type((nti - 2) * reso_x); + // global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed. + definit_cell_bords(nti, ntj); // appel de la fonction + definit_lim_infsup(nti, ntj); // appel de la fonction } -void PMap::reg_indices(CParam& param) +void PMap::reg_indices(CParam ¶m) { - //Regional indices + // Regional indices const int nb_regs = param.nb_region; - regimin.allocate(0,nb_regs-1); - regjmin.allocate(0,nb_regs-1); - regimax.allocate(0,nb_regs-1); - regjmax.allocate(0,nb_regs-1); - for (int a =0; a < param.nb_region; a++){ - - double lonmin = param.area[a]->lgmin; double lonmax = param.area[a]->lgmax; - double latmin = param.area[a]->ltmax; double latmax = param.area[a]->ltmin; - regimin[a] = param.lontoi(lonmin); regimax[a] = param.lontoi(lonmax); - regjmin[a] = param.lattoj(latmin); regjmax[a] = param.lattoj(latmax); - + regimin.allocate(0, nb_regs - 1); + regjmin.allocate(0, nb_regs - 1); + regimax.allocate(0, nb_regs - 1); + regjmax.allocate(0, nb_regs - 1); + for (int a = 0; a < param.nb_region; a++) + { -//cout << a+1 << "\t" << param.area[a]->lgmin << "\t" << param.area[a]->lgmax << "\t" << param.area[a]->ltmin << "\t" << param.area[a]->ltmax << endl; -//cout << a+1 << "\t" << regimin[a] << "\t" << regimax[a] << "\t" << regjmin[a] << "\t" << regjmax[a] << endl; + double lonmin = param.area[a]->lgmin; + double lonmax = param.area[a]->lgmax; + double latmin = param.area[a]->ltmax; + double latmax = param.area[a]->ltmin; + regimin[a] = param.lontoi(lonmin); + regimax[a] = param.lontoi(lonmax); + regjmin[a] = param.lattoj(latmin); + regjmax[a] = param.lattoj(latmax); + + // cout << a+1 << "\t" << param.area[a]->lgmin << "\t" << param.area[a]->lgmax << "\t" << param.area[a]->ltmin << "\t" << param.area[a]->ltmax << endl; + // cout << a+1 << "\t" << regimin[a] << "\t" << regimax[a] << "\t" << regjmin[a] << "\t" << regjmax[a] << endl; } } /* @@ -167,587 +189,648 @@ char PMap::get_bord_layer_x(const int i, const int j){ bord_layer.b = nbl_bord_cell[i][j]; char pos = bord_layer.cotex(); - return pos; + return pos; } char PMap::get_bord_layer_y(const int i, const int j){ bord_layer.b = nbl_bord_cell[i][j]; char pos = bord_layer.cotey(); - return pos; + return pos; } */ void PMap::definit_cell_bords(const int nti, const int ntj) // définition des bords des cellules de la zone { - CBord bordures; //bordures, objet de la classe Cbords - - int i, j; - - // Initialisation des cellules avec tous les cotes ouverts - // SANS = sans frontiere, tous cotes ouverts - // G_FERME = ferme a gauche - // D_FERME = ferme a droite - - bordures.b = SANS; - for (i = 1; i < nti-1; i++) - { - for (j = 1; j < ntj-1; j++) - { - bord_cell[i][j] = bordures.b; - nbl_bord_cell[i][j] = bordures.b; - } - } - //*********************** - // Cadre de la grille - // Definition des cotés des cellules situées aux bords de la grille - //*********************** - // bord horizontal haut de la grille, donc ferme a gauche - i = 1; - bordures.cote.y = SANS; - bordures.cote.x = G_FERME; -//cout << "i=1 " << (int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; - for (j = 1; j < ntj-1; j++) - { - bord_cell[i][j] = bordures.b; - nbl_bord_cell[i][j] = bordures.b; - } - - // bord horizontal bas de la grille donc ferme a droite - i = nti-2; - bordures.cote.y = SANS; - bordures.cote.x = D_FERME; -//cout << "i=nbi-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " " << bordures.b << endl; - for (j = 1; j < ntj-1; j++) - { - bord_cell[i][j] = bordures.b; - nbl_bord_cell[i][j] = bordures.b; - } - - // bord vertical gauche de la grille donc ferme a gauche - j = 1; - bordures.cote.x = SANS; - bordures.cote.y = G_FERME; -//cout << "j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; - for (i = 1; i < nti-1; i++) - { - bord_cell[i][j] = bordures.b; - nbl_bord_cell[i][j] = bordures.b; - } - - // bord vertical droit de la grille donc ferme a droite - j = ntj-2; - bordures.cote.x = SANS; - bordures.cote.y = D_FERME; -//cout << "j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; - for (i = 1; i < nti-1; i++) - { - bord_cell[i][j] = bordures.b; - nbl_bord_cell[i][j] = bordures.b; - } - - // finir les coins du cadre de la grille - - //coin en haut a gauche - bordures.cote.x = G_FERME; - bordures.cote.y = G_FERME; - bord_cell[1][1] = bordures.b; - nbl_bord_cell[1][1] = bordures.b; -//cout << "i=1; j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; - - //coin en haut a droite - bordures.cote.x = G_FERME; - bordures.cote.y = D_FERME; - bord_cell[1][ntj-2]=bordures.b; - nbl_bord_cell[1][ntj-2]=bordures.b; -//cout << "i=1; j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; - - // coin en bas a gauche - bordures.cote.x = D_FERME; - bordures.cote.y = G_FERME; - bord_cell[nti-2][1] = bordures.b; - nbl_bord_cell[nti-2][1] = bordures.b; -//cout << "i=nbi-2; j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; - - // coin en bas a droite - bordures.cote.x = D_FERME; - bordures.cote.y = D_FERME; - bord_cell[nti-2][ntj-2] = bordures.b; - nbl_bord_cell[nti-2][ntj-2] = bordures.b; -//cout << "i=nbi-2; j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; -//exit(1); - //************************ - // Interieur de la grille - // cas ou il y a des terres isolees - //************************ - for (i = 1; i < nti-1; i++) { - for (j = 1; j < ntj-1; j++){ - if (!carte[i][j]) {// si la valeur de la carte = 0 (donc terre) - if (i>1){ // si cellule n'est pas sur le bord horizontal haut - - bordures.b = bord_cell[i-1][j]; - if (bordures.cotex() == SANS) - { - bordures.cote.x = D_FERME; - bord_cell[i-1][j] = bordures.b; - } + CBord bordures; // bordures, objet de la classe Cbords + + int i, j; + + // Initialisation des cellules avec tous les cotes ouverts + // SANS = sans frontiere, tous cotes ouverts + // G_FERME = ferme a gauche + // D_FERME = ferme a droite + + bordures.b = SANS; + for (i = 1; i < nti - 1; i++) + { + for (j = 1; j < ntj - 1; j++) + { + bord_cell[i][j] = bordures.b; + nbl_bord_cell[i][j] = bordures.b; } - if (i 1) + { // si cellule n'est pas sur le bord horizontal haut + + bordures.b = bord_cell[i - 1][j]; + if (bordures.cotex() == SANS) + { + bordures.cote.x = D_FERME; + bord_cell[i - 1][j] = bordures.b; + } + } + if (i < nti - 2) + { // si cellule n'est pas sur le bord horizontal bas + + bordures.b = bord_cell[i + 1][j]; + if (bordures.cotex() == SANS) + { + bordures.cote.x = G_FERME; + bord_cell[i + 1][j] = bordures.b; + } + } + if (j > 1) + { // si cellule n'est pas sur le bord vertical gauche + + bordures.b = bord_cell[i][j - 1]; + if (bordures.cotey() == SANS) + { + bordures.cote.y = D_FERME; + bord_cell[i][j - 1] = bordures.b; + } + } + if (j < ntj - 2) + { // si cellule n'est pas sur le bord vertical droit + + bordures.b = bord_cell[i][j + 1]; + if (bordures.cotey() == SANS) + { + bordures.cote.y = G_FERME; + bord_cell[i][j + 1] = bordures.b; + } + } + bordures.cote.x = TERRE; + bordures.cote.y = TERRE; + bord_cell[i][j] = bordures.b; } - } - if (j>1){ // si cellule n'est pas sur le bord vertical gauche - - bordures.b = bord_cell[i][j-1]; - if (bordures.cotey() == SANS) - { - bordures.cote.y = D_FERME; - bord_cell[i][j-1] = bordures.b; + } // fin de boucle j + } // fin de boucle i + + // set the boundaries for nbl-layer cells + for (i = 1; i < nti - 1; i++) + { + for (j = 1; j < ntj - 1; j++) + { + if (carte[i][j] < nbl) + { // if number of layers in the cell < total nbl + if (i > 1) + { // si cellule n'est pas sur le bord horizontal haut + + bordures.b = nbl_bord_cell[i - 1][j]; + if (bordures.cotex() == SANS) + { + bordures.cote.x = D_FERME; + nbl_bord_cell[i - 1][j] = bordures.b; + } + } + if (i < nti - 2) + { // si cellule n'est pas sur le bord horizontal bas + + bordures.b = nbl_bord_cell[i + 1][j]; + if (bordures.cotex() == SANS) + { + bordures.cote.x = G_FERME; + nbl_bord_cell[i + 1][j] = bordures.b; + } + } + if (j > 1) + { // si cellule n'est pas sur le bord vertical gauche + + bordures.b = nbl_bord_cell[i][j - 1]; + if (bordures.cotey() == SANS) + { + bordures.cote.y = D_FERME; + nbl_bord_cell[i][j - 1] = bordures.b; + } + } + if (j < ntj - 2) + { // si cellule n'est pas sur le bord vertical droit + + bordures.b = nbl_bord_cell[i][j + 1]; + if (bordures.cotey() == SANS) + { + bordures.cote.y = G_FERME; + nbl_bord_cell[i][j + 1] = bordures.b; + } + } + bordures.cote.x = TERRE; + bordures.cote.y = TERRE; + nbl_bord_cell[i][j] = bordures.b; } + } // fin de boucle j + } // fin de boucle i + + //****************** + // cas ou il ya des Terres dans les bords de la grille + //****************** + i = 1; + for (j = 1; j < ntj - 1; j++) + { + bordures.b = bord_cell[i][j]; + if (carte[i + 1][j] == 0) + { + bordures.cote.x = TERRE; + bord_cell[i][j] = bordures.b; } - if (j1){ // si cellule n'est pas sur le bord horizontal haut - - bordures.b = nbl_bord_cell[i-1][j]; - if (bordures.cotex() == SANS){ - bordures.cote.x = D_FERME; - nbl_bord_cell[i-1][j] = bordures.b; - } + } + + j = 1; + for (i = 1; i < nti - 1; i++) + { + bordures.b = bord_cell[i][j]; + if (carte[i][j + 1] == 0) + { + bordures.cote.y = TERRE; + bord_cell[i][j] = bordures.b; } - if (i1){ // si cellule n'est pas sur le bord vertical gauche - - bordures.b = nbl_bord_cell[i][j-1]; - if (bordures.cotey() == SANS){ - bordures.cote.y = D_FERME; - nbl_bord_cell[i][j-1] = bordures.b; - } + } + + // once more, similarly let's set boundaries for nb_layer cells + i = 1; + for (j = 1; j < ntj - 1; j++) + { + bordures.b = nbl_bord_cell[i][j]; + if (carte[i + 1][j] < nbl) + { + bordures.cote.x = TERRE; + nbl_bord_cell[i][j] = bordures.b; } - if (j= colmin; j--) - { - int i = rowmax; - while (i >= rowmin && carte[i][j] == 0) - { i--; } - isuptmp[j] = i; - } - -//cout << "isuptmp" << isuptmp << endl; - - jmin = colmin; - while (iinftmp(jmin) > isuptmp(jmin)) jmin++; - - jmax = colmax; - while (iinftmp(jmax) > isuptmp(jmax)) jmax--; - -//cout << "jmin,jmax " << jmin << ' ' << jmax << endl; - - for (int j = jmin; j <= jmax; j++) - { - int i = iinftmp[j]; - if (!(carte[i - 1][j] == 0 && carte[i][j] != 0)) - { - //cout << carte[i][j] << ' ' << carte[i + 1][j] << endl; - //cout << i << ' ' << j << endl; - //cout << __FILE__ << ':' << __LINE__ << endl; - exit(1); - } - - i = isuptmp[j]; - if (!(carte[i][j] != 0 && carte[i + 1][j] == 0)) - { - //cout << carte[i][j] << ' ' << carte[i + 1][j] << endl; - //cout << i << ' ' << j << endl; - //cout << __FILE__ << ':' << __LINE__ << endl; - exit(1); - } - - } - - iinf.allocate(jmin, jmax); - isup.allocate(jmin, jmax); - iinf.initialize(); - isup.initialize(); - - for (int j = jmin; j <= jmax; j++) - { - iinf[j] = iinftmp[j]; - isup[j] = isuptmp[j]; - } -//cout << "iinf " << iinf << endl; -//cout << "isup " << isup << endl; -//exit(1); -//===================================== - - IVECTOR jinftmp(rowmin, rowmax); - IVECTOR jsuptmp(rowmin, rowmax); - jinftmp.initialize(); - jsuptmp.initialize(); -//cout << "jinftmp " << jinftmp << endl; -//cout << "jsuptmp " << jsuptmp << endl; - - for (int i = rowmin; i <= rowmax; i++) - { - int j = colmin; - while (j <= colmax && carte[i][j] == 0) - { j++; } - jinftmp[i] = j; - } -//cout << "jinftmp " << jinftmp << endl; -//exit(1); - - for (int i = rowmax; i >= rowmin; i--) - { - int j = colmax; - while (j >= colmin && carte[i][j] == 0) - { j--; } - jsuptmp[i] = j; - } -//cout << "jsuptmp " << jsuptmp << endl; -//exit(1); - - imin = rowmin; - while (jinftmp(imin) > jsuptmp(imin)) imin++; - imax = rowmax; - while (jinftmp(imax) > jsuptmp(imax)) imax--; -//cout << imin << endl; -//cout << imax << endl; -//exit(1); - - jinf.allocate(imin, imax); - jsup.allocate(imin, imax); - for (int i = imin; i <= imax; i++) - { - jinf[i] = jinftmp[i]; - jsup[i] = jsuptmp[i]; - } -//cout << "jinf " << jinf << endl; -//cout << "jsup " << jsup << endl; - - for (int j=jmin; j <= jmax; j++) - { - int lb = iinf[j]; - int ub = isup[j]; - for (int i = lb; i <= ub; i++) - { - if (jinf[i] > j) - { - jinf[i] = j; - } - if (jsup[i] < j) - { - jsup[i] = j; - } - } - } - -//cout << "jinf " << jinf << endl; -//cout << "jsup" << jsup << endl; -//exit(1); - - for (int i = imin; i <= imax; i++) - { - int lb = jinf[i]; - int ub = jsup[i]; - for (int j = lb; j <= ub; j++) - { - if (iinf[j] > i) - { - iinf[j] = i; - } - if (isup[j] < i) - { - isup[j] = i; - - } - } - } - - // determining extended mask edges for matrices + const int rowmin = carte.rowmin(); + const int rowmax = carte.rowmax(); + const int colmin = carte.colmin(); + const int colmax = carte.colmax(); + + // cout << "rowmin " << rowmin << endl; + // cout << "rowmax " << rowmax << endl; + // cout << "colmin " << colmin << endl; + // cout << "colmax " << colmax << endl; + + IVECTOR iinftmp(colmin, colmax); + IVECTOR isuptmp(colmin, colmax); + iinftmp.initialize(); + isuptmp.initialize(); + + // cout << "iinftmp " << iinftmp << endl; + // cout << "isuptmp " << isuptmp << endl; + + for (int j = colmin; j <= colmax; j++) + { + int i = rowmin; + while (i <= rowmax && carte[i][j] == 0) + { + i++; + } + iinftmp[j] = i; + } + + // cout << "iinftmp " << iinftmp << endl; + + for (int j = colmax; j >= colmin; j--) + { + int i = rowmax; + while (i >= rowmin && carte[i][j] == 0) + { + i--; + } + isuptmp[j] = i; + } + + // cout << "isuptmp" << isuptmp << endl; + + jmin = colmin; + while (iinftmp(jmin) > isuptmp(jmin)) + jmin++; + + jmax = colmax; + while (iinftmp(jmax) > isuptmp(jmax)) + jmax--; + + // cout << "jmin,jmax " << jmin << ' ' << jmax << endl; + + for (int j = jmin; j <= jmax; j++) + { + int i = iinftmp[j]; + if (!(carte[i - 1][j] == 0 && carte[i][j] != 0)) + { + // cout << carte[i][j] << ' ' << carte[i + 1][j] << endl; + // cout << i << ' ' << j << endl; + // cout << __FILE__ << ':' << __LINE__ << endl; + exit(1); + } + + i = isuptmp[j]; + if (!(carte[i][j] != 0 && carte[i + 1][j] == 0)) + { + // cout << carte[i][j] << ' ' << carte[i + 1][j] << endl; + // cout << i << ' ' << j << endl; + // cout << __FILE__ << ':' << __LINE__ << endl; + exit(1); + } + } + + iinf.allocate(jmin, jmax); + isup.allocate(jmin, jmax); + iinf.initialize(); + isup.initialize(); + + for (int j = jmin; j <= jmax; j++) + { + iinf[j] = iinftmp[j]; + isup[j] = isuptmp[j]; + } + // cout << "iinf " << iinf << endl; + // cout << "isup " << isup << endl; + // exit(1); + //===================================== + + IVECTOR jinftmp(rowmin, rowmax); + IVECTOR jsuptmp(rowmin, rowmax); + jinftmp.initialize(); + jsuptmp.initialize(); + // cout << "jinftmp " << jinftmp << endl; + // cout << "jsuptmp " << jsuptmp << endl; + + for (int i = rowmin; i <= rowmax; i++) + { + int j = colmin; + while (j <= colmax && carte[i][j] == 0) + { + j++; + } + jinftmp[i] = j; + } + // cout << "jinftmp " << jinftmp << endl; + // exit(1); + + for (int i = rowmax; i >= rowmin; i--) + { + int j = colmax; + while (j >= colmin && carte[i][j] == 0) + { + j--; + } + jsuptmp[i] = j; + } + // cout << "jsuptmp " << jsuptmp << endl; + // exit(1); + + imin = rowmin; + while (jinftmp(imin) > jsuptmp(imin)) + imin++; + imax = rowmax; + while (jinftmp(imax) > jsuptmp(imax)) + imax--; + // cout << imin << endl; + // cout << imax << endl; + // exit(1); + + jinf.allocate(imin, imax); + jsup.allocate(imin, imax); + for (int i = imin; i <= imax; i++) + { + jinf[i] = jinftmp[i]; + jsup[i] = jsuptmp[i]; + } + // cout << "jinf " << jinf << endl; + // cout << "jsup " << jsup << endl; + + for (int j = jmin; j <= jmax; j++) + { + int lb = iinf[j]; + int ub = isup[j]; + for (int i = lb; i <= ub; i++) + { + if (jinf[i] > j) + { + jinf[i] = j; + } + if (jsup[i] < j) + { + jsup[i] = j; + } + } + } + + // cout << "jinf " << jinf << endl; + // cout << "jsup" << jsup << endl; + // exit(1); + + for (int i = imin; i <= imax; i++) + { + int lb = jinf[i]; + int ub = jsup[i]; + for (int j = lb; j <= ub; j++) + { + if (iinf[j] > i) + { + iinf[j] = i; + } + if (isup[j] < i) + { + isup[j] = i; + } + } + } + + // determining extended mask edges for matrices // used in computing finite differences on bounds - imin1 = imin-1; - int jmin1 = jmin-1; - if (imin1 < rowmin) imin1 = rowmin; - if (jmin1 < colmin) jmin1 = colmin; - imax1 = imax+1; - int jmax1 = jmax+1; - if (imax1 > rowmax) imax1 = rowmax; - if (jmax1 > colmax) jmax1 = colmax; - - jinf1.allocate(imin1,imax1); jinf1.initialize(); - jsup1.allocate(imin1,imax1); jsup1.initialize(); - - for (int j = colmin; j <= colmax; j++){ + imin1 = imin - 1; + int jmin1 = jmin - 1; + if (imin1 < rowmin) + imin1 = rowmin; + if (jmin1 < colmin) + jmin1 = colmin; + imax1 = imax + 1; + int jmax1 = jmax + 1; + if (imax1 > rowmax) + imax1 = rowmax; + if (jmax1 > colmax) + jmax1 = colmax; + + jinf1.allocate(imin1, imax1); + jinf1.initialize(); + jsup1.allocate(imin1, imax1); + jsup1.initialize(); + + for (int j = colmin; j <= colmax; j++) + { int i = rowmin; - while (i < rowmax && carte[i+1][j] == 0) i++; + while (i < rowmax && carte[i + 1][j] == 0) + i++; iinftmp[j] = i; } - for (int j = colmax; j >= colmin; j--){ + for (int j = colmax; j >= colmin; j--) + { int i = rowmax; - while (i > rowmin && carte[i-1][j] == 0) i--; + while (i > rowmin && carte[i - 1][j] == 0) + i--; isuptmp[j] = i; } - for (int i = imin1; i <= imax1; i++){ + for (int i = imin1; i <= imax1; i++) + { jinf1[i] = jinftmp[i]; jsup1[i] = jsuptmp[i]; } - for (int j = jmin1; j <= jmax1; j++) { + for (int j = jmin1; j <= jmax1; j++) + { int lb = iinftmp[j]; int ub = isuptmp[j]; - for (int i = lb; i <= ub; i++){ - if (jinf1[i] > j-1){ - jinf1[i] = j-1; - } - if (jsup1[i] < j+1){ - jsup1[i] = j+1; - } - } - } + for (int i = lb; i <= ub; i++) + { + if (jinf1[i] > j - 1) + { + jinf1[i] = j - 1; + } + if (jsup1[i] < j + 1) + { + jsup1[i] = j + 1; + } + } + } ///////////////////////////////////////////////// -//cout << rowmin << " " << rowmax << " " << colmin << " " << colmax << endl; -//cout << imin << " " << imax << " " << imin1 << " " << imax1 << endl; -//cout << "jinf " << jinf << endl; -//cout << "jsup" << jsup << endl; -//cout << "jinf1 " << jinf1 << endl; -//cout << "jsup1" << jsup1 << endl; -//exit(1); - -/* - DMATRIX dmi(imin, imax, jinf, jsup); - dmi.initialize(); - int count = 0; - for (int i = imin; i <= imax; i++) - { - const int nonzeros = jsup[i] - jinf[i] + 1; - count += jmax - nonzeros; -cout << i << ' ' << nonzeros << ' ' << count << endl; - for (int j = jinf[i]; j <= jsup[i]; j++) - { - if (iinf[j] <= i && i <= isup[j]) - { - dmi[i][j] += 1; - } - } - } -//cout << dmi << endl; -cout << imin << endl; -cout << imax << endl; -cout << jmin << endl; -cout << jmax << endl; -cout << count << endl; -exit(1); - DMATRIX dmj(jmin, jmax, iinf, isup); - dmj.initialize(); - for (int j = jmin; j <= jmax; j++) - { - for (int i = iinf[j]; i <= isup[j]; i++) - { - if (jinf[i] <= j && j <= jsup[i]) - { - dmj[j][i] += 1; - } - } - } -cout << dmj << endl; -exit(1); -*/ + // cout << rowmin << " " << rowmax << " " << colmin << " " << colmax << endl; + // cout << imin << " " << imax << " " << imin1 << " " << imax1 << endl; + // cout << "jinf " << jinf << endl; + // cout << "jsup" << jsup << endl; + // cout << "jinf1 " << jinf1 << endl; + // cout << "jsup1" << jsup1 << endl; + // exit(1); + + /* + DMATRIX dmi(imin, imax, jinf, jsup); + dmi.initialize(); + int count = 0; + for (int i = imin; i <= imax; i++) + { + const int nonzeros = jsup[i] - jinf[i] + 1; + count += jmax - nonzeros; + cout << i << ' ' << nonzeros << ' ' << count << endl; + for (int j = jinf[i]; j <= jsup[i]; j++) + { + if (iinf[j] <= i && i <= isup[j]) + { + dmi[i][j] += 1; + } + } + } + //cout << dmi << endl; + cout << imin << endl; + cout << imax << endl; + cout << jmin << endl; + cout << jmax << endl; + cout << count << endl; + exit(1); + DMATRIX dmj(jmin, jmax, iinf, isup); + dmj.initialize(); + for (int j = jmin; j <= jmax; j++) + { + for (int i = iinf[j]; i <= isup[j]; i++) + { + if (jinf[i] <= j && j <= jsup[i]) + { + dmj[j][i] += 1; + } + } + } + cout << dmj << endl; + exit(1); + */ } -void PMap::domain_type(const int nlon) { +void PMap::domain_type(const int nlon) +{ global = 0; - if (nlon==360) + if (nlon == 360) global = 1; cout << "Is global domain? "; - if (global) cout << "YES" << endl; - else cout << "NO" << endl; + if (global) + cout << "YES" << endl; + else + cout << "NO" << endl; } - From c4189d45cea499509faa600305719d1d9289320d Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Tue, 16 Apr 2024 13:40:01 +1100 Subject: [PATCH 04/10] Security against mask holes --- src/Map.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/Map.cpp b/src/Map.cpp index 81344c8..ebcd520 100755 --- a/src/Map.cpp +++ b/src/Map.cpp @@ -586,6 +586,14 @@ void PMap::definit_lim_infsup(const int nti, const int ntj) while (iinftmp(jmax) > isuptmp(jmax)) jmax--; + for (int j = jmin; j <= jmax; j++) + { + if (iinftmp[j]>isuptmp[j]){ + cerr << "Error: There is a vertical hole in the mask" << endl; + exit(2); + } + } + // cout << "jmin,jmax " << jmin << ' ' << jmax << endl; for (int j = jmin; j <= jmax; j++) @@ -665,6 +673,14 @@ void PMap::definit_lim_infsup(const int nti, const int ntj) // cout << imax << endl; // exit(1); + for (int i = imin; i <= imax; i++) + { + if (jinftmp[i]>jsuptmp[i]){ + cerr << "Error: There is a horizontal hole in the mask" << endl; + exit(2); + } + } + jinf.allocate(imin, imax); jsup.allocate(imin, imax); for (int i = imin; i <= imax; i++) From 6cfada9ae24367d5e6d515b2f29f4ec47218dda8 Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Wed, 17 Apr 2024 09:33:57 +1100 Subject: [PATCH 05/10] Add better error message --- src/calrec_adre.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/calrec_adre.cpp b/src/calrec_adre.cpp index 83a83fc..7207649 100755 --- a/src/calrec_adre.cpp +++ b/src/calrec_adre.cpp @@ -35,7 +35,8 @@ void CCalpop::calrec1(const PMap& map, dvar_matrix& uu, const dmatrix& mortality rhs[i]=-d[i][j]*value(uu(i,j-1)) + (2*iterationNumber-e[i][j])*value(uu(i,j)) - f[i][j]*value(uu(i,j+1)); - } else { + } else { + cerr << "Error: cf file " << __FILE__ << ", line " << __LINE__ << ": itr = " << itr << ", i=" << i << ",j=" << j << endl; cout << __LINE__ << endl; exit(1); //rhs[i]=(2*iterationNumber)*value(uu(i,j)); } From 9fc40c002fedf65a595898e96b7fd7f88f738534 Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Wed, 17 Apr 2024 16:42:40 +1100 Subject: [PATCH 06/10] Fix location of saved tag files (to match master c ommit) --- src/like.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/like.cpp b/src/like.cpp index 659a936..f2960a2 100755 --- a/src/like.cpp +++ b/src/like.cpp @@ -246,7 +246,7 @@ double SeapodymCoupled::get_tag_like(dvariable& likelihood, bool writeoutputs) if (month>9) ostr << month <<15; else ostr << 0 << month <<15; - string file_out = "./" + param->sp_name[0] + "_tags_pred_" + ostr.str() + ".txt"; + string file_out = param->strout_tags + param->sp_name[0] + "_tags_pred_" + ostr.str() + ".txt"; wtxt.open(file_out.c_str(), ios::out); if (wtxt){ @@ -256,7 +256,7 @@ double SeapodymCoupled::get_tag_like(dvariable& likelihood, bool writeoutputs) } wtxt.close(); - file_out = "./" + param->sp_name[0] + "_tags_obs_" + ostr.str() + ".txt"; + file_out = param->strout_tags + param->sp_name[0] + "_tags_obs_" + ostr.str() + ".txt"; //file_out = "./" + param->sp_name[0] + "_releases_obs_" + ostr.str() + ".txt"; wtxt.open(file_out.c_str(), ios::out); if (wtxt){ From 8af8411d1a904d7730d66804e66bba400ed0e5af Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Thu, 18 Apr 2024 11:45:55 +1100 Subject: [PATCH 07/10] Add template parfile for New fields: use_tag_masks strdir_tagmasks --- template_parfile.xml | 712 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 712 insertions(+) create mode 100644 template_parfile.xml diff --git a/template_parfile.xml b/template_parfile.xml new file mode 100644 index 0000000..ac2eeec --- /dev/null +++ b/template_parfile.xml @@ -0,0 +1,712 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 3 1 3 3 3 1 1 1 1 1 1 1 1 1 3 + + + +P1 P21 P22 P23 P3 S4 S5 S6 S7 S10 S11 S12 S13 S14 P15 + + +1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + 0 1 0 0 0 1 1 1 1 1 1 1 1 1 0 + + + + + + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +epi meso mmeso bathy mbathy hmbathy + + + + + + + + + + + + + + + +skj + + + + larvae juvenile adult + + + + 1 2 34 + + + + + + + 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 600 + + + + 0 1 2 3 4 5 6 7 8 9 10 11 12 12 13 13 14 14 15 15 16 16 17 17 17 + 18 18 18 19 19 19 19 20 20 20 20 21 + + + + + 3 4.5 7.52 11.65 16.91 21.83 26.43 30.72 34.73 38.49 41.99 45.27 48.33 51.19 53.86 56.36 58.70 60.88 62.92 64.83 66.61 68.27 69.83 71.28 72.64 73.91 75.10 76.21 77.25 78.22 79.12 79.97 80.76 81.50 82.19 82.83 87.96 + + + + 0.0003 0.00109 0.00571 0.023 0.08 0.18 0.32 0.53 0.78 1.09 1.44 1.84 2.27 2.73 3.21 3.72 4.23 4.76 5.30 5.83 6.36 6.89 7.40 7.91 8.41 8.89 9.36 9.81 10.25 10.66 11.07 11.45 11.82 12.17 12.51 12.83 15.56 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1 2 3 4 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From d0a008b9f340e1cca123c604abca0f8f809c88e4 Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Thu, 25 Apr 2024 10:17:08 +1100 Subject: [PATCH 08/10] Fix SeapodymCoupled::CalcSums() --- src/SeapodymCoupled_Funcs.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SeapodymCoupled_Funcs.cpp b/src/SeapodymCoupled_Funcs.cpp index 296ae92..c5f9491 100644 --- a/src/SeapodymCoupled_Funcs.cpp +++ b/src/SeapodymCoupled_Funcs.cpp @@ -256,7 +256,7 @@ void SeapodymCoupled::CalcSums() double total_tags = 0.0; for (int pop=1; pop<=param->nb_tag_files;pop++){ - if (!param->use_tag_masks | (i > tagmaps[pop].imin1 && i < tagmaps[pop].imax1 && j > tagmaps[pop].jinf[i] && j < tagmaps[pop].jsup[i])){ + if (!param->use_tag_masks | (i > tagmaps[pop-1].imin1 && i < tagmaps[pop-1].imax1 && j > tagmaps[pop-1].jinf[i] && j < tagmaps[pop-1].jsup[i])){ total_tags += value(mat.dvarDensity(pop,aa,i,j)); if (pop==1 && value(mat.dvarDensity(pop,aa,i,j))<0) cerr << "NEGATIVE BIOMASS for " << aa << " " << value(mat.dvarDensity(pop,aa,i,j)) << endl; From 828dfd1f0c671951f8f049c19fffda04161bca4c Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Fri, 21 Jun 2024 11:45:40 +1100 Subject: [PATCH 09/10] Fix access to tagmaps when use_tag_masks is false --- src/SeapodymCoupled_Funcs.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/SeapodymCoupled_Funcs.cpp b/src/SeapodymCoupled_Funcs.cpp index c5f9491..d63ac07 100644 --- a/src/SeapodymCoupled_Funcs.cpp +++ b/src/SeapodymCoupled_Funcs.cpp @@ -256,7 +256,14 @@ void SeapodymCoupled::CalcSums() double total_tags = 0.0; for (int pop=1; pop<=param->nb_tag_files;pop++){ - if (!param->use_tag_masks | (i > tagmaps[pop-1].imin1 && i < tagmaps[pop-1].imax1 && j > tagmaps[pop-1].jinf[i] && j < tagmaps[pop-1].jsup[i])){ + int add_to_total_tags = 0; + if (!param->use_tag_masks){ + add_to_total_tags = 1; + }else if ((i > tagmaps[pop-1].imin1 && i < tagmaps[pop-1].imax1 && j > tagmaps[pop-1].jinf[i] && j < tagmaps[pop-1].jsup[i])) + { + add_to_total_tags = 1; + } + if (add_to_total_tags){ total_tags += value(mat.dvarDensity(pop,aa,i,j)); if (pop==1 && value(mat.dvarDensity(pop,aa,i,j))<0) cerr << "NEGATIVE BIOMASS for " << aa << " " << value(mat.dvarDensity(pop,aa,i,j)) << endl; From 332ad54db684fc4a91b9d6ceda54ebf6031d41d7 Mon Sep 17 00:00:00 2001 From: LucasBnnn Date: Thu, 22 Aug 2024 14:56:14 +1100 Subject: [PATCH 10/10] Add fields to example-configs/skipjack/ --- example-configs/skipjack/skipjack.xml | 2 ++ example-configs/skipjack/skipjack_F0.xml | 2 ++ 2 files changed, 4 insertions(+) diff --git a/example-configs/skipjack/skipjack.xml b/example-configs/skipjack/skipjack.xml index 0561453..f887e11 100755 --- a/example-configs/skipjack/skipjack.xml +++ b/example-configs/skipjack/skipjack.xml @@ -144,6 +144,8 @@ + + diff --git a/example-configs/skipjack/skipjack_F0.xml b/example-configs/skipjack/skipjack_F0.xml index 4132140..a08e268 100755 --- a/example-configs/skipjack/skipjack_F0.xml +++ b/example-configs/skipjack/skipjack_F0.xml @@ -144,6 +144,8 @@ + +