@@ -11,7 +11,7 @@ vignette: >
1111
1212
1313
14- ``` {r setup, include =FALSE,eval=TRUE}
14+ ``` {r setup, echo =FALSE, include=TRUE, eval=TRUE}
1515knitr::opts_chunk$set(echo = TRUE
1616 ,eval = TRUE
1717 ,cache = TRUE
@@ -23,11 +23,113 @@ library(rlist)
2323library(data.table)
2424library(foreach)
2525
26- source("gravity_model_functions.R")
26+ ```
27+
28+ ``` {r, echo=FALSE}
29+
30+ # The good old original R code to run SIM that was used before cppSim
31+
32+ ### this file contains the functions used to run the gravity model
33+
34+ cost_function <- function(d, beta, type = "exp") {
35+
36+ # d is the od distance matrix, taking the exponential
37+ # means doing this operation for each individual element
38+ # with an exponent beta
39+ # the type parameter allows to set it to either exp or pow.
40+ # pow means we use a power function as cost, rather than exponential
41+
42+ if (type == "exp") {
43+ exp(-beta*d)
44+ } else if (type == "pow") {
45+ d^(-beta)
46+ } else {
47+ print("provide a type of functino to compute")
48+ }
49+ }
50+
51+ r_2 <- function(d, f) cor(d |> as.numeric()
52+ ,f |> as.numeric())^2
53+
54+ r <- function(d, f) cor(d |> as.numeric()
55+ ,f |> as.numeric())
56+
57+ rmse <- function(d,f) sum((d-f)^2)
58+
59+ calibration <- function(cost_fun,O,D,delta = 0.05) {
60+ B <- rep_len(1,nrow(cost_fun))
61+ eps <- abs(sum(B))
62+ e <- NULL
63+ i <- 0
64+ while((eps > delta) & (i<50)) {
65+ A_new <- 1/(apply(cost_fun,function (x) sum(B*D*x),MARGIN = 1))
66+ B_new <- 1/(apply(cost_fun,function (x) sum(A_new*O*x),MARGIN = 2))
67+ eps <- abs(sum(B_new-B))
68+ e <- append(e,eps)
69+ A <- A_new
70+ B <- B_new
71+ i <- i+1
72+ }
73+ list(
74+ "A"= A
75+ ,"B" = B
76+ ,"e" = e
77+ )
78+ }
79+
80+ run_model <- function(flows
81+ ,distance
82+ ,beta = 0.25
83+ ,type = "exp"
84+ #,cores = 3
85+ ) {
86+
87+ F_c <- cost_function(d = {{distance}},beta = {{beta}},type = type)
88+ print("cost function computed")
89+ O <- apply(flows,sum, MARGIN = 1) |> as.integer()
90+ D <- apply(flows,sum, MARGIN = 2) |> as.integer()
91+ A_B <- calibration(cost_fun = F_c
92+ ,O=O
93+ ,D=D
94+ ,delta = .001)
95+ print("calibration: over")
96+ A <- A_B$A
97+ B <- A_B$B
98+
99+ flows_model <- foreach(j = c(1:nrow(F_c))
100+ ,.combine = rbind) %do% {
101+ round(A[j]*B*O[j]*D*F_c[j,])
102+ }
103+
104+ print("model run: over")
105+ e_sor <- e_sorensen(flows,flows_model) |> as.numeric()
106+ print(paste0("E_sor = ",e_sor))
107+ r2 <- r_2(flows_model,flows) |> as.numeric()
108+ print(paste0("r2 = ",r2))
109+ RMSE <- rmse(flows_model,flows) |> as.numeric()
110+ print(paste0("RMSE = ",RMSE))
111+
112+ list("values" = flows_model
113+ ,"r2" = r2
114+ ,"rmse" = RMSE
115+ ,"calib" = A_B$e
116+ ,"e_sor" = e_sor
117+ )
118+
119+ }
120+
121+ ### Validation
122+
123+ e_sorensen <- function(data, fit) {
124+ 2*sum(apply(cbind(data |> c()
125+ ,fit |> c()), MARGIN = 1, FUN = min))/(sum(data) + sum(fit))
126+ }
127+
27128
28129```
29130
30131
132+
31133``` {r load_data}
32134
33135flows_matrix <- list.load("flows_matrix.rds")
@@ -59,9 +161,9 @@ with the derivation made from a recursive chain with initial values 1. Let's ref
59161
60162``` {r}
61163# creating the O, D vectors.
62- O <- apply(flows_matrix,reduce , MARGIN = 2, sum ) |> c()
164+ O <- apply(flows_matrix, sum , MARGIN = 2) |> c()
63165
64- D <- apply(flows_matrix,reduce , MARGIN = 1, sum ) |> c()
166+ D <- apply(flows_matrix, sum , MARGIN = 1) |> c()
65167```
66168
67169
0 commit comments