1
#include <Rcpp.h>
2
// [[Rcpp::plugins(cpp11)]]
3
// [[Rcpp::depends(RcppEigen)]]
4
// [[Rcpp::depends(BH)]]
5

6
using namespace Rcpp;
7

8

9
#include <vector>
10
#include <array>
11
#include <string>
12
#include <algorithm>
13
#include <unordered_map>
14
#include <iomanip>
15
#include <iostream>
16
#include <RcppEigen.h>
17
#include "functions.h"
18
#include <boost/functional/hash.hpp>
19

20
// [[Rcpp::export]]
21 4
Rcpp::List rcpp_generate_model_object(Rcpp::S4 opts, bool unreliable_formulation, Rcpp::S4 data, bool verbose) {
22
  //// Initialization
23
  if (verbose) Rcpp::Rcout << "Initialization" << std::endl;
24
  // create variables
25 4
  std::size_t counter=0;
26
  std::size_t n_pu_INT;
27
  std::size_t n_edges_INT;
28 4
  std::size_t n_edges2_INT=0;
29
  std::size_t n_species_INT;
30
  bool boundary;
31 4
  double boundary_threshold=1.0e-10;
32 4
  double zero_adjust=0.0;
33

34
  //// Preliminary processing
35
  if (verbose) Rcpp::Rcout << "Preliminary processing" << std::endl;
36 4
  Rcpp::checkUserInterrupt();
37
  /// extract parameters from Rcpp::S4 opts
38 4
  double blm_DBL=Rcpp::as<double>(opts.slot("BLM"));
39 4
  boundary = blm_DBL > boundary_threshold;
40 4
  double failure_multiplier_DBL=0.0;
41 4
  std::size_t maxrlevel_INT=0;
42 4
  if (!unreliable_formulation) {
43 4
    failure_multiplier_DBL=Rcpp::as<double>(opts.slot("failure.multiplier"));
44 4
    maxrlevel_INT=Rcpp::as<double>(opts.slot("max.r.level"));
45
  }
46

47
  /// extract data from Rcpp::S4 data
48
  // species data
49
  if (verbose) Rcpp::Rcout << "\tSpecies data" << std::endl;
50 4
  Rcpp::DataFrame species_DF=Rcpp::as<Rcpp::DataFrame>(data.slot("species"));
51 4
  n_species_INT=species_DF.nrows();
52
  // planning unit data
53
  if (verbose) Rcpp::Rcout << "\tpu data" << std::endl;
54 4
  Rcpp::DataFrame pu_DF=Rcpp::as<Rcpp::DataFrame>(data.slot("pu"));
55 4
  std::vector<double> pu_DF_area = pu_DF["area"];
56 4
  std::vector<double> pu_DF_cost = pu_DF["cost"];
57 4
  std::vector<std::size_t> pu_DF_status = pu_DF["status"];
58 4
  n_pu_INT=pu_DF_area.size();
59

60
  // pu.species.probabilities
61
  if (verbose) Rcpp::Rcout << "\tpuvspecies data" << std::endl;
62 4
  Rcpp::DataFrame puvspecies_DF=Rcpp::as<Rcpp::DataFrame>(data.slot("pu.species.probabilities"));
63 4
  std::vector<std::size_t> puvspecies_DF_pu = puvspecies_DF["pu"];
64 4
  std::vector<std::size_t> puvspecies_DF_species = puvspecies_DF["species"];
65 4
  std::vector<double> puvspecies_DF_value = puvspecies_DF["value"];
66 4
  for (std::size_t i=0; i<puvspecies_DF_pu.size(); ++i) {
67 4
    puvspecies_DF_pu[i]-=1;
68 4
    puvspecies_DF_species[i]-=1;
69
  }
70

71
  // boundary
72
  if (verbose) Rcpp::Rcout << "\tboundary data" << std::endl;
73 4
  Rcpp::DataFrame boundary_DF=Rcpp::as<Rcpp::DataFrame>(data.slot("boundary"));
74 4
  std::vector<std::size_t> boundary_DF_id1 = boundary_DF["id1"];
75 4
  std::vector<std::size_t> boundary_DF_id2 = boundary_DF["id2"];
76 4
  std::vector<double> boundary_DF_boundary = boundary_DF["boundary"];
77 4
  n_edges_INT=boundary_DF_boundary.size();
78 4
  std::vector<std::size_t> edge_pos_INT;
79 4
  std::vector<std::pair<std::size_t, std::size_t>> edge_id_PAIR;
80 4
  if (boundary) {
81 0
    edge_pos_INT.reserve(n_edges_INT);
82 0
    for (std::size_t i=0; i<n_edges_INT; ++i) {
83
      // convert to base-0 indexing
84 0
      boundary_DF_id1[i]-=1;
85 0
      boundary_DF_id2[i]-=1;
86
      // main processing
87 0
      if (boundary_DF_id1[i]!=boundary_DF_id2[i]) {
88
        /// if boundary_DF_id1[i] != boundary_DF_id2[i]
89
        // store quadratic variable
90 0
        edge_id_PAIR.push_back(std::pair<std::size_t,std::size_t>(boundary_DF_id1[i], boundary_DF_id2[i]));
91
        // cache location of quadratic variable
92 0
        edge_pos_INT.push_back(i);
93
      }
94
      // increase cost variable with boundary
95 0
      pu_DF_cost[boundary_DF_id1[i]] += (blm_DBL * boundary_DF_boundary[i]);
96
    }
97 0
    edge_id_PAIR.shrink_to_fit();
98 0
    n_edges2_INT=edge_id_PAIR.size();
99
  }
100

101
  /// attribute.space
102
  if (verbose) Rcpp::Rcout << "\tattribute space data" << std::endl;
103 4
  Rcpp::List attributespaces_LST=Rcpp::as<Rcpp::List>(data.slot("attribute.spaces"));
104 4
  std::size_t n_attribute_spaces_INT=attributespaces_LST.size();
105
  
106
  /// extract objects from S4 AttributeSpaces objects
107
  if (verbose) Rcpp::Rcout << "\tAttributeSpaces data" << std::endl;
108
  // important variables
109 4
  std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> species_space_dpcoords_MTX;
110 4
  std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> species_space_pucoords_MTX;
111 4
  std::vector<Eigen::VectorXd> species_space_dpweights_VXD;
112 4
  std::vector<Rcpp::IntegerVector> species_space_puids_RIV;
113 4
  std::vector<std::size_t> species_space_ndp_INT;
114 4
  std::vector<std::size_t> species_space_npu_INT;
115 4
  std::vector<std::size_t> species_attributespace_species_INT;
116 4
  std::vector<std::size_t> species_attributespace_space_INT;
117
  
118 4
  std::size_t preallocate_INT=n_species_INT*n_attribute_spaces_INT;
119 4
  species_space_dpcoords_MTX.reserve(preallocate_INT);
120 4
  species_space_pucoords_MTX.reserve(preallocate_INT);
121 4
  species_space_dpweights_VXD.reserve(preallocate_INT);
122 4
  species_space_puids_RIV.reserve(preallocate_INT);
123 4
  species_space_ndp_INT.reserve(preallocate_INT);
124 4
  species_space_npu_INT.reserve(preallocate_INT);
125 4
  species_attributespace_species_INT.reserve(preallocate_INT);
126 4
  species_attributespace_space_INT.reserve(preallocate_INT);
127
  {
128
    // temporary variables
129
    std::size_t curr_spp;
130 4
    Rcpp::S4 tmp_dp_S4;
131 4
    Rcpp::S4 tmp_pu_S4;
132 4
    Rcpp::S4 tmp1_S4;
133 4
    Rcpp::S4 tmp2_S4;
134 4
    Rcpp::List tmp_LST;
135 4
    Rcpp::NumericVector currWeights_DBL;
136 4
    Rcpp::IntegerVector currIds_INT;
137 4
    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp1_MTX;
138 4
    Rcpp::NumericMatrix tmp2_MTX;
139
    double *pcv;
140
    
141 4
    for (std::size_t i=0; i<n_attribute_spaces_INT; ++i) {
142
      if (verbose) Rcpp::Rcout << "\t\tstarting AttributeSpaces " << i << std::endl;
143
      // init
144 4
      tmp_LST=Rcpp::as<Rcpp::List>(Rcpp::as<Rcpp::S4>(attributespaces_LST[i]).slot("spaces"));
145 4
      for (std::size_t j=0; j<tmp_LST.size(); ++j) {
146
        // init
147
        if (verbose) Rcpp::Rcout << "\t\t\tstarting species " << j << std::endl;
148 4
        tmp1_S4=Rcpp::as<Rcpp::S4>(tmp_LST[j]);
149 4
        tmp_pu_S4=Rcpp::as<Rcpp::S4>(tmp1_S4.slot("planning.unit.points"));
150 4
        tmp_dp_S4=Rcpp::as<Rcpp::S4>(tmp1_S4.slot("demand.points"));
151 4
        curr_spp=Rcpp::as<std::size_t>(tmp1_S4.slot("species"));
152 4
        species_attributespace_species_INT.push_back(curr_spp-1);
153 4
        species_attributespace_space_INT.push_back(i);
154

155
        // planning unit points
156
        if (verbose) Rcpp::Rcout << "\t\t\tpu points coordinates" << std::endl;
157 4
        tmp2_MTX=Rcpp::as<Rcpp::NumericMatrix>(tmp_pu_S4.slot("coords"));
158 4
        pcv = &tmp2_MTX(0,0);
159 4
        Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>> tmp1_MTX(pcv, tmp2_MTX.nrow(), tmp2_MTX.ncol());
160 4
        species_space_pucoords_MTX.push_back(tmp1_MTX);
161

162
        if (verbose) Rcpp::Rcout << "\t\t\tpu points ids " << std::endl;
163 4
        Rcpp::IntegerVector tmpvec(tmp_pu_S4.slot("ids"));
164 4
        species_space_puids_RIV.push_back(tmpvec-1);
165

166
        if (verbose) Rcpp::Rcout << "\t\t\tnumber of planning units" << std::endl;
167 4
        species_space_npu_INT.push_back(tmpvec.size());
168
        
169
        // demand points
170
        if (verbose) Rcpp::Rcout << "\t\t\tdemand point coordinates" << std::endl;
171 4
        tmp2_MTX=Rcpp::as<Rcpp::NumericMatrix>(tmp_dp_S4.slot("coords"));
172 4
        pcv = &tmp2_MTX(0,0);
173 4
        Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>> tmp3_MTX(pcv, tmp2_MTX.nrow(), tmp2_MTX.ncol());
174 4
        species_space_dpcoords_MTX.push_back(tmp3_MTX);
175
        
176
        if (verbose) Rcpp::Rcout << "\t\t\tdemand point weights" << std::endl;
177 4
        Eigen::Map<Eigen::VectorXd> tmpvec2(Rcpp::as<Eigen::Map<Eigen::VectorXd>>(tmp_dp_S4.slot("weights")));
178 4
        species_space_dpweights_VXD.push_back(tmpvec2);
179
        
180
        if (verbose) Rcpp::Rcout << "\t\t\tnumber of demand points" << std::endl;
181 4
        species_space_ndp_INT.push_back(tmpvec2.size());
182
        
183
      }
184
    }
185
  }
186
  
187 4
  species_space_dpcoords_MTX.shrink_to_fit();
188 4
  species_space_pucoords_MTX.shrink_to_fit();
189 4
  species_space_dpweights_VXD.shrink_to_fit();
190 4
  species_space_puids_RIV.shrink_to_fit();
191 4
  species_space_ndp_INT.shrink_to_fit();
192 4
  species_space_npu_INT.shrink_to_fit();
193 4
  species_attributespace_species_INT.shrink_to_fit();
194 4
  species_attributespace_space_INT.shrink_to_fit();
195 4
  std::size_t n_species_attributespace_INT=species_attributespace_species_INT.size();
196
  
197
  // target data
198
  if (verbose) Rcpp::Rcout << "\ttarget data" << std::endl;
199 4
  Rcpp::DataFrame target_DF=Rcpp::as<Rcpp::DataFrame>(data.slot("targets"));
200 4
  std::vector<std::size_t> target_DF_species = target_DF["species"];
201 4
  std::vector<std::size_t> target_DF_target = target_DF["target"];
202 4
  std::vector<double> target_DF_value = target_DF["proportion"];
203 4
  std::vector<double> species_space_proptargets_DBL(n_species_attributespace_INT, NAN);
204 4
  std::vector<double> species_areatargets_DBL(n_species_INT, 0.0);
205 4
  for (std::size_t i=0; i<target_DF_species.size(); ++i) {
206 4
    if (target_DF_target[i]==0) {
207 4
      species_areatargets_DBL[target_DF_species[i]-1] = target_DF_value[i];
208
    } else {
209 4
      if (!std::isnan(target_DF_value[i])) {
210 4
        for (std::size_t a=0, ii=species_attributespace_species_INT[a], jj=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a, ii=species_attributespace_species_INT[a], jj=species_attributespace_space_INT[a]) {
211 4
          if ((ii==(target_DF_species[i]-1)) & (jj==(target_DF_target[i]-1))) {
212 4
            species_space_proptargets_DBL[a]=target_DF_value[i];
213 4
            break;
214
          }
215
        }
216
      }
217
    }
218
  }
219
  
220
  /// preprocessing for area targets
221
  if (verbose) Rcpp::Rcout << "\tpre-process area targets" << std::endl;
222
  // ini
223
  if (verbose) Rcpp::Rcout << "\t\tinitializing" << std::endl;
224
  double currArea;
225 4
  std::vector<double> species_totalarea_DBL(n_species_INT, 0.0);
226 4
  std::vector<Rcpp::IntegerVector> species_pu_ids_INT(n_species_INT);
227 4
  std::vector<std::vector<double>> species_pu_probs_DBL(n_species_INT);
228 4
  std::vector<std::size_t> species_npu_INT(n_species_INT);
229 4
  for (std::size_t i=0; i<n_species_INT; ++i) {
230 4
    species_pu_probs_DBL[i].reserve(n_pu_INT);
231
  }
232

233
  // calcs
234
  if (verbose) Rcpp::Rcout << "\t\tcalculate numbers" << std::endl;
235 4
  for (std::size_t i=0; i<puvspecies_DF_pu.size(); ++i) {
236
    // calculate area
237 4
    currArea=puvspecies_DF_value[i] * pu_DF_area[puvspecies_DF_pu[i]];
238
    // sum area occupied by species
239 4
    species_totalarea_DBL[puvspecies_DF_species[i]]+=currArea;
240
    // assign pu to species vector
241 4
    species_pu_ids_INT[puvspecies_DF_species[i]].push_back(puvspecies_DF_pu[i]);
242
    // assign prob to species vector
243 4
    species_pu_probs_DBL[puvspecies_DF_species[i]].push_back(puvspecies_DF_value[i]);
244
  }
245
  if (verbose) Rcpp::Rcout << "\t\tresize vectors" << std::endl;
246 4
  for (std::size_t i=0; i<n_species_INT; ++i) {
247 4
    species_pu_probs_DBL[i].shrink_to_fit();
248
  }
249
  
250
  /// create centroid variables
251
  if (verbose) Rcpp::Rcout << "\tcentroid positions" << std::endl;
252 4
  std::vector<Eigen::RowVectorXd> species_space_centroid_VXD(n_species_attributespace_INT);
253 4
  for (std::size_t a=0; a<n_species_attributespace_INT; ++a) {
254 4
    species_space_centroid_VXD[a]=species_space_dpcoords_MTX[a].colwise().sum();
255 4
    species_space_centroid_VXD[a]/=static_cast<double>(species_space_ndp_INT[a]);
256
  }
257

258
  if (verbose) Rcpp::Rcout << "\tdistance vars" << std::endl;
259 4
  std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> species_space_weightdist_MTX(n_species_attributespace_INT);
260 4
  std::vector<double> species_space_tss_DBL(n_species_attributespace_INT);
261
  {
262 4
    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp1_MTX;
263 4
    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp2_MTX;
264 4
    Eigen::ArrayXd tmp_AXD;
265
    double currFailDist;
266 4
    for (std::size_t a=0; a<n_species_attributespace_INT; ++a) {
267
      /// initialization
268
      /// preliminary processing
269 4
      if (unreliable_formulation) {
270 4
        species_space_weightdist_MTX[a].resize(species_space_ndp_INT[a], species_space_npu_INT[a]);
271
      } else {
272 4
        species_space_weightdist_MTX[a].resize(species_space_ndp_INT[a], species_space_npu_INT[a]+1);
273
      }
274 4
      species_space_weightdist_MTX[a].setZero();
275
      // calculate all distances between pus and dps
276 4
      for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
277 4
        for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
278 4
          tmp_AXD=species_space_pucoords_MTX[a].row(l) - species_space_dpcoords_MTX[a].row(k);
279 4
          species_space_weightdist_MTX[a](k,l) = tmp_AXD.square().sum() * species_space_dpweights_VXD[a][k];
280
        }
281
      }
282 4
      species_space_weightdist_MTX[a].array()+=zero_adjust; 
283
      // calculate distances for centroid
284 4
      species_space_tss_DBL[a]=zero_adjust;
285 4
      for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
286 4
        tmp_AXD=species_space_centroid_VXD[a] - species_space_dpcoords_MTX[a].row(k);
287 4
        species_space_tss_DBL[a]+=(tmp_AXD.square().sum() * species_space_dpweights_VXD[a][k]);
288
      }
289
      // failure pu if reliable formulation
290 4
      if (!unreliable_formulation) {
291
        // failure pu
292 4
        currFailDist=species_space_weightdist_MTX[a].maxCoeff()*failure_multiplier_DBL;
293 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
294 4
          species_space_weightdist_MTX[a](k,species_space_npu_INT[a])+=currFailDist;
295
        }
296
      }
297
    }
298
  }
299

300
  // resize vectors
301 4
  std::vector<std::size_t> species_space_rlevel_INT(n_species_attributespace_INT);
302 4
  if (!unreliable_formulation) {
303 4
    for (std::size_t a=0; a<n_species_attributespace_INT; ++a) {
304 4
      species_space_rlevel_INT[a]=std::min(maxrlevel_INT, species_space_npu_INT[a]-1);
305
    }
306
  }
307
  
308
  /// species space probs and transitional probabilities
309
  if (verbose) Rcpp::Rcout << "\ttransitional probs" << std::endl;
310 4
  std::vector<Rcpp::NumericVector> species_space_pu_probs_RDV(n_species_attributespace_INT);
311 4
  std::vector<Rcpp::NumericVector> species_space_pu_tprobs_RDV(n_species_attributespace_INT);
312 4
  std::vector<Rcpp::IntegerVector> species_space_pupos_RIV(n_species_attributespace_INT);
313 4
  if (!unreliable_formulation) {
314 4
    for (std::size_t a=0; a<n_species_attributespace_INT; ++a) {
315
      /// initialization
316 4
      species_space_pupos_RIV[a]=Rcpp::match(species_space_puids_RIV[a],species_pu_ids_INT[species_attributespace_species_INT[a]])-1;
317
      // add in real pus
318 4
      for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
319 4
        species_space_pu_probs_RDV[a].push_back(species_pu_probs_DBL[species_attributespace_species_INT[a]][species_space_pupos_RIV[a][l]]);
320 4
        species_space_pu_tprobs_RDV[a].push_back((1.0 - species_space_pu_probs_RDV[a][l]) / species_space_pu_probs_RDV[a][l]);
321
      }
322
      // add in failure pu
323 4
      species_space_pu_probs_RDV[a].push_back(1.0);
324
    }
325
  }
326
  
327
  // calculate targets
328
  if (verbose) Rcpp::Rcout << "\tspace targets" << std::endl;
329 4
  std::vector<double> species_space_rawtargets_DBL(n_species_attributespace_INT);
330 4
  for (std::size_t a=0; a<n_species_attributespace_INT; ++a)
331 4
    species_space_rawtargets_DBL[a] = (species_space_proptargets_DBL[a] - 1.0) * species_space_tss_DBL[a] * -1.0;
332
  
333
  /// create unordered map with variable names
334
  if (verbose) Rcout << "\tcreating undordered_maps with variable names" << std::endl;
335 4
  std::unordered_map<std::size_t, std::size_t> pu_MAP;
336 4
  std::unordered_map<std::pair<std::size_t,std::size_t>, std::size_t, boost::hash<std::pair<std::size_t,std::size_t>>> edge_MAP;
337
  std::array<std::size_t,4> curr_uARRAY;
338
  std::array<std::size_t,5> curr_rARRAY;
339 4
  std::unordered_map<std::array<std::size_t, 4>, std::size_t, boost::hash<std::array<std::size_t,4>>> uY_MAP;
340 4
  std::unordered_map<std::array<std::size_t, 5>, std::size_t, boost::hash<std::array<std::size_t,5>>> rY_MAP;
341 4
  std::unordered_map<std::array<std::size_t, 5>, std::size_t, boost::hash<std::array<std::size_t,5>>> rP_MAP;
342 4
  std::unordered_map<std::array<std::size_t, 5>, std::size_t, boost::hash<std::array<std::size_t,5>>> rW_MAP;
343
  
344
  // pu vars
345 4
  for (std::size_t i=0; i<n_pu_INT; ++i) {
346 4
    pu_MAP[i] = counter;
347 4
    ++counter;
348
  }
349

350
  // pu_pu boundary vars
351 4
  if (boundary) {
352 0
    for (std::size_t i=0; i<n_edges2_INT; ++i) {
353 0
      edge_MAP[edge_id_PAIR[i]] = counter;
354 0
      ++counter;
355
    }
356
  }
357
  
358
  // space vars
359 4
  if (unreliable_formulation) {
360 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
361 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
362 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
363 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
364 4
          for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
365
            // Y_var
366 4
            curr_uARRAY = {{i, j, k, l}};
367 4
            uY_MAP[curr_uARRAY] = counter;
368 4
            ++counter;
369
          }
370
        }
371
      }
372
    }
373
  } else {
374 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
375 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
376 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
377 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
378 4
          for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
379 4
            for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
380
              // create array
381 4
              curr_rARRAY = {{i , j, k, l, r}};
382
              
383
              // Y_var
384 4
              rY_MAP[curr_rARRAY] = counter;
385 4
              ++counter;
386

387
              // P_var
388 4
              rP_MAP[curr_rARRAY] = counter;
389 4
              ++counter;
390

391
              // W_var
392 4
              rW_MAP[curr_rARRAY] = counter;
393 4
              ++counter;
394
            }
395
          }
396
        }
397
      }
398
    }
399
  }
400
  
401
  //// Main processing
402
  if (verbose) Rcpp::Rcout << "Main processing" << std::endl;
403 4
  Rcpp::checkUserInterrupt();
404

405
  /// objective function
406 4
  Rcpp::checkUserInterrupt();
407
  if (verbose) Rcpp::Rcout << "\tobjective function" << std::endl;
408
  std::size_t problem_size_INT;
409 4
  if (unreliable_formulation) {
410 4
    problem_size_INT = pu_MAP.size() + edge_MAP.size() + uY_MAP.size();
411
  } else {
412 4
    problem_size_INT = pu_MAP.size() + edge_MAP.size() + rY_MAP.size() + rP_MAP.size() + rW_MAP.size();
413
  }
414 4
  std::vector<double> obj_DBL(problem_size_INT);
415

416
  // cost variables
417
  if (verbose) Rcpp::Rcout << "\t\tcost variables" << std::endl;
418 4
  for (std::size_t i=0; i<n_pu_INT; ++i) {
419 4
    obj_DBL[pu_MAP[i]] = pu_DF_cost[i];
420
  }
421

422
  // boundary variables
423
  if (verbose) Rcpp::Rcout << "\t\tboundary variables" << std::endl;
424 4
  if (boundary) {
425 0
    for (std::size_t i=0; i<n_edges2_INT; ++i) {
426 0
      obj_DBL[edge_MAP[edge_id_PAIR[i]]] = -(blm_DBL * boundary_DF_boundary[edge_pos_INT[i]]);
427
    }
428
  }
429

430
  /// constraints
431 4
  Rcpp::checkUserInterrupt();
432
  if (verbose) Rcpp::Rcout << "\tconstraints" << std::endl;
433 4
  std::vector<std::size_t> model_rows_INT;
434 4
  std::vector<std::size_t> model_cols_INT;
435 4
  std::vector<double> model_vals_DBL;
436 4
  std::vector<std::string> sense_STR;
437 4
  std::vector<double> rhs_DBL;
438 4
  counter=0;
439

440
  // preallocate variables storing model matrix
441 4
  preallocate_INT=n_pu_INT * n_species_INT * n_attribute_spaces_INT * 
442 4
    static_cast<std::size_t>(std::accumulate(species_space_ndp_INT.cbegin(), species_space_ndp_INT.cend(), 0));
443 4
  if (!unreliable_formulation) {
444 4
    preallocate_INT*=maxrlevel_INT;
445
  }
446 4
  model_rows_INT.reserve(preallocate_INT);
447 4
  model_cols_INT.reserve(preallocate_INT);
448 4
  model_vals_DBL.reserve(preallocate_INT);
449 4
  sense_STR.reserve(preallocate_INT);
450 4
  rhs_DBL.reserve(preallocate_INT);
451

452
  // area target constraints
453 4
  Rcpp::checkUserInterrupt();
454
  if (verbose) Rcpp::Rcout << "\t\tarea constraints" << std::endl;
455 4
  for (std::size_t i=0; i<n_species_INT; ++i) {
456 4
    for (std::size_t j=0; j<species_pu_ids_INT[i].size(); ++j) {
457 4
      model_rows_INT.push_back(counter);
458 4
      model_cols_INT.push_back(pu_MAP[species_pu_ids_INT[i][j]]);
459 4
      model_vals_DBL.push_back(species_pu_probs_DBL[i][j] * pu_DF_area[species_pu_ids_INT[i][j]]);
460
    }
461 4
    sense_STR.push_back(">=");
462 4
    rhs_DBL.push_back(species_totalarea_DBL[i] * species_areatargets_DBL[i]);
463 4
    ++counter;
464
  }
465

466
  // space target constraints
467 4
  Rcpp::checkUserInterrupt();
468
  if (verbose) Rcpp::Rcout << "\t\tspace constraints" << std::endl;
469 4
  if (unreliable_formulation) {
470 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
471 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
472 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
473 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
474 4
          for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
475 4
            curr_uARRAY={{i, j, k, l}};
476 4
            model_rows_INT.push_back(counter);
477 4
            model_cols_INT.push_back(uY_MAP[curr_uARRAY]);
478 4
            model_vals_DBL.push_back(species_space_weightdist_MTX[a](k,l));
479
          }
480
        }
481 4
        sense_STR.push_back("<=");
482 4
        rhs_DBL.push_back(species_space_rawtargets_DBL[a]);
483 4
        ++counter;
484
      }
485
    }
486
  } else {
487 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
488 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
489 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
490 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
491 4
          for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
492 4
            for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
493 4
              curr_rARRAY = {{i, j, k, l, r}};
494 4
              model_rows_INT.push_back(counter);
495 4
              model_cols_INT.push_back(rW_MAP[curr_rARRAY]);
496 4
              model_vals_DBL.push_back(species_space_weightdist_MTX[a](k,l));
497
            }
498
          }
499
        }
500 4
        sense_STR.push_back("<=");
501 4
        rhs_DBL.push_back(species_space_rawtargets_DBL[a]);
502 4
        ++counter;
503
      }
504
    }
505
  }
506
  
507
  // pu status constraints
508
  if (verbose) Rcpp::Rcout << "\t\tpu status constraints" << std::endl;
509 4
  for (std::size_t i=0; i<n_pu_INT; ++i) {
510 4
    if (pu_DF_status[i]==2) {
511
      // locked in
512 0
      model_rows_INT.push_back(counter);
513 0
      model_cols_INT.push_back(pu_MAP[i]);
514 0
      model_vals_DBL.push_back(1);
515 0
      sense_STR.push_back("=");
516 0
      rhs_DBL.push_back(1.0);
517 0
      ++counter;
518 4
    } else if (pu_DF_status[i]==3) {
519
      // locked  out
520 0
      model_rows_INT.push_back(counter);
521 0
      model_cols_INT.push_back(pu_MAP[i]);
522 0
      model_vals_DBL.push_back(1);
523 0
      sense_STR.push_back("=");
524 0
      rhs_DBL.push_back(0.0);
525 0
      ++counter;
526
    }
527
  }
528

529
  // boundary constraints
530 4
  Rcpp::checkUserInterrupt();
531
  if (verbose) Rcpp::Rcout << "\t\tboundary constraints" << std::endl;
532 4
  if (boundary) {
533 0
    for (std::size_t i=0; i<n_edges2_INT; ++i) {
534
      /// constraints if boundary_DF_id1[i] != boundary_DF_id2[i]
535
      // model+=boundary_DF_idpair_STR[i] + " - " + pu_DF_id_STR[boundary_DF_id1[i]] + " <= 0\n";
536 0
      model_rows_INT.push_back(counter);
537 0
      model_rows_INT.push_back(counter);
538 0
      model_cols_INT.push_back(edge_MAP[edge_id_PAIR[i]]);
539 0
      model_cols_INT.push_back(pu_MAP[boundary_DF_id1[edge_pos_INT[i]]]);
540 0
      model_vals_DBL.push_back(1.0);
541 0
      model_vals_DBL.push_back(-1.0);
542 0
      sense_STR.push_back("<=");
543 0
      rhs_DBL.push_back(0.0);
544 0
      ++counter;
545

546
      // model+=boundary_DF_idpair_STR[i] + " - " + pu_DF_id_STR[boundary_DF_id2[i]] + " <= 0\n";
547 0
      model_rows_INT.push_back(counter);
548 0
      model_rows_INT.push_back(counter);
549 0
      model_cols_INT.push_back(edge_MAP[edge_id_PAIR[i]]);
550 0
      model_cols_INT.push_back(pu_MAP[boundary_DF_id2[edge_pos_INT[i]]]);
551 0
      model_vals_DBL.push_back(1.0);
552 0
      model_vals_DBL.push_back(-1.0);
553 0
      sense_STR.push_back("<=");
554 0
      rhs_DBL.push_back(0.0);
555 0
      ++counter;
556

557
      // model+=boundary_DF_idpair_STR[i] + " - " +
558
        // pu_DF_id_STR[boundary_DF_id1[i]] + " - " +
559
        // pu_DF_id_STR[boundary_DF_id2[i]] + " >= -1\n";
560
      
561
      // constraints not strictly needed since decision variables are binary
562
//       model_rows_INT.push_back(counter);
563
//       model_rows_INT.push_back(counter);
564
//       model_rows_INT.push_back(counter);
565
//       model_cols_INT.push_back(variable_MAP[boundary_DF_idpair_STR[i]]);
566
//       model_cols_INT.push_back(variable_MAP[pu_DF_id_STR[boundary_DF_id1[edge_pos_INT[i]]]]);
567
//       model_cols_INT.push_back(variable_MAP[pu_DF_id_STR[boundary_DF_id2[edge_pos_INT[i]]]]);
568
//       model_vals_DBL.push_back(1);
569
//       model_vals_DBL.push_back(-1);
570
//       model_vals_DBL.push_back(-1);
571
//       sense_STR.push_back(">=");
572
//       rhs_DBL.push_back(-1.0);
573
//       ++counter;
574
    }
575
  }
576

577 4
  if (unreliable_formulation) {
578
    // 1b
579 4
    Rcpp::checkUserInterrupt();
580
    if (verbose) Rcpp::Rcout << "\t\teqn. 1b constraints" << std::endl;
581 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
582 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
583 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
584 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
585 4
          for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
586 4
            curr_uARRAY = {{i, j, k, l}};
587 4
            model_rows_INT.push_back(counter);
588 4
            model_cols_INT.push_back(uY_MAP[curr_uARRAY]);
589 4
            model_vals_DBL.push_back(1.0);
590
          }
591 4
          sense_STR.push_back("=");
592 4
          rhs_DBL.push_back(1.0);
593 4
          ++counter;
594
        }
595
      }
596
    }
597

598
    // 1c
599 4
    Rcpp::checkUserInterrupt();
600
    if (verbose) Rcpp::Rcout << "\t\teqn. 1c constraints" << std::endl;
601 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
602 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
603 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
604 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
605 4
          for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
606 4
            curr_uARRAY = {{i, j, k, l}};
607 4
            model_rows_INT.push_back(counter);
608 4
            model_cols_INT.push_back(uY_MAP[curr_uARRAY]);
609 4
            model_vals_DBL.push_back(1.0);
610

611 4
            model_rows_INT.push_back(counter);
612 4
            model_cols_INT.push_back(pu_MAP[species_space_puids_RIV[a][l]]);
613 4
            model_vals_DBL.push_back(-1.0);
614

615 4
            sense_STR.push_back("<=");
616 4
            rhs_DBL.push_back(0.0);
617 4
            ++counter;
618
          }
619
        }
620
      }
621
    }
622
  } else {
623
    if (verbose) Rcpp::Rcout << "\t\teqn. 1b constraints (rows " << counter << ", ";
624 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
625 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
626 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
627 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
628 4
          for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
629
            
630
            // original formulation
631
//               for (std::size_t l=0; l<(species_space_npu_INT(i,j)+1); ++l) {
632
//                 curr_STR="Y_"+int_STR[i]+"_"+int_STR[j]+"_"+int_STR[k]+"_"+int_STR[l]+"_"+int_STR[r];
633
//                 model_rows_INT.push_back(counter);
634
//                 model_cols_INT.push_back(variable_MAP[curr_STR]);
635
//                 model_vals_DBL.push_back(1.0);
636
//               }
637
//               for (std::size_t r2=0; r2<r; ++r2) {
638
//                 curr_STR="Y_"+int_STR[i]+"_"+int_STR[j]+"_"+int_STR[k]+"_"+int_STR[species_space_npu_INT(i,j)]+"_"+int_STR[r2];
639
//                 model_rows_INT.push_back(counter);
640
//                 model_cols_INT.push_back(variable_MAP[curr_STR]);
641
//                 model_vals_DBL.push_back(1.0);
642
//               }
643
//               sense_STR.push_back("=");
644
//               rhs_DBL.push_back(1.0);
645
//               ++counter;
646
//             }
647

648
            // force each R-level to be assigned to a pu
649 4
            for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
650 4
              curr_rARRAY = {{i, j, k, l, r}};
651 4
              model_rows_INT.push_back(counter);
652 4
              model_cols_INT.push_back(rY_MAP[curr_rARRAY]);
653 4
              model_vals_DBL.push_back(1.0);
654
            }
655 4
            sense_STR.push_back("=");
656 4
            rhs_DBL.push_back(1.0);
657 4
            ++counter;
658
          }
659
        }
660
      }
661
    }
662
    if (verbose) Rcpp::Rcout << counter << ")" << std::endl;
663
    
664
    // 1b extra - force each pu to be assigned to only 1 r-level
665
    if (verbose) Rcpp::Rcout << "\t\teqn. 1b (extra) constraints (rows " << counter << ", ";
666 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
667 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
668 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
669 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
670 4
          for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
671 4
            for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
672 4
              curr_rARRAY = {{i, j, k, l, r}};
673 4
              model_rows_INT.push_back(counter);
674 4
              model_cols_INT.push_back(rY_MAP[curr_rARRAY]);
675 4
              model_vals_DBL.push_back(1.0);
676
            }
677 4
            sense_STR.push_back("<=");
678 4
            rhs_DBL.push_back(1.0);
679 4
            ++counter;
680
          }
681
        }
682
      }
683
    }
684
    if (verbose) Rcpp::Rcout << counter << ")" << std::endl;
685
    
686
    // 1c
687 4
    Rcpp::checkUserInterrupt();
688
    if (verbose) Rcpp::Rcout << "\t\teqn. 1c constraints (rows " << counter << ", ";
689 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
690 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
691 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
692 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
693 4
          for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
694 4
            for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
695 4
              curr_rARRAY = {{i, j, k, l, r}};
696 4
              model_rows_INT.push_back(counter);
697 4
              model_cols_INT.push_back(rY_MAP[curr_rARRAY]);
698 4
              model_vals_DBL.push_back(1.0);
699
            }
700 4
            model_rows_INT.push_back(counter);
701 4
            model_cols_INT.push_back(pu_MAP[species_space_puids_RIV[a][l]]);
702 4
            model_vals_DBL.push_back(-1.0);
703 4
            sense_STR.push_back("<=");
704 4
            rhs_DBL.push_back(0.0);
705 4
            ++counter;
706
          }
707
        }
708
      }
709
    }
710
    if (verbose) Rcpp::Rcout << counter << ")" << std::endl;
711
    
712
    // 1d
713 4
    Rcpp::checkUserInterrupt();
714
    if (verbose) Rcpp::Rcout << "\t\teqn. 1d constraints (rows " << counter << ", ";
715 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
716 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
717 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
718 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
719
            
720
            // original formulation
721
  //             for (std::size_t r=0; r<(species_space_rlevel_INT(i,j)+1); ++r) {
722
  //               curr_STR="Y_"+int_STR[i]+"_"+int_STR[j]+"_"+int_STR[k]+"_"+int_STR[species_space_npu_INT(i,j)]+"_"+int_STR[r];
723
  //               model_rows_INT.push_back(counter);
724
  //               model_cols_INT.push_back(variable_MAP[curr_STR]);
725
  //               model_vals_DBL.push_back(1.0);
726
  //             }
727
  //             sense_STR.push_back("=");
728
  //             rhs_DBL.push_back(1.0);
729
  //             ++counter;
730
            
731
          // ensure that failure planning unit is assigned to last r-level
732 4
          curr_rARRAY = {{i, j, k, species_space_npu_INT[a], species_space_rlevel_INT[a]}};
733 4
          model_rows_INT.push_back(counter);
734 4
          model_cols_INT.push_back(rY_MAP[curr_rARRAY]);
735 4
          model_vals_DBL.push_back(1.0);
736 4
          sense_STR.push_back("=");
737 4
          rhs_DBL.push_back(1.0);
738 4
          ++counter;
739
        }
740
      }
741
    }
742
    if (verbose) Rcpp::Rcout << counter << ")" << std::endl;
743

744
    // 1e
745 4
    Rcpp::checkUserInterrupt();
746
    if (verbose) Rcpp::Rcout << "\t\teqn. 1e constraints (rows " << counter << ", ";
747 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
748 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
749 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
750 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
751 4
          for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
752
            // real pu
753 4
            curr_rARRAY = {{i, j, k, l, 0}};
754 4
            model_rows_INT.push_back(counter);
755 4
            model_cols_INT.push_back(rP_MAP[curr_rARRAY]);
756 4
            model_vals_DBL.push_back(1.0);
757

758 4
            sense_STR.push_back("=");
759 4
            rhs_DBL.push_back(species_space_pu_probs_RDV[a][l]);
760 4
            ++counter;
761
          }
762

763
          // failure pu
764 4
          curr_rARRAY= {{i, j, k, species_space_npu_INT[a], 0}};
765 4
          model_rows_INT.push_back(counter);
766 4
          model_cols_INT.push_back(rP_MAP[curr_rARRAY]);
767 4
          model_vals_DBL.push_back(1.0);
768

769 4
          sense_STR.push_back("=");
770 4
          rhs_DBL.push_back(1.0);
771 4
          ++counter;
772
        }
773
      }
774
    }
775
    if (verbose) Rcpp::Rcout << counter << ")" << std::endl;
776

777
    // 1f
778 4
    Rcpp::checkUserInterrupt();
779
    if (verbose) Rcpp::Rcout << "\t\teqn. 1f constraints (rows " << counter << ", ";
780 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
781 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
782 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
783 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
784 4
          for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
785 4
            for (std::size_t r=1, r2=0; r<(species_space_rlevel_INT[a]+1); ++r, ++r2) {
786 4
              curr_rARRAY = {{i, j, k, l, r}};
787 4
              model_rows_INT.push_back(counter);
788 4
              model_cols_INT.push_back(rP_MAP[curr_rARRAY]);
789 4
              model_vals_DBL.push_back(1.0);
790

791 4
              for (std::size_t l2=0; l2<species_space_npu_INT[a]; ++l2) {
792 4
                curr_rARRAY= {{i, j, k, l2, r2}};
793 4
                model_rows_INT.push_back(counter);
794 4
                model_cols_INT.push_back(rW_MAP[curr_rARRAY]);
795 4
                model_vals_DBL.push_back(-(species_space_pu_probs_RDV[a][l] * species_space_pu_tprobs_RDV[a][l2]));
796
              }
797

798 4
              sense_STR.push_back("=");
799 4
              rhs_DBL.push_back(0.0);
800 4
              ++counter;
801
            }
802
          }
803
        }
804
      }
805
    }
806
    if (verbose) Rcpp::Rcout << counter << ")" << std::endl;
807

808
    /// W variables
809 4
    Rcpp::checkUserInterrupt();
810
    if (verbose) Rcpp::Rcout << "\t\teqn 2. constraints (rows " << counter << ", ";
811 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
812 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
813 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
814 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
815 4
          for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
816 4
            for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
817
              // init
818 4
              curr_rARRAY = {{i, j, k, l, r}};
819
              // 2a
820
              //model+=currW + " - " + currP + " <= 0\n";
821 4
              model_rows_INT.push_back(counter);
822 4
              model_rows_INT.push_back(counter);
823 4
              model_cols_INT.push_back(rW_MAP[curr_rARRAY]);
824 4
              model_cols_INT.push_back(rP_MAP[curr_rARRAY]);
825 4
              model_vals_DBL.push_back(1.0);
826 4
              model_vals_DBL.push_back(-1.0);
827 4
              sense_STR.push_back("<=");
828 4
              rhs_DBL.push_back(0.0);
829 4
              ++counter;
830
              // 2b
831
              //model+=currW + " - " + currY + " <= 0\n";
832 4
              model_rows_INT.push_back(counter);
833 4
              model_rows_INT.push_back(counter);
834 4
              model_cols_INT.push_back(rW_MAP[curr_rARRAY]);
835 4
              model_cols_INT.push_back(rY_MAP[curr_rARRAY]);
836 4
              model_vals_DBL.push_back(1.0);
837 4
              model_vals_DBL.push_back(-1.0);
838 4
              sense_STR.push_back("<=");
839 4
              rhs_DBL.push_back(0.0);
840 4
              ++counter;
841
              // 2c
842
              //model+=currW + " >= 0\n";
843 4
              model_rows_INT.push_back(counter);
844 4
              model_cols_INT.push_back(rW_MAP[curr_rARRAY]);
845 4
              model_vals_DBL.push_back(1.0);
846 4
              sense_STR.push_back(">=");
847 4
              rhs_DBL.push_back(0.0);
848 4
              ++counter;
849
              // 2d
850
              //model+= currW + " - " + currP + " - " + currY + " >= -1\n";
851
              //model+= currW + " >= " + currP + " + " + currY + " - 1\n";
852 4
              model_rows_INT.push_back(counter);
853 4
              model_rows_INT.push_back(counter);
854 4
              model_rows_INT.push_back(counter);
855 4
              model_cols_INT.push_back(rW_MAP[curr_rARRAY]);
856 4
              model_cols_INT.push_back(rP_MAP[curr_rARRAY]);
857 4
              model_cols_INT.push_back(rY_MAP[curr_rARRAY]);
858 4
              model_vals_DBL.push_back(1.0);
859 4
              model_vals_DBL.push_back(-1.0);
860 4
              model_vals_DBL.push_back(-1.0);
861 4
              sense_STR.push_back(">=");
862 4
              rhs_DBL.push_back(-1.0);
863 4
              ++counter;
864
            }
865
          }
866
        }
867
      }
868
    }
869
    if (verbose) Rcpp::Rcout << counter << ")" << std::endl;
870
  }
871
  
872
  ////// variable types and bounds
873 4
  Rcpp::checkUserInterrupt();
874
  if (verbose) Rcpp::Rcout << "\tvariable types" << std::endl;
875 4
  std::vector<std::string> vtype_STR(obj_DBL.size());
876 4
  std::vector<std::size_t> lb_INT(obj_DBL.size(), 0);
877 4
  std::vector<std::size_t> ub_INT(obj_DBL.size(), 1);
878

879
  ///// binary
880
  if (verbose) Rcpp::Rcout << "\t\t binary vars" << std::endl;
881

882
  // 1g
883 4
  Rcpp::checkUserInterrupt();
884 4
  for (std::size_t i=0; i<n_pu_INT; ++i) {
885 4
    vtype_STR[pu_MAP[i]]="B";
886
  }
887

888
  // boundary variables
889 4
  Rcpp::checkUserInterrupt();
890 4
  if (boundary) {
891 0
    for (std::size_t i=0; i<n_edges2_INT; ++i)
892 0
      vtype_STR[edge_MAP[edge_id_PAIR[i]]]="B";
893
  }
894

895
  // 1g constraints
896 4
  Rcpp::checkUserInterrupt();
897 4
  if (unreliable_formulation) {
898
    // 1g
899 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
900 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
901 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
902 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
903 4
          for (std::size_t l=0; l<species_space_npu_INT[a]; ++l) {
904 4
            curr_uARRAY = {{i, j, k, l}};
905 4
            vtype_STR[uY_MAP[curr_uARRAY]]="B";
906
          }
907
        }
908
      }
909
    }
910
  } else {
911
    // 1g
912 4
    Rcpp::checkUserInterrupt();
913 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
914 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
915 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
916 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
917 4
          for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
918 4
            for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
919 4
              curr_rARRAY={{i, j, k, l, r}};
920 4
              vtype_STR[rY_MAP[curr_rARRAY]]="B";
921
            }
922
          }
923
        }
924
      }
925
    }
926
    
927
    ///// semi-continuous
928 4
    Rcpp::checkUserInterrupt();
929
    if (verbose) Rcpp::Rcout << "\t\tsemi-continuous vars" << std::endl;
930 4
    for (std::size_t a=0, i=species_attributespace_species_INT[a], j=species_attributespace_space_INT[a]; a<n_species_attributespace_INT; ++a) {
931 4
      i=species_attributespace_species_INT[a]; j=species_attributespace_space_INT[a];
932 4
      if (!std::isnan(species_space_proptargets_DBL[a])) {
933 4
        for (std::size_t k=0; k<species_space_ndp_INT[a]; ++k) {
934 4
          for (std::size_t l=0; l<(species_space_npu_INT[a]+1); ++l) {
935 4
            for (std::size_t r=0; r<(species_space_rlevel_INT[a]+1); ++r) {
936 4
              curr_rARRAY = {{i, j, k, l, r}};
937
              // w variables
938 4
              vtype_STR[rW_MAP[curr_rARRAY]]="S";
939
              // p variables
940 4
              vtype_STR[rP_MAP[curr_rARRAY]]="S";
941
            }
942
          }
943
        }
944
      }
945
    }
946
  }
947
  
948
  // calculate best space proportions
949
  if (verbose) Rcpp::Rcout << "\tcalculating best possible space targets" << std::endl;
950 4
  std::vector<double> species_space_best_DBL(n_species_attributespace_INT);
951 4
  if (unreliable_formulation) {
952 4
    for (std::size_t a=0; a<n_species_attributespace_INT; ++a) {
953 4
      species_space_best_DBL[a]=1.0-(unreliable_space_value(species_space_weightdist_MTX[a]) / species_space_tss_DBL[a]);
954
    }
955
  } else {
956 4
    for (std::size_t a=0; a<n_species_attributespace_INT; ++a)
957 4
      species_space_best_DBL[a]=1.0-(reliable_space_value(species_space_weightdist_MTX[a],species_space_pu_probs_RDV[a],species_space_rlevel_INT[a]) / species_space_tss_DBL[a]);
958
  }
959
  
960
  // dump all variables to string
961 4
  std::vector<std::string> variables_STR(problem_size_INT);
962 4
  for (auto i : pu_MAP)
963 4
    variables_STR[i.second] = i.first;
964
  for (auto i : edge_MAP)
965 0
    variables_STR[i.second] = "pu_" + num2str<std::size_t>(i.first.first) + "_" + num2str<std::size_t>(i.first.second);
966 4
  if (unreliable_formulation) {
967 4
    for (auto i : uY_MAP)
968 4
      variables_STR[i.second] = "Y_" + num2str<std::size_t>(i.first[0]) + "_" + num2str<std::size_t>(i.first[1]) + "_" + num2str<std::size_t>(i.first[2]) + "_" + num2str<std::size_t>(i.first[3]);
969
  } else {
970 4
    for (auto i : rY_MAP)
971 4
      variables_STR[i.second] = "Y_" + num2str<std::size_t>(i.first[0]) + "_" + num2str<std::size_t>(i.first[1]) + "_" + num2str<std::size_t>(i.first[2]) + "_" + num2str<std::size_t>(i.first[3]) + "_" + num2str<std::size_t>(i.first[4]);
972 4
    for (auto i : rP_MAP)
973 4
      variables_STR[i.second] = "P_" + num2str<std::size_t>(i.first[0]) + "_" + num2str<std::size_t>(i.first[1]) + "_" + num2str<std::size_t>(i.first[2]) + "_" + num2str<std::size_t>(i.first[3]) + "_" + num2str<std::size_t>(i.first[4]);
974 4
    for (auto i : rW_MAP)
975 4
      variables_STR[i.second] = "W_" + num2str<std::size_t>(i.first[0]) + "_" + num2str<std::size_t>(i.first[1]) + "_" + num2str<std::size_t>(i.first[2]) + "_" + num2str<std::size_t>(i.first[3]) + "_" + num2str<std::size_t>(i.first[4]);
976
  }
977
  
978
  // create pointers to store cache
979 4
  std::vector<std::size_t>* species_attributespace_species_INT_PTR = new std::vector<std::size_t>;
980 4
  std::vector<std::size_t>* species_attributespace_space_INT_PTR = new std::vector<std::size_t>;
981 4
  std::vector<Rcpp::NumericVector>* species_space_pu_probs_RDV_PTR = new std::vector<Rcpp::NumericVector>;
982 4
  std::vector<Rcpp::IntegerVector>* species_space_puids_RIV_PTR = new std::vector<Rcpp::IntegerVector>;
983 4
  std::vector<Rcpp::IntegerVector>* species_space_pupos_RIV_PTR = new std::vector<Rcpp::IntegerVector>;
984 4
  std::vector<double>* species_areatargets_DBL_PTR = new std::vector<double>;
985 4
  std::vector<double>* species_totalarea_DBL_PTR = new std::vector<double>;
986 4
  std::vector<double>* species_space_rawtargets_DBL_PTR = new std::vector<double>;
987 4
  std::vector<double>* species_space_tss_DBL_PTR = new std::vector<double>;
988 4
  std::vector<double>* species_space_best_DBL_PTR = new std::vector<double>;
989 4
  std::vector<double>* species_space_proptargets_DBL_PTR = new std::vector<double>;
990 4
  std::vector<std::size_t>* species_space_rlevel_INT_PTR = new std::vector<std::size_t>;
991 4
  std::vector<std::string>* variables_STR_PTR = new std::vector<std::string>;
992 4
  std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>* species_space_weightdist_MTX_PTR  = new std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>;
993
  
994
  // copy data for preservation
995 4
  *species_attributespace_species_INT_PTR = species_attributespace_species_INT;
996 4
  *species_attributespace_space_INT_PTR = species_attributespace_space_INT;
997 4
  *species_space_puids_RIV_PTR = species_space_puids_RIV;
998 4
  *species_space_pu_probs_RDV_PTR = species_space_pu_probs_RDV;
999 4
  *species_space_pupos_RIV_PTR = species_space_pupos_RIV;
1000 4
  *species_areatargets_DBL_PTR = species_areatargets_DBL;
1001 4
  *species_totalarea_DBL_PTR = species_totalarea_DBL;
1002 4
  *species_space_rawtargets_DBL_PTR = species_space_rawtargets_DBL;
1003 4
  *species_space_tss_DBL_PTR = species_space_tss_DBL;
1004 4
  *species_space_best_DBL_PTR = species_space_best_DBL;
1005 4
  *species_space_proptargets_DBL_PTR = species_space_proptargets_DBL;
1006 4
  *species_space_rlevel_INT_PTR = species_space_rlevel_INT;
1007 4
  *variables_STR_PTR = variables_STR;
1008 4
  *species_space_weightdist_MTX_PTR = species_space_weightdist_MTX;
1009
  
1010
  // convert pointers to Rcpp::XPters
1011 4
  Rcpp::XPtr<std::vector<std::size_t>> species_attributespace_species_INT_XPTR(species_attributespace_species_INT_PTR);
1012 4
  Rcpp::XPtr<std::vector<std::size_t>> species_attributespace_space_INT_XPTR(species_attributespace_space_INT_PTR);
1013 4
  Rcpp::XPtr<std::vector<Rcpp::NumericVector>> species_space_pu_probs_RDV_XPTR(species_space_pu_probs_RDV_PTR);
1014 4
  Rcpp::XPtr<std::vector<Rcpp::IntegerVector>> species_space_puids_RIV_XPTR(species_space_puids_RIV_PTR);
1015 4
  Rcpp::XPtr<std::vector<Rcpp::IntegerVector>> species_space_pupos_RIV_XPTR(species_space_pupos_RIV_PTR);
1016 4
  Rcpp::XPtr<std::vector<double>> species_areatargets_DBL_XPTR(species_areatargets_DBL_PTR);
1017 4
  Rcpp::XPtr<std::vector<double>> species_totalarea_DBL_XPTR(species_totalarea_DBL_PTR);
1018 4
  Rcpp::XPtr<std::vector<double>> species_space_rawtargets_DBL_XPTR(species_space_rawtargets_DBL_PTR);
1019 4
  Rcpp::XPtr<std::vector<double>> species_space_tss_DBL_XPTR(species_space_tss_DBL_PTR);
1020 4
  Rcpp::XPtr<std::vector<double>> species_space_best_DBL_XPTR(species_space_best_DBL_PTR);
1021 4
  Rcpp::XPtr<std::vector<double>> species_space_proptargets_DBL_XPTR(species_space_proptargets_DBL_PTR);
1022 4
  Rcpp::XPtr<std::vector<std::size_t>> species_space_rlevel_INT_XPTR(species_space_rlevel_INT_PTR);
1023 4
  Rcpp::XPtr<std::vector<std::string>> variables_STR_XPTR(variables_STR_PTR);
1024 4
  Rcpp::XPtr<std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>> species_space_weightdist_MTX_XPTR(species_space_weightdist_MTX_PTR);
1025
  
1026
  //// Exports
1027 4
  Rcpp::checkUserInterrupt();
1028
  if (verbose) Rcpp::Rcout << "Sending data to R" << std::endl;
1029
  return(
1030
    Rcpp::List::create(
1031 4
      Rcpp::Named("Ar") = Rcpp::List::create(
1032 4
        Rcpp::Named("row") = Rcpp::wrap(model_rows_INT),
1033 4
        Rcpp::Named("col") = Rcpp::wrap(model_cols_INT),
1034 4
        Rcpp::Named("value") = Rcpp::wrap(model_vals_DBL)
1035
      ),
1036 4
      Rcpp::Named("obj") = Rcpp::wrap(obj_DBL),
1037 4
      Rcpp::Named("sense") = Rcpp::wrap(sense_STR),
1038 4
      Rcpp::Named("rhs") = Rcpp::wrap(rhs_DBL),
1039 4
      Rcpp::Named("vtype") = Rcpp::wrap(vtype_STR),
1040 4
      Rcpp::Named("lb") = Rcpp::wrap(lb_INT),
1041 4
      Rcpp::Named("ub") = Rcpp::wrap(ub_INT),
1042 4
      Rcpp::Named("cache") = Rcpp::List::create(
1043 4
        Rcpp::Named("species_attributespace_species_INT") = species_attributespace_species_INT_XPTR,
1044 4
        Rcpp::Named("species_attributespace_space_INT") = species_attributespace_space_INT_XPTR,
1045 4
        Rcpp::Named("species_space_pu_probs_RDV") = species_space_pu_probs_RDV_XPTR,
1046 4
        Rcpp::Named("species_space_puids_RIV") = species_space_puids_RIV_XPTR,
1047 4
        Rcpp::Named("species_space_pupos_RIV") = species_space_pupos_RIV_XPTR,
1048 4
        Rcpp::Named("species_areatargets_DBL") = species_areatargets_DBL_XPTR,
1049 4
        Rcpp::Named("species_totalarea_DBL") = species_totalarea_DBL_XPTR,
1050 4
        Rcpp::Named("species_space_rawtargets_DBL") = species_space_rawtargets_DBL_XPTR,
1051 4
        Rcpp::Named("species_space_tss_DBL") = species_space_tss_DBL_XPTR,
1052 4
        Rcpp::Named("species_space_best_DBL") = species_space_best_DBL_XPTR,
1053 4
        Rcpp::Named("species_space_proptargets_DBL") = species_space_proptargets_DBL_XPTR,
1054 4
        Rcpp::Named("species_space_rlevel_INT") = species_space_rlevel_INT_XPTR,
1055 4
        Rcpp::Named("variables") = variables_STR_XPTR,
1056 4
        Rcpp::Named("wdist") = species_space_weightdist_MTX_XPTR
1057
      ),
1058 4
      Rcpp::Named("modelsense") = "min"
1059
    )
1060
  );
1061
}

Read our documentation on viewing source code .

Loading