From 1c5a2b8571bf5a932cd399b7c64bc66c04a15acf Mon Sep 17 00:00:00 2001 From: akone42 Date: Wed, 28 May 2025 16:24:58 -0400 Subject: [PATCH 1/9] City Bike command line --- week1/citibike.sh | 43 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 40 insertions(+), 3 deletions(-) diff --git a/week1/citibike.sh b/week1/citibike.sh index 25604f545..71e8b7331 100755 --- a/week1/citibike.sh +++ b/week1/citibike.sh @@ -4,15 +4,52 @@ # # count the number of unique stations - +# $ cut -d, -f4 201402-citibike-tripdata.csv | sort -u | wc -l + # 330 # count the number of unique bikes - +# $ cut -d, -f12 201402-citibike-tripdata.csv | sort -u | wc -l + #5700 # count the number of trips per day +# $ cut -d, -f2 201402-citibike-tripdata.csv | cut -d" " -f1 | sort | uniq -c +# 12771 2014-02-01 +# 13816 2014-02-02 +# 2600 2014-02-03 +# 8709 2014-02-04 +# 2746 2014-02-05 +# 7196 2014-02-06 +# 8495 2014-02-07 +# 5986 2014-02-08 +# 4996 2014-02-09 +# 6846 2014-02-10 +# 8343 2014-02-11 +# 8580 2014-02-12 +# 876 2014-02-13 +# 3609 2014-02-14 +# 2261 2014-02-15 +# 3003 2014-02-16 +# 4854 2014-02-17 +# 5140 2014-02-18 +# 8506 2014-02-19 +# 11792 2014-02-20 +# 8680 2014-02-21 +# 13044 2014-02-22 +# 13324 2014-02-23 +# 12922 2014-02-24 +# 12830 2014-02-25 +# 11188 2014-02-26 +# 12036 2014-02-27 +# 9587 2014-02-28 + # find the day with the most rides -# find the day with the fewest rides +# $ cut -d, -f2 201402-citibike-tripdata.csv | cut -d" " -f1 | sort | uniq -c| sort -r | head -n1 +# 13816 2014-02-02 +# find the day with the fewest rides +# $ cut -d, -f2 201402-citibike-tripdata.csv | cut -d" " -f1| sort | uniq -c | sort | head -n2 +# 1 starttime +# 876 2014-02-13 # find the id of the bike with the most rides # count the number of rides by gender and birth year From b288a4b7ee65b8eb8dd558ccb8d6fa9e3568fe38 Mon Sep 17 00:00:00 2001 From: akone42 Date: Wed, 28 May 2025 16:39:55 -0400 Subject: [PATCH 2/9] City Bike More updates --- week1/citibike.sh | 150 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) diff --git a/week1/citibike.sh b/week1/citibike.sh index 71e8b7331..d9ea9609a 100755 --- a/week1/citibike.sh +++ b/week1/citibike.sh @@ -51,10 +51,160 @@ # 1 starttime # 876 2014-02-13 # find the id of the bike with the most rides +# $ cut -d, -f12 201402-citibike-tripdata.csv | sort | uniq -c | head -n1 +# 10 14529 # count the number of rides by gender and birth year +# $ cut -d, -f15,14 201402-citibike-tripdata.csv | sort | uniq -c +# 6717 \N,0 +# 9 1899,1 +# 68 1900,1 +# 11 1901,1 +# 5 1907,1 +# 4 1910,1 +# 1 1913,1 +# 3 1917,1 +# 1 1921,1 +# 32 1922,1 +# 5 1926,2 +# 2 1927,1 +# 1 1932,1 +# 7 1932,2 +# 10 1933,1 +# 21 1934,1 +# 14 1935,1 +# 31 1936,1 +# 24 1937,1 +# 70 1938,1 +# 5 1938,2 +# 24 1939,1 +# 19 1939,2 +# 83 1940,1 +# 1 1940,2 +# 148 1941,1 +# 16 1941,2 +# 173 1942,1 +# 9 1942,2 +# 108 1943,1 +# 22 1943,2 +# 277 1944,1 +# 34 1944,2 +# 171 1945,1 +# 43 1945,2 +# 424 1946,1 +# 30 1946,2 +# 391 1947,1 +# 60 1947,2 +# 664 1948,1 +# 143 1948,2 +# 624 1949,1 +# 101 1949,2 +# 738 1950,1 +# 152 1950,2 +# 6 1951,0 +# 1006 1951,1 +# 146 1951,2 +# 1040 1952,1 +# 143 1952,2 +# 1474 1953,1 +# 301 1953,2 +# 1636 1954,1 +# 306 1954,2 +# 1568 1955,1 +# 349 1955,2 +# 1777 1956,1 +# 542 1956,2 +# 1676 1957,1 +# 562 1957,2 +# 2333 1958,1 +# 643 1958,2 +# 2281 1959,1 +# 539 1959,2 +# 2679 1960,1 +# 776 1960,2 +# 2315 1961,1 +# 432 1961,2 +# 2808 1962,1 +# 833 1962,2 +# 3514 1963,1 +# 715 1963,2 +# 3679 1964,1 +# 570 1964,2 +# 2957 1965,1 +# 687 1965,2 +# 3440 1966,1 +# 565 1966,2 +# 4016 1967,1 +# 634 1967,2 +# 3931 1968,1 +# 545 1968,2 +# 4557 1969,1 +# 898 1969,2 +# 4657 1970,1 +# 1079 1970,2 +# 4132 1971,1 +# 791 1971,2 +# 4066 1972,1 +# 962 1972,2 +# 4097 1973,1 +# 877 1973,2 +# 4957 1974,1 +# 891 1974,2 +# 4185 1975,1 +# 699 1975,2 +# 4557 1976,1 +# 1022 1976,2 +# 4817 1977,1 +# 1140 1977,2 +# 5645 1978,1 +# 1231 1978,2 +# 6433 1979,1 +# 1338 1979,2 +# 6173 1980,1 +# 1488 1980,2 +# 6620 1981,1 +# 1588 1981,2 +# 6244 1982,1 +# 1724 1982,2 +# 6890 1983,1 +# 1889 1983,2 +# 7348 1984,1 +# 1791 1984,2 +# 7043 1985,1 +# 2262 1985,2 +# 6147 1986,1 +# 1962 1986,2 +# 5776 1987,1 +# 1696 1987,2 +# 6449 1988,1 +# 1599 1988,2 +# 5408 1989,1 +# 1435 1989,2 +# 4541 1990,1 +# 1156 1990,2 +# 8 1991,0 +# 2377 1991,1 +# 689 1991,2 +# 1758 1992,1 +# 410 1992,2 +# 1398 1993,1 +# 289 1993,2 +# 927 1994,1 +# 288 1994,2 +# 664 1995,1 +# 163 1995,2 +# 234 1996,1 +# 100 1996,2 +# 164 1997,1 +# 87 1997,2 +# 1 birth year,gender # count the number of trips that start on cross streets that both contain numbers (e.g., "1 Ave & E 15 St", "E 39 St & 2 Ave", ...) +# $ cut -d, -f5 201402-citibike-tripdata.csv | grep '[1-9].*[1-9]'| wc -l +# 122406 # compute the average trip duration + +# $ awk -F, '{sum+=$1; count++} END {print sum/count}' 201402-citibike-tripdata.csv +# 874.516 \ No newline at end of file From 1915f6489bfbcac7b3cb4e381392d999feb06120 Mon Sep 17 00:00:00 2001 From: akone42 Date: Fri, 30 May 2025 10:20:13 -0400 Subject: [PATCH 3/9] Day 3 --- week1/citibike.R | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/week1/citibike.R b/week1/citibike.R index ad01de1d3..110b7adf4 100644 --- a/week1/citibike.R +++ b/week1/citibike.R @@ -6,7 +6,7 @@ library(lubridate) ######################################## # read one month of data -trips <- read_csv('201402-citibike-tripdata.csv') +trips <- read_csv('./week1/201402-citibike-tripdata.csv') # replace spaces in column names with underscores names(trips) <- gsub(' ', '_', names(trips)) @@ -23,26 +23,54 @@ trips <- mutate(trips, gender = factor(gender, levels=c(0,1,2), labels = c("Unkn ######################################## # count the number of trips (= rows in the data frame) +print(nrow(trips)) + # find the earliest and latest birth years (see help for max and min to deal with NAs) + birth_year_col <- as.numeric(birth_year_col) + birth_year_col <- na.omit(birth_year_col) + min(birth_year_col) + max(birth_year_col) # use filter and grepl to find all trips that either start or end on broadway +filter(trips, grepl("Broadway", start_station_name) | grepl("Broadway", end_station_name)) # do the same, but find all trips that both start and end on broadway +filter(trips, grepl("Broadway", start_station_name) & grepl("Broadway", end_station_name)) # find all unique station names + stations <- paste(trips$start_station_name, trips$end_station_name) + unique_stations <- unique(stations) + # count the number of trips by gender, the average trip time by gender, and the standard deviation in trip time by gender # do this all at once, by using summarize() with multiple arguments +trips %>% group_by(gender) %>% summarize(count =n(), mean = mean(tripduration), std = sd(tripduration)) + + # find the 10 most frequent station-to-station trips +trips %>% +group_by(start_station_name, end_station_name) +%>% summarize(frequency = n()) %>% arrange(desc(frequency)) %>% slice_head(n = 10) # find the top 3 end stations for trips starting from each start station - + trips %>% group_by(start_station_name, end_station_name) %>% summarize(count = n()) %>% group_by(start_station_name) %>% arrange(desc(count)) %>% slice_head(n=3) # find the top 3 most common station-to-station trips by gender - +trips %>% group_by(start_station_name, end_station_name) %>% +group_by(gender) %>% summarize(frequency = n()) %>% arrange(desc(frequency)) %>% + slice_head(n =3) # find the day with the most trips # tip: first add a column for year/month/day without time of day (use as.Date or floor_date from the lubridate package) +trips %>% group_by(as.Date(starttime)) %>% summarize(frequency = n()) %>% ar$ +slice_head(n =1) + # compute the average number of trips taken during each of the 24 hours of the day across the entire month # what time(s) of day tend to be peak hour(s)? + +trips %>% mutate( hours = hour(starttime)) %>% + group_by(hours) %>% summarise(count = n(), day_in_month_count = days_in_month(starttime), avg = count / day_in_month_count) + +trips %>% mutate( hours = hour(starttime)) %>% + group_by(hours) %>% summarise(count = n()) %>% arrange(desc(count)) %>% slice(1) From 448669f99b428c3814d166df3cf3018e5f3e7dc2 Mon Sep 17 00:00:00 2001 From: akone42 Date: Mon, 2 Jun 2025 10:33:28 -0400 Subject: [PATCH 4/9] Day 4 plotting --- week1/plot_trips.R | 38 +++++++++++++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/week1/plot_trips.R b/week1/plot_trips.R index 4f25437ba..9829d088b 100644 --- a/week1/plot_trips.R +++ b/week1/plot_trips.R @@ -11,20 +11,48 @@ theme_set(theme_bw()) # load RData file output by load_trips.R load('trips.RData') - - ######################################## # plot trip data ######################################## # plot the distribution of trip times across all rides (compare a histogram vs. a density plot) - +ggplot(trips, aes(x = tripduration))+ + geom_histogram(fill = 'blue', bins = 30) + + scale_x_log10(labels=comma) + + ylab('Number of trip per Time') + + ggplot(trips, aes(x = tripduration))+ + geom_density(fill = 'grey') + + scale_x_log10(labels=comma) # plot the distribution of trip times by rider type indicated using color and fill (compare a histogram vs. a density plot) +ggplot(trips, aes(x = tripduration, , color = usertype, fill = usertype))+ + geom_histogram(bins = 30) + + scale_x_log10(labels=comma) + + ylab('Number of trip per Time')+ + facet_grid( ~usertype) + +ggplot(trips, aes(x = tripduration, , color = usertype, fill = usertype))+ + geom_density(fill = 'grey') + + scale_x_log10(labels=comma) + + ylab('Number of trip per Time')+ + facet_grid( ~usertype) + # plot the total number of trips on each day in the dataset +trips %>% mutate(date = as.Date(starttime)) %>% +ggplot(aes(x = date))+ + geom_histogram(fill = 'blue', bins = 30)+ + scale_y_continuous(label = comma)+ + scale_x_date(date_labels = "%Y-%m-%d") + # plot the total number of trips (on the y axis) by age (on the x axis) and gender (indicated with color) - +trips %>% mutate(age = sys.Date() - birth_year) %>% +ggplot(aes(x = age, color = usertype, fill = usertype))+ + geom_histogram(bins = 30) + + scale_x_log10(labels=comma) + + ylab()+ + facet_grid(~usertype) # plot the ratio of male to female trips (on the y axis) by age (on the x axis) # hint: use the pivot_wider() function to reshape things to make it easier to compute this ratio # (you can skip this and come back to it tomorrow if we haven't covered pivot_wider() yet) @@ -43,7 +71,7 @@ load('trips.RData') ######################################## # join trips and weather -trips_with_weather <- inner_join(trips, weather, by="ymd") + # plot the number of trips as a function of the minimum temperature, where each point represents a day # you'll need to summarize the trips and join to the weather data to do this From 3476cb3a4e127a1ffb0dc067a2e58944db534efa Mon Sep 17 00:00:00 2001 From: akone42 Date: Mon, 2 Jun 2025 10:39:18 -0400 Subject: [PATCH 5/9] Day 4 plotting update --- week1/plot_trips.R | 48 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 41 insertions(+), 7 deletions(-) diff --git a/week1/plot_trips.R b/week1/plot_trips.R index 9829d088b..d893a3d4c 100644 --- a/week1/plot_trips.R +++ b/week1/plot_trips.R @@ -47,12 +47,11 @@ ggplot(aes(x = date))+ scale_x_date(date_labels = "%Y-%m-%d") # plot the total number of trips (on the y axis) by age (on the x axis) and gender (indicated with color) -trips %>% mutate(age = sys.Date() - birth_year) %>% -ggplot(aes(x = age, color = usertype, fill = usertype))+ - geom_histogram(bins = 30) + - scale_x_log10(labels=comma) + - ylab()+ - facet_grid(~usertype) +trips %>% mutate(age = as.numeric(format(ymd, '%Y')) - as.numeric(birth_year)) %>% + ggplot(aes(x = age), fill = 'gender')+ + geom_histogram(bins = 30) + + # plot the ratio of male to female trips (on the y axis) by age (on the x axis) # hint: use the pivot_wider() function to reshape things to make it easier to compute this ratio # (you can skip this and come back to it tomorrow if we haven't covered pivot_wider() yet) @@ -61,6 +60,8 @@ ggplot(aes(x = age, color = usertype, fill = usertype))+ # plot weather data ######################################## # plot the minimum temperature (on the y axis) over each day (on the x axis) +ggplot(weather, aes(x =date, y = tmin)) + + geom_point() # plot the minimum temperature and maximum temperature (on the y axis, with different colors) over each day (on the x axis) # hint: try using the pivot_longer() function for this to reshape things before plotting @@ -71,16 +72,49 @@ ggplot(aes(x = age, color = usertype, fill = usertype))+ ######################################## # join trips and weather - +trips_with_weather <- inner_join(trips, weather, by = "ymd") # plot the number of trips as a function of the minimum temperature, where each point represents a day # you'll need to summarize the trips and join to the weather data to do this +trips_with_weather %>% + group_by(ymd,tmin) %>% + summarise(count = n(), .groups = "drop") %>% + ggplot(aes(x = tmin, y = count)) + + geom_point() + + ylab("minimum temperature") + + xlab("number of trip") + # repeat this, splitting results by whether there was substantial precipitation or not # you'll need to decide what constitutes "substantial precipitation" and create a new T/F column to indicate this +trips_with_weather %>% + mutate(prcp_T_F = ifelse(prcp == 0, 'F', 'T')) %>% + group_by(ymd,tmin,prcp_T_F) %>% + summarise(count = n(),.groups = "drop") %>% + ungroup() %>% + ggplot(aes(x = tmin, y = count)) + + geom_point() + + ylab("minimum temperature") + + xlab("number of trip") + + facet_wrap(~prcp_T_F) + # add a smoothed fit on top of the previous plot, using geom_smooth + + +trips_with_weather %>% + mutate(prcp_T_F = ifelse(prcp > mean(prcp), 'T', 'F')) %>% + group_by(ymd,tmin,prcp_T_F) %>% + summarise(count = n(),.groups = "drop") %>% + ungroup() %>% + ggplot(aes(x = tmin, y = count)) + + geom_point() + + geom_smooth() + + ylab("minimum temperature") + + xlab("number of trip") + + facet_wrap(~prcp_T_F) + # compute the average number of trips and standard deviation in number of trips by hour of the day # hint: use the hour() function from the lubridate package From e9b21fb353e6d1e1e409cb15b123ec3a02d96833 Mon Sep 17 00:00:00 2001 From: akone42 Date: Mon, 2 Jun 2025 17:35:42 -0400 Subject: [PATCH 6/9] Day 5 with pivot --- week1/plot_trips.R | 45 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 39 insertions(+), 6 deletions(-) diff --git a/week1/plot_trips.R b/week1/plot_trips.R index d893a3d4c..0cd8935c3 100644 --- a/week1/plot_trips.R +++ b/week1/plot_trips.R @@ -56,6 +56,17 @@ trips %>% mutate(age = as.numeric(format(ymd, '%Y')) - as.numeric(birth_year)) % # hint: use the pivot_wider() function to reshape things to make it easier to compute this ratio # (you can skip this and come back to it tomorrow if we haven't covered pivot_wider() yet) +view(head(trips,100)) + +view(head(trips,100)) + +trips %>% mutate(age = year(ymd) - as.numeric(birth_year)) %>% + group_by(age, gender) %>% + summarise(count = n()) %>% + pivot_wider(names_from = gender, values_from = count) %>% + ggplot(aes(x = age, y = Male/Female)) + + geom_point() + + scale_y_log10() ######################################## # plot weather data ######################################## @@ -66,6 +77,10 @@ ggplot(weather, aes(x =date, y = tmin)) + # plot the minimum temperature and maximum temperature (on the y axis, with different colors) over each day (on the x axis) # hint: try using the pivot_longer() function for this to reshape things before plotting # (you can skip this and come back to it tomorrow if we haven't covered reshaping data yet) +weather %>% pivot_longer(names_to = "min_max", values_to = "temp", -c(date, prcp, snwd, snow, ymd)) %>% + ggplot(aes(x = date, y = temp, color= min_max)) + + geom_point() + ######################################## # plot trip and weather data @@ -76,27 +91,29 @@ trips_with_weather <- inner_join(trips, weather, by = "ymd") # plot the number of trips as a function of the minimum temperature, where each point represents a day # you'll need to summarize the trips and join to the weather data to do this +weather %>% head() +trips_with_weathe %>% head() trips_with_weather %>% group_by(ymd,tmin) %>% summarise(count = n(), .groups = "drop") %>% ggplot(aes(x = tmin, y = count)) + geom_point() + - ylab("minimum temperature") + - xlab("number of trip") + xlab("minimum temperature") + + ylab("number of trip") # repeat this, splitting results by whether there was substantial precipitation or not # you'll need to decide what constitutes "substantial precipitation" and create a new T/F column to indicate this trips_with_weather %>% - mutate(prcp_T_F = ifelse(prcp == 0, 'F', 'T')) %>% + mutate(prcp_T_F = ifelse(prcp > mean(prcp), 'T', 'F')) %>% group_by(ymd,tmin,prcp_T_F) %>% summarise(count = n(),.groups = "drop") %>% ungroup() %>% ggplot(aes(x = tmin, y = count)) + geom_point() + - ylab("minimum temperature") + - xlab("number of trip") + + xlab("minimum temperature") + + ylab("number of trip") + facet_wrap(~prcp_T_F) # add a smoothed fit on top of the previous plot, using geom_smooth @@ -117,8 +134,24 @@ trips_with_weather %>% # compute the average number of trips and standard deviation in number of trips by hour of the day # hint: use the hour() function from the lubridate package - +trips %>% head() + trips %>% group_by(hour = hour(starttime),ymd) %>% + summarise(n_trips = n(), .groups = "drop")%>% + group_by(hour) %>% summarise(avg = mean(n_trips),std = sd(n_trips)) # plot the above +trips %>% head() + trips %>% group_by( hour = hour(starttime),ymd) %>% + summarise(n_trips = n(), .groups = "drop")%>% + group_by(hour) %>% summarise(avg = mean(n_trips),std = sd(n_trips)) %>% + ggplot(aes(x= hour, y= avg))+ + geom_line(color = "red")+ + geom_ribbon(aes(ymin = (avg - std), ymax = (avg + std), alpha = 0.25)) + # repeat this, but now split the results by day of the week (Monday, Tuesday, ...) or weekday vs. weekend days # hint: use the wday() function from the lubridate package + trips %>% + mutate(day_week = wday(ymd, label = TRUE)) %>% + group_by( hour = hour(starttime), ymd,day_week) %>% + summarise(n_trips = n() , .groups = "drop")%>% + group_by(day_week,hour) %>% summarise(avg = mean(n_trips),std = sd(n_trips), .groups = "drop") From db3ce25bb9e65b7046072f23e6a5978cdab9da61 Mon Sep 17 00:00:00 2001 From: akone42 Date: Tue, 3 Jun 2025 18:20:39 -0400 Subject: [PATCH 7/9] last update Alou Kone --- week1/plot_trips.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/week1/plot_trips.R b/week1/plot_trips.R index 0cd8935c3..cb8fce9d5 100644 --- a/week1/plot_trips.R +++ b/week1/plot_trips.R @@ -154,4 +154,6 @@ trips %>% head() mutate(day_week = wday(ymd, label = TRUE)) %>% group_by( hour = hour(starttime), ymd,day_week) %>% summarise(n_trips = n() , .groups = "drop")%>% - group_by(day_week,hour) %>% summarise(avg = mean(n_trips),std = sd(n_trips), .groups = "drop") + group_by(day_week,hour) %>% summarise(avg = mean(n_trips),std = sd(n_trips), .groups = "drop") %>% + ggplot(aes(hour, avg)) + geom_line() + geom_ribbon(aes(ymin= avg - std, ymax = avg + std, alpha = 0.5)) + facet_wrap(~day_week) + \ No newline at end of file From b5a99bee9289251a5def12a11c8073ea9f0c1542 Mon Sep 17 00:00:00 2001 From: akone42 Date: Mon, 16 Jun 2025 16:30:13 -0400 Subject: [PATCH 8/9] best_model --- week4/best_model.RData | Bin 0 -> 52916 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 week4/best_model.RData diff --git a/week4/best_model.RData b/week4/best_model.RData new file mode 100644 index 0000000000000000000000000000000000000000..69fb97f552bfd0ed023bc3ef09fc139caff897a3 GIT binary patch literal 52916 zcmZ6xz|akY0@B?KFm(4Y+5bEk$M@dSx~t1vxQwh7(FLqu3SPuUmq&-~{|wb) zP>TH)8k#~FOGmf5Y7%#PN|CD0maL%lvX4YY7T9@aqc-ui_hRYnt{ZhtdwdZenBYF( zZFVEL;Ng5g`I>Zyco6H=tLRrgr-mpf9PlE{#G_u;7KD#~%w!Bbrdw2&ZtkYx^k>t!k_`nepYDA745(OdF( zG2LNp!O1ib8Wcn=^404gxO%GT_lMpa(bMP+3fb-V4X>?FoHv#?B~C*(--gjFu{B(M zdp@?S^EWx3o6hAMo=eX6<)~Qsh>%~lhy;=!O#q$(i6*Qllf=0shJtdclRjS!ym~pkpvICT-TQG4s8Z%s8iF@JN)kvW*Ww~M=^0g$v<## zuHWV2$J9r%8x(UQY97bpd2}kG!SP^5`v$EZ2>+z zsLd^u*XEfN7ZUz*;puJnM_8xeyjH{lk-kvwq(1_86buyzVoHCgKVThE(*M%O>5Y(p(q9L% z?!Aicn=E!|9y!0&6@X6<#iaV*Xf~#v?`%%7U46d3#nmXrijgGoUoWZ_j_sJ|C{~Bz zyEpC6ZC9De30ObSHbHL`Qx_~C%ia5Qb-AFQAhYhT=^+NU>E2gw4Baj*KOE_EG0~fm z#de=6sQ4xiJGQ&iD;Il=ndo*!XYVlRehg;Fk^xVbkGjN*_}yoI3aODzVfd(+7&BNz zLge@od^1LWOZanJJ{`2L=$|9qp0 z0Va@QaAk&6B`|wf@Vq>kwmUvJ3R4eDwfptBI@1@Sj}LA>o4_p^RbtEX+nP^R7tr#| zIy^L>bI@z16)rc4ZSPu0lv?~1)wi^cTg=@+_ovTz<9j(3;lpX%=RcolEXMUp;A?cG z&iEB=x!zUyF&AbaPT4zVTqnV+ob%v& zOhl}5-&C**Tqd8-kZ;7tl|VnLbk;8Sn(azvbGbNfM-G2>q3Sd?bzYb0+lu&%AukPSYTiet@0u%JUzMXmGWuO%Mnj3jrlE`)|ImEAf_kdN9@sE( zJ5S`i#~kH58;ra{&wV7YYyS#Ntn#8ysk~J{{#??ST?;;(GI}0*U)NngEfGNIo&bO@ z_|ntFL{nbW8wAD_&*ilH?C5c2xrKb%7SgBR1l4<AA# zU+=^vp?yS7uCF{AT+u=6F*Jp>*P1p zm1%c_tYVm=?Irlq+rI6xn;ifO{b;<68DjzEq!Az9oCd4N)YE?FQD*|z%8?GWn3$}& zX2kiKiabO_8I!1y=Gm{4mC9~$k|{Lys%fSUif<&Yt$YZh7>I(Uw-HQijavP^ma&T= z?Kf)Do{9EfYO;>m6@E276pX1zD;@^%eh7RBKuJP12Nq1{rH{lE{^l{CDUP~Va@Z#B zF(zZpv?$Yn=-8-&mN4wUD1`)o0(e2ZhqT>(b%*V|6yS9~yNOHnz8Jncp?vRmTE%1X zAypqKmCcjpy6PIg6>F$GHp_wOKOfGUb%?tZuG#1Xo!#Kme_0#){Qb>8rS`LqedyPB z(GSFJTN?UizlK;Wk0l#Aljw}W3dz?3n`7@zIi#p!-|e zzTJO^4=I8YIDa)S4*dq6_pHN<33MUp(cBGBUCGtm&P;`dTG@tM*5P=8QB zz!*0tsVroNTSj`9x=kB!$51CO-ngYZl$@Dm`$(QaF$GQV9ND3yQCuI8tC!X zs>d8#ZZWWZryr$|%xT*mF-gzn(!d>Kt3OFrD23ghM|M@WfgX{7;orBd+^ol{&F}C* zHW=QG%ebj{|E=53>ls@?ba!J=CpCP{^CLEyKEvi(4wdRKVe}Rj=(T;l)a=N)venUz zQF+k*Q9(v#pAv;jf#M_!5xmB*E+tA}-;(tD??L#zUFfbXcja{Ax9cfQ9&4cDrRgz! zK&l&sw{sdh4oo}}!$Y$Pf5}X)=+X=vGQ^Sab|sLTqZktQ($M{Fm2b$N6&Y%F6^)dt~}JH{?1a##%#B!q^uXi94!6wh1D6L@xfG zKGeS2{mnzs@B6CzJj!s{zxPQHb-Iu`f7Jp1vZdDjES)>?(J`})?N?n0&ExRqoYMU1 zdf;zBe;ezX>N2y;8U+rGWEM<&1Q7hFK6@?e-iaUtB1tOPj`2|Ge{M*Wo zc@EDz5doILw42!s`7K@^By~DfkGmWL;TYaLAwKGG#iG;CE5se2Lf_pnw|K5j%GZNi zGPxS9>9{m(_K%!DOg&2{E**?3BWO9JsMJc)pW5tlr<=*7Y}!NE-1$L82xa7Hs4bD$ zSRs|!(BD3U7wx66`fDEQWzH|Xvo}6?6wTt`h^uKdrUI8{vB~;3ZycyE>S&zhAIou- zR^2hm2Xp7rHY1P9&sMb4uO21Ji-dg3DK3lz`V#KgzDp5=6-{<{CV$~3oN~8vzZh^~ zs`S{S{N7A25O;Obt$9t-x;)h8A!n7hdlOU;}DK@ zd~{=(EW*G#q0P%aN|Ax#Ml~9sC=^w#?AeSlbCV66ivk0FOyHYAlGr)$fJ}SwvEf6{ zy|Sn_1nP5hQ0wf7zJlUn5tDr8NPdWWcDRbLU76-VK=jti#MBZ=e&50 zXn8cFy#_ij3Uw6C@(A!~E%;Sd!fp?UIJbb3t9OB^QIeE2$Rlw=0t>4vGqc85J+@&e z=T(t#a5&VGHnI{gNN?>(do4!Y=-5qrktnh;h_s+-#8Wob$U3iIII!>yi@+s@&0O-? z_xNy{!-+|(I{|Par`TZ6Z%ge4?e}(%Ys-oNhO;70Q@B}8?rl6tw85q+r-ZbjqVior z#Lr)?jC+ECCvC|kEey9p~E=*IO zlta!;?3peYl^?%TYDfDT7H^h}_)gvysZ5*z{}uK}+0tODg~+_fQ2ZyR=#b>G)I0l; zoZRkQcKDdIKeznI{&W61ZF9MxG4rb0X56~Yka>kA&n8FgF-M#cX`FcFDkD?bC&tXU z1CZt2u#F?fR2(%6$q9_P6>)Bj%48_`91?Cb`K$ z6NywxIvzJ=G9`c;W;e~HGD(NeRWxkDMLtEjyCg~?>y@5+cLf?7b5s1}uzxWAPB>>l zt?jE}Dr=Fs&GnjVx4@_Uv7IQ3zg#|J$l+Rpoe89^-s-)BL&zax$$cY+J|wY8xSBQa z6=MvOJg;O*-lEK)qOs6(68~o;g!|#y#SY6M4fuu-%STyU^}-JKqPtc>jIKE5->G%A zudbKEW%1rL%XLLjjp`oD-E8&L3iIB)qnTuKU$kp(Q2kQ*0Ci&L%s{T@-aQ$xfjKqI0!jHZjr ziIY+V%OTk=(?I6Jn3z&rlbd+}qWL_ct{Wt+I{zX8uQrx3Pf{k>_JDe-W81Vo_JyUA zPT8Fn=NKR1<*aR)%Y=&nAWLv7VvEN?N<~lHm;QOod%XICI+%i`hpvyK*Ly~J#W7cX zd;7v0yso%9zC`&{yA>65TJ$32;^%yC&hmeiYDo|Pvb2&yC;kO(F#z!61`7hC0#qhb zGJNh^8u9zQeN(3AW~rD8=&wnK?xRYXDqIygj~eh3(*qAbI(AHDfxsE#y9e|G6SZF5 z^9}MRv}MOQ1Rs({^_E)rFJ+ZNNS1?k|EdZ6P~mbF@qd)t`hMO|zrzxpYD|Dm3?_LFqO}X`IBwFHzD5Q+UJN8u$O(}Kz3VztjFFa*_s{-{~#P21|Q=3f` z*PPQ=RVjm3&~L;8e(#epC56fU{*{}!ROaeiBF1zzjJ(~#HS-7>#2)X*zbU`d~X>^^C?+_ zn^~ki0wQqY{Z9)HzDp7(a(ya1n3@V87jjZ6W&q}~TT>pYe)tnIpG|hL7qGEd>xW-g z@44I+*Z4AK!6_bbGMachjeW6J^B7?@zdNGJOhVtp;|bmyvMeu)$@NnKy``W$@iRAW zdp5Kv1H8N(U{ng4Owj%P)n+!I?>emKjjVu0R^EfZEgprN^GQl z+NHqCV@xDA-U)+6v%;}yqw&S{>4%(gVNWsrh{ESvH4pmWR|UMbV56=ixha8seCQRS z2Gfl%pkYi+udm#q@`ud{2E}Q&Fk9C@iVNOXX|(eO3k2c3VEme*a>X2^To%#mS6rAs z;jt?+ZXKfE^8+s#WHgLHy2bOqRvpey(bjCs!O^x|ZITo_$Y9cnLr`-Q@$D5vV(W5pv|ic|II)1e`}zp+2J; zkhYJg-g_{LgIw}jn;sG3htC>G&zK;JMl~O#1oC#Ev|n;HZn8kFp;)#zkVB3B{ zc+<35I_8RNWjC#C=dMCV&T^VA5)y}QKU9kX<-bne2!$r}?ZqVV65k!QhCi*VDGzOk zJ6dtbXR)eHqxaK3*^-9Uf5`1ucs(x|@)`Ji${wsSidgvGS(;cIcEnyK8@vlPfy*l7 z;H6B@W7GG_EbP)GY%^MFTVLRn+9f^r`nc#Zz+uIsH;3p-71r_u+1dO!6z5A{c-vyz zPe^`ua*=DN`;t6F?k5ND6%F|jGJu^AOmh&Yw9960D;G&N7_PRvoQ$NoxPzWyCx?vb z4+IXQ3nR`HF)%vZW1r&;-u?-#;HY0n-uQ;biO3`ztR=TqzKStZf-Qeelldh#nEl$- zyO5o)SfPnO`7#av@>0kb);h<)r*q27U}Kb^L!ODFbR3*`1&4Si+qW}zse1G+PZtiB zizSpKst_FPDY}|5`bT?lVXrv`Kai7e9R~5hABFh7HAIU>IvTZkcCVnIPGkUTXSF^gu4NALY&mrI9yXmGJZDNHaVpqORi z&qZfX$hP&VE(h&bpW=(P8Sx;}d3VI3!)0C?fjxj)S-_hQjeP#CR`hmvbZ7T&>#zt* zj8!FZK33_m{J$KXVh7(K>BuG*Q|GR#(8fZ|{=kpijFx&cA3rY;CQn%fc$cOgO+=U> z(*5_t<6F?oiybXhqA)={S6~>xGR~m1e!jUJTW=LQIfY|2-AC zN|Vz3``_~a=yYX0!vtQmWqU8FgJ2u20>QH>oP0tbHb{v+M0y>O_>!-8Gq@gpj=zhz zN^)_vpV#(Pr`zE>8xX|Mpl`aGK-jzF%TqMzR3F56CG&Q8beGYdG0NEs7%q>fC`DEB z-0YWVD*>}fH78cci}hb${X`eWFKWJao@P3J6^gK z?-F7c&s|7k6yF0+h6)wiE4S8?%7z(x|EOW?B)(0eS^Y3S_}*h%cepq1`FL5VTIOJAD8LlE>DD7a71&u#9V8|U+6rkMSpLbsc2~dc|#J^$e(paEMuIVTxzUe-Z zGt^tF%TWSH3!Eg(LOa$_GCSH)cxk`sB*_((e--0tSsI2C$&NFCBBmVTrlgCL<#rZ3 ze-TdmTc3P-Op@f|8r4}kY97y6`b4F#xFg1$#C^3)HE6JI-Sy{!Rtcr{kpFtS6P4B5 z{=7K5of?;Kh7pv#eMf6bEN}XpaP{IKTZuAm7>LqYTD9{!m3xl6&Z-U@P@bMasQ;e0#~~DB1wnV|X#g{5#_G zHN9o7Nf6OyiPHCO!kb}q7ptpt*fsadv*Zk?GS_k_j$AyBMcA`F+v1SC0z{cLE#F7S z{aypPgid6qtHx~w)k{J8^L6oPIkEh*yVRe@@)#5YSw9A>lYAN~nb&&MtL~EbjV|d(OV;}_kBJABcryi>Y69>Ajj5J zBa(-x*hr>82JR=wN9H;9o>v=(9ELoDz7{?QNk`f+%1cRTNMrbNHI3X8G4#HS`glCr zm#jCb_P${-THrUCuU9e?oh5g-fHnypouyXViN?4$JsRO73T{L{eN$DlHkwUUVS>~# z2Hl_^c)5Vf7JP91x0ovm!AEbV2Py8}6-MCOCnN3E(nJ%P@NUQ?X5Y#_*^kZICt>&9j^nLGS>-2I{oRu#&=0M0sLph39LXY2~*}j>6S75Q)>k<~9BmQ1C z@jYYFE(HL7P{UAdakxMO>r#yD-O=OTFix&J8yE;M?n|1=dQ_&!I~~k!2&(FBb5yFo zr4$_e95+kj2LZlMcwz4XMxNLB^Y96{I4N(C#CG0!q#u>^zDUtL=kK_X9%f?~_CD`PGPFdXVK~{$BqnlH@X$Tw&AZ7rq z7~WY>ZzuD@kSngk1IhQa9BuiLAq@S@9ah$vLCM7me_Ux550Zo(6JL6Hy3)xF++sUm zaWSG{toq)yX)z)6LBcBS8p$zE2sxrYiT6(?&fAL5y^|?_f!vIr3Kw zMd@bUzKZ#Z|B~pxuHsw_MpyZaV#U2Jc4gg~Eb;dZ;{*^BIwJyFcJw9K!#Kfcr^~0u zs4eaglyecZ_GIx&v2~1g>K}9*j5t}oQtZ4yZLYvU6IxgDtb03TZ0lAp>mNZA)LI?a z{|{YEw6DgV4sQmxxiOH6YV}dO=g8yDpX%W$=7PdmpstQfGj3r$X2)7oAHpHkwVzNB z%1glp&WYfKQhPrlJCqNtQSfPyT8KRR)5RfA*{e)6+Y~Vat1bO$;EQahv&q1T6ISZc zL~eDShp4>qBM+%Kcd64q>vmb1{N34as+TIjbQk7t1a{$|8-oa98MTP-{oK!FQ@FPr zS=$-!<%cnFqhfnKdBZT@AYGgL;2FuM*X4|1{+NOxMYh^VUL<5Awcj1UwXXLMdXHi0 zV}G4Py_mS_2uEBfaI|a!DQ-YV3c0{BfYy}UIAc6sFR|;MP44G=-Ai|D)OldXPEf-3 zZV?&u!Vhrn2z*QvJOte#w52w;Fp8zY8$*CQr*Yv5=-lUoAXm^l?jT8{ zTj2Sa1mt>4RQV<*{UY zFRoMQmoblUk0wT`5T7NNW$*+;CNW28>Kn{tevV@+SaYAya!JC4^PE0kAes|~^ot6v z!-jAdpf%RVkd!^Z3AZV#`AF4!>W02IYTC%_EdOwE*M9(_J*Q|m-#rggn3+-EF8a> zH_8+LXv}e@vhe^3c52d;2yiZ2xHOxPe1O}7b+B&$l6aL7rtVG2mfp`MwjJx!AMR1^ zrniAFcRd6OpBQS;`NZ)RV+w#cA9VB%VlYc?HQV(7L zI=(P*xu?Wkzk~EBUh{u5+H9v;<#kWGB4)j4g3lBMohM#Exzslx86uup5a1EzVCyh& zJwgaYix9}Qu2g@G(AYwHua?6=cV+$J{0FromaUs7-AXlRGfjI&w^6{FWcThtsjq|+ zz7lIaO>10Le7X;=cL>X;feUb8JxiUh={e4ZjdcVRcr{j%ST##AMJwHE-vDw+{~I}4 z6UeaFR}YBTWKYEc_Hp~u#1(w#V`bvF6l}&jD0rV+fV@YooV}6Y0Ya{53rY$58j0>v z+9oeV#~93doe;tc-xt(Ac^xRf$lfYBxo?zbxqv0N!={HG%*s*AW#4ZhdA%^JR9v)y z4bQG^#MB~-9;@An9);hyWreP?Gf;WX6+JOsizEj{5EdTt|+Ur%;JjKKUWy?C_9 ztJO{JT?7B&1w1z3<#hn*ulG8LsQiXMBHuOq<6Lxew|+uI)ch_ya0`hYjm&3SOmFU+(hHL0r(?NSCR{OWSwa zhAWFNCf0|*s{;r-Mb*8C3N;iwN3BlG2@Kqz%?C>VbCW|=;5KV`|MnkhsJzgC{XI=H z!Ser7J73>VQ#NX!RX*0zwQ2eqJHZQ@s!avuTC)m!f0-Xf=>wB z%TU5p&&DfnraZ6L$J~EehnR`=Mzcz5gN2%{j1gsiE|i!i&by*|CKn6cy0PrVtTk+w z6K1C&-hl@%H&6nnW((iNP+7qK@QVGuufcxCyn**yD``%31KCGC*BdvT>&A`Ld1ro| zG9Gr1u>Nih;K*M7B|+TL@G?9Z;lgYK&Ec-xvGaoVv-1Hi)|3W`{Wh}o^9`~Nr3OPrf&RFt2=ynn%%U=GsO;{42`dwN z_{DYgcnvpYUQawQg}0plC@deKrQt+QlTB)&r=hbIHwEVzost5Rk*B`DUXWOC>J@!x zJ{;ZkW~k43>yEMbC=t*U84}_iTIz_N&%Ym9K3=VT%o_~p_Vw^8%b(m840|1SwPKaL zSB}4GS73FGyWJK5exVvBV6mjV3F0~oUweGpwde_=mfXln_WVwVew)Bb-3@#g-6cL; z*Mdg1Nm3iPMt;1Hu^@=Ztr1c6X>O-Rh)%;$OHdRCt9p4bn;QF6WB!x zcAxW*Go7`9^dc<#<)w=xfoUoSE-<a}rDM#E)) z?@IO!PxDF5`PPmRu=lymab`31_3$bcsnVwhe0>i1-`nVk7lW{EVEmqsyNA6xb@l8C z!Vx;m(aGGA##p;o!BTNryd8X>|L^grU31S2#*I>zQx@KjA0Ten!~ERZ&BeK7jM=O) z3`k*A@WuLX-|;MXiv0<&JAK-B+ENzoCsifo77C2K54%CC(P1G69jt#m;^DeMB}2(D}Ed4SZQ zO0 zw1~I-NgWfz@r6)%FobttPy-^<(Wr` zh^+Y6&nKXLTD$mY>C9n6b1ptyW1*9~GQCKv5YY1cV-n-nujrZY{`hzqoCVHb5E)V}(2G1_j2y+x%b~tK>Kb!AkNdH8 zW~S&=r?~bb1DJCX3$=tt*4~C4)l)<6<4zaY2U%(V&O1-!wmNvp7Y2+--S?|2V5Y90 zxS=KU@^V|uC-jyV*|0?DC%D8=j$EfnE%xJLWF-LugiUZ#d>#1_2+=#&tuQmRrd3_D z;_b(3R~`dR3wL%?qJ}a68SCf5)O+8puM|(fbQ{=@j0wxz{ltAJR+9Y zjV=9rksSX>T^l%})~=i?(12i0B)_+KFb?aipL@m3Fm|ta%J;~k;Jt@xcj3#WpO_M~ zGN)gp>a#xGgEwt63`}z_5E}VHdz2rpcJVPW32@dL6uM__`APQj$MtnML(-KKVgi`16xr8mc6ugo#HneXn0%bth}PosBwJDlhGBHQ!kK+F8TezT!8u>(fN36@WI|2qSl-)m?iqLbM4Tk9 zJB#%lAVnUOA5fCjSL54~i&Kbf8E$5im_+ZDvm9*?$4J~R-n7M+=y`R?IDB+6C+OpJ zqp>-^hIC48(m6fs)EU+TqvX;(mawRLjt0<*Sa~#PW1dj+>^ceBg1d@KDQ4EZ7bp7U z`osk7E_uH3AhZ(_7pK+iBWOIg9wtr3c5aJ;s9IyoEw=WF`<~x2yv*g&s?L>mw{B0t zS8ZbMd-csu;p+vp=`2oUT*+QKO@CR>H^%glLk{tlIpMgCC>*>?ZS$BCFGv~tF95<; zPnol!69bJBZ<>wREK-@}MbITTzsnOJOsrZzw$0*W2>QlgUy>Ez*;>zNkV0jYqK2QC zHmb0!NTFM)#YCDbh(_s;iI9wXydPy4y0*LVKs7qOZ|B0=^K#~a z#Z8y|wuH?CpyUKhEKsp=B%DJL&j)E2QPBs77 zc-0!Nm+0l?8bb-MY+&Lxe@%w+*!OcVfBqY~taH`yglnzq3z*3hn=ZLO5Ju3@-Kd>A z&-S*TlHsd84X#>;(C?T>3XF&(+~XqCf^5&IWiLw`Eknp^Zaqys z`p(`M*Ze0PK$9Okh-3kavWkJU#1%Slel1u&Q4UG&U_Moly7ljY5zM2r&39ztqtb3f za*HGbPvb4G$L@xgcJYa6>Yg1UNzWj`LjN`qGk342OYtcgKS3>g%koIyP;IUmV@ArI zrh0M&l`%ELJKFFROzKJNJ2>mQu!p8%O&t}SpjsA{T7UM!#%6ZA1-;_NBaUtFd!ybU z&i$Nj!*u~$TE}=`NQ#CMrHLY&Ciy&S3@Nbl5iW)EQqo8_j6r`~ykI$%mhCqCnwDJG|Ey(i|FI3bQ0wl$2v?K6N1gq9XejmL z>u{r8EB@|VvGdkq(^ehe(g&Y7K*z_F_f2=DfZ)g$W2+s;y5P=#)cj#+j?!C=f2xbn z+zrGfP%?9UtNE2A8|eUyms*G4z8-EV zdDT67Gw&`4!Hmd)P0@q6F8c{-L3`PwkP==p-_q7~MeD_Mm`0S8gr>LQngz5|>T`32 z{#N&E3btIiD#|icfRu`Wxe&_kee(kKd_Ln_*xi%${Kiw=B_OhC;k>Rj!fI3NjaZ~H z1D(~&nDsU6pB#x?aFJM-|Cv7-C&v`V`U$~uTFPck(u6eMy2Kg9mHzEQXW%pUlP&j5 zyPrvKoiMkzwz4O!#by%8!iN(HaIhs@PaPnihE@9hefYa7hF6>1Ud!P~ROAlnc4D*IDZoCf2~AnL`)1wXe2l-JCe6(rP7$&x2}V zB11H~%G1d&-h86ptRWt@jNYhTV4KP0DQITqlFI+?-*aowVkx2|TZ9DN>Jk}KR2LEU zKjlV2WPJxbR$gebX(Frel&beX@$gYF1~LTy^^IG7_#EG2KDUxCb}r(zwsONc#o=xQ z&VDd&*_-x<=ejQ<2@=bKhPa#(%u_rrC5vP;Z! zT=M`d4nQ5ktniOeyL;L0kesLue|456QT8N@`boE&1^-W??}9ke?a9c3JcVCJ<@9J3 z=o{+4_ksKu0D^$DvhG=K{yUtAl^{G?pm^t)!|5!Fow66Pb)x3|@rPNtO5NSmrnpYZ zO)tqz%bVa(wgFPbtAH1w^P3%#|47ZAJEujkb^paS@|$+5yKmkM*>3iCI;Z@H_NAW* z!`7V8Lz(U8&33wk@EE)q_I&|-a!VyuH}dmahF#q0C%{^J%?cFYe=Uf6UIw)_C_9{K1lU3U3b_|HEgqtgsbW>E+NOZmuYHpIG-<+? z&B1390e3@2K8&^{IWVM7D8Nc z;U}4VQlP0eDBC^hBS2&VnE?Y6n{dSl$y7gCAI8|6cYYRIzabMzc~=)$H~8hZBxjN} z4ixLshf;^X7^C+d?+{BfWkEVk1BN9Nm0>gv`Fv%ozM$n)=W6O)rZVcBG7oSaI?D?n z$wMx{WF(+3Y4<;W355Or_tIGI3$Uxj-&@^_wsxkiq5qVdd7c+x!`%!mCvE+gjQKHV z3et_6IZgo(bZ}Zd&1rh{!l>CEOYi#5{aa@YvEtMEXIX2CkJJjwAaZ|d5%CsExX^9) zc>%g3UQ3pe+2Ng!L88ImAVy7i+v5HQRBE$OMh@(}H-h$TrD3I6vJdbCuHb3|m-CK+ z;-=S^d6ekXE)dCq!-0jiRZ*H0!fDE7$)58`dE)TA5FCAkC_e|uvWenu%j{P$AV^`a z>QxuTpfffn=o>eYX#+RYC1K9!Le2Ax0Kl`_cr zjTFoT5`9}*Vk%UltNy|`r))od$JwKu2@fNRe}jgNPWQf@`D){KVl#AWXHqr`IO*d# zNsY8YQjbKHnla-0FunqJ)tGGU+IPix=vfCTZ28cpgo5=!;M~RQE=vc8l2f<(XCkjL z<8%_ueoOIn`ChGBv5CDOADj61JG48dC-_n=Fe-#&EKr&J!=T2B7EQpPkwU|sTv<^z z(`pKTSkXF~jJKdCjwLGrtjIz3X2Ee_4wrtFB5GGLsQW7^HBA(lfc?c*#8l59minZ0 z&<(nM`*Wws*zKvy8T8K}Y*F3#)jH2xU1&>*OU#Px6UkLt*M;*%hG16JVINeNG2b2C zC0TTHhXs7Rl@xGM!SWC3MuBdGKog$eT6W4CmpaTrt^E1K^!{r~c+F7{9A*zHZHAXI zNL@$_@t{A|-~(JE0za5b>aBrMCpC2z|Au{9WxzOKm_+R9C-~`wVh7x@xlgb<|H~_o zV|e1{F7}XrZ^Al9mQ!jRWl%y~w|;s8q%X$N?NCjw>LI0lqM6t~CFj{}Kd|d>Koc^0 zK%^D3b(3?vohS8Tx51LMLrdiA zat$AVUYl6L0HqZ2@80Mhpplu2=H`FgC_|E--8=)KkD$mo>!b=S|F$H&ceJkRzjf0f zxNkK9HxfG#7xG29b6v{K0C0Nwg+>Pq&+y8)=96qv_nA%PG8n9Sq*vD)X3S?cN+ zOg0>FY;s#mEq1NbDgwiiX1TdhJe}%Z{+o*M3_lp$DsjI1w&$DGU#H?2i#@r4=;qU} z+~o9k$vSf0qXYpV7oyQiOFtjJPLr%Yl(pCUuU+^&q~{%3SKqYEg0ak;`nL#0lgd4L z4o^sc&1l^t0bH1MT|nQs*uVUwPG&7uMKno%Bx3&HMU5t&`hMUEYXa)>!~);b=JqP4 zcdxm+%x2;7$I`?0q||Qe=Q1s78zB7|$IX`2^G1F5TZOI>f<6aR@%y{Gu9G$K^Sc!# zjO=_uo*=29zG3SUL-PFoMHce|8&oLU@Ru9GTqyiLeQJVlr|LL7aKGnSm+vj7-;ODjC zj+W^Jia0M8ssDNUz5Jrq9yGKRx>_5ePH$)vO}xlEYu1{per~&S$o>WKqh+q!f5-8z zTuA}PGFQjl$AJPU={jxS_*O0tHIR3Y$UTb~1P?5A+GL6=clj<_2m?``wAIcM|C!Sw zwPs88`*6rZ$DR49FV!wCZR4#a-+iv}{ARR;Bz6PdRJUTRquulaWtJGJ~&x19ypI@7HG>g z5EG6+a_<2ooRJ^q6W}F#ztf>*EqTCQdu{0Y(%%O@J7SxG^WLqx4N!BGhPT$mNqSJ# z`|Ks!$|TE2=rNGrB^&03Kb&HZRTDOEXvR6HP;Z`_B3hW#W}3(-+~nLh4)@1rCvD_~ zHft=l|4Lm9nuWtN5q}HAPx3MtfF7-M7SvaoUae2O(P~ zMtP8ilqbVYzU>um4hP>&NFW7!8x-Uvsu~$@!mtg^1j{i5B2lyu*!@(1)*& zW6_wrmw>0%$xthU75I&PXBLgGu$Pq?2hKig9y&HnFIT|)T-{Vm(LLixN`+v*%p)&` z@ICR~w>2L+jNadepnxB8Nz_hD-0DhDx#-i;L5VV0gsRRL=M~}pNzc25>kzAxmlKfA z`OYG~FBYLoYaHlDXC;*H)%B^*o>QV^mdEEjoCN8a7v`IOLZ9Lh1#3y2mQilx_`#L2 zIMW7ny!Y6zqhcT1X>skFgEMKO=c?J3({-wkUy}{1pG}AinNNL(?WL^7F%{mJd`!g{;)RDg_FBAMK1&8>n4)7$sT9Qgo??@1f5) zCDM#GvxIf9%&zRBu!_!+5>UEV=5Oc4eRprbao8l>r(A$+t-wLE2a;ot2Wd|`d>?K_ z?G>@hY!sN9JKk~dUwpjz|8*v@&&N$euoIsB{ep_ZKtLo1vaKTYUx~J%nb5*lNzjA$ z1+rU7X8nQv%lg(Nw#3?Gfq~vs8^A%;SmKlt7P`miMPm4QDccaPk4Sli!*_U-_2sE% zP#F8LJ0xHhMcyK|6^vdPvEoWbVBz#wS+k&ulW`RN67eJQ_Yj&fZs@*Bh9-NSmclT zK%8(ZF$oRj7Qe7FRQmrUr2(UUyo}PvlP?LG>o@(8gk&t!SIfLHuPU3oGLUXRqukvC ztHVkHnvaN@vmAD$yx1jT;KFo(w4#4^>oy@SL)@DrV{Ka_fzADpLmYm{R%+N!;77na zF+ttn>-9a%_6DCUC{@t+|ENzI1ZvA)g9wLV%qGM)wozcfYs+|2!ONpWKU?>N`LQe8 zPJ%YzLMSY9`nd}Zp5=y7zCoH!dB%GhO6 zEzHn&AVLcJTD;v`e=`oYcKA1w0v?L8*uQjkp4 zJ3P-qDBBO4jV>PAGtwU;#ms@-@~=f}WNa6GdP%+EwP2a&Nw6a|kRH z=xy((89!T?%O{LEnF%}F;Sja|TFZ__Np7<}0qqBv}@S2y{##O~? zR!v6=Ct;5OkeB?|If^pc%x)>UKmVD4NT!tz-i=et13qU8qgoCIeMft#9aytx*$B^t zc6{;SY%{DMyS);{aZ~yC>Fh6g^`A@h=a;`m2QvzuQk9lz#j}H=1Lu6#*+h^h z#B-_{YpQK0Dzxm6J=)xW{8Oj5r!nCXpH0c{v_Tt^?|Zv7@C_0eIzpZr?#X#}?uW*x z)%DCF@0!#7qIGNUOL!uk=vbWz+-}|MqIcsZwBzRR%boN10(8(m`a>W3EgdNvxu4{@ zJMtjt`BH8Ehy&=Y9Fat!9{W3p2r_s}4zR@wtv-qJtUdhwjQc(x98 zurIb4WwYvsn|t=EQknJvWq265ebs8`3daoy2zoDNixO>2?>=e2-y9Rfbky^ zx?StnZSovN$bH^T1T*KoFnw>r&ZfyC1bfLOauYGV&FebC(z!N2R=e#U;){-CcDwKzK2UO{7|9w<|Z_a zG9#r5Yi*)V;y?!pbhF^~%pyiu-Y%1Tnci3onv>(_DLwfHM}*kp8_4Dib6CCNj4R?u zbXa{k;d8B#Q=0Y7{)PQP3b*&KmW@qs`vD$vRSN%o?{EA{YQ*KA>2pI{OV%$3vJwmU zVa4Bg<=&6#|0IN)_=Vlw1Uy!2_`auxMWwE0IPIjr48vOicx9(!pV8h1r%RSrnTD-S z5NmZ&Q@CTDVSXUcGBQ$?iY)2YsQG*>fmbQ;VmsFPQO0ispE5}2r9>Fo*%h0EH0BR> z!WhvPo&@W(>CHHWx@tgRR@-uTm$zJ-%YMC6E_j>Q8-882_ z>AQ>egoh5p47qd3BCLrZJ-YY(tO^~>&TbmEwmk%@s=pgj7{tM%r|Gp$Rz0A!A+L^- zOb6I?hDOdVeg#hS9zuk#Wrp{GVmAp2 zh?jLnDkZVL=dBgLANM^As+&H&?LJQj6FJiRQvBF3yzo4dtEnRkU15gytH%Ug((ePm z3u@QpHCI8rY|Sdk_rqBJ8}Z>?9NU=YPy)Tpneyn3?xMmB~^ zp3RXuT46_!|Lp{^M%*y;r&%%N=D!8ii?3-KM9G0vPpMB|2#2uzQaG8j?>=D%c|lhn zj5rwsyBU&wv}zNWU*YrYDNjuDv3!yiO`b$tv4EkEPfv26{{;gwG?)q1*7NAeawe-B z#`6%*7d+k=Kdi>Xn#O`!u{p%^)pb8E2mK}ta^7F~G<)e2sN$kk;ohMRP?-6sm*^+- zV@DpGfA$cc7s$+2dN*FI0*haK_Rw(vz4h8*d%In|ngpAV3C zx8KG0tHr~#9$ce(2& zO2C3saSuPYssfzAb%MQNgD{@R*jss-3>G|KY4keo{|Eoi55GJ4@WHT?JM!s!ufd!h z>@UnOf5P}*rMLSFb+KKt|rQ!_uT4#u#A81E^Xpf5y4X7ksT96ce0sRij6YmzEgtivz zL&=&3O4oXKeon3dIT}NbH$`LNL!l34r7V)Lhl#clpUVVgjIvJWL$Wa5B=C9M&^7pl=O2vv$kS-EvmK(j&Z_AC3=Ai{UX;YI zhLuT(!Q$694hyMFWAl&K&ob>8x!1bKVd(|xlPm|1*abw`G_+Fg*%Lc2p!6_y)4Ca`*@ByorBrU_6UELw4rts1fw@?xe~S-2nqs zc7@RWc!lZH%{^s2*H1G4(bKHewer#!j2}1(YE3!d0c{;?cSQ1dVRk(KB4W!AOdnDK zIj2-lKeF=2^eXX@@DhQ9B8E$_$`qF%kpPTCqMnbA-2gSW{Z4iYr-B0GOx=4EY%q=N zT82%O7}RayQt>=ah3RYF*!;f4L#z1uVbMXuCx-|wLdWFu51p8XQKzRDC$rk>(NhM? zX7`&!P)?u0P+e02^n!lM$S1)j^x`|YSNRKUh-P(*i2x@xdY7j|di+^FA|bKHyFj6Z z@=D{bo{We=`2?oM%YUCoxkBarjq7v~!RE^!Gu61z<6pi=XR*2?Pe=?Uv<7t%5%WAw z53*~(c1T6{Y4s5hzP6%jeT5(7vS`mgdu$HUbWY@O*Q9_fiW5qnFC3uxDTjvYU9BML zi(I^@zX8Y@%=>xyv=~I|J6nK>Gzj9^aZFK^AEY$xd|}ml8fI33q);G>_YZ?tdTY5D zit+wHg4bt8`k9j;dqKE{IPZ`D>{i+n(Wl7p_hJ0+h@U3?K3JaV=LQ*EvdJL)-NjiT z=?>Edb&~EI3xU`y!PW@FA0S&#cBDc`4CFj*Ni_d<1y(&VqIi14A7rPjlIqo*hXtp7 zm%8k(u5I~U;u=}ezu$8Mdzwk^HrqnU7ii~MOrJKJGgf??D%>gFdMi+;|LEL! zV;4*l_cXPuy#$MXEQp@F7!I8-CWWMNhk_j9kPG)(BtZ6!{_2UHCeVxh?P%ohD12Nz zF8H`GdR&l?U2+xQ2mBthXS4b;){jtvxOkBzvt6#3-FsihNte5yK}IwYn+UlNEN;Z{ z%TL!~_So9M2Au6R%wI5vL627CPY@c)Z2FC|4(ogLQFHf4cAsGR2|B|Rg;!8Mm6P)@ zV>BKwj7OQ|`A4q@KLVv;x89gIdx9LhW4D&>`hhHGCI2Q(>RbT zn@=os!eSuuN5HT5dV68+ zv+q`%7TFm7{^1NZ@;Y50(LQb2QOp-QHYp2ySLwj`Ul;H4lKZL=h&}V3wXZ@M<56G` zZn&bV8P9jv8LRQ4d7ui~C&e()wi|+EpZ5Hk&xydE)#=e}q9*3AOyV1y?C-?S56E3z zt66%y2JBA0$#sps3({m5_=;8wK*<~AU3KXm;K0l{&#E=H{`o5|C1jcir1#`sHB(@Q zx%GR`lC2PfJfgne(Qm1?=mCDtqFCB=S8f~tDS^?e=PNG(U#)=3$YWI)Z+veK%2Q>|Xc@H!$&-#u zeydlp_XeIv0g}20Wi1Jzht)F;w2BQRygkSix#)}OQ$EYN<5e>+Fn?v0G{qgRYgj%t z+V-pS0rx=ebIZfmn{_eX;Q3d0)Z4c1`)&{waE6q_auDNfX2NT-1#vP=7xVRjFXc%A z5Z;^|Rmj^8tIq2g`Eqms-0vXg82WS2WM0^zhoHuYIMMz9twX>+ zuk58|4+Rid7_2b$vmebj?c**@pGDRCgtUkEEu$JuIb9vApHYoe?n}WbyU~ZeUbyn* z7f7%Z!*oa=0qXi{`iFuXFJk*<$DR^SPSohri^@iEJ=AcrR`qJFDEjbrhso%BdgN)w z?D&3WM)Zy=v2uNDB$9lg`0AoaC*rPBaCM|&5m8FeKi$2ra9y9|zV9-i0@1r2UJfc0 zfb{dZbzyB9APJ|_zGLDY2)(Fw?8Z$4kP*zLwtS~~9p}1q_Pv^*+V#}F%PIsQC4-)u z_tim=K|Ql@L}dMYa!T#KCObP2?;F*{a!v=Nvp&p>+h>LO%TX{WbaW60;pd!2zWz9g z`Hyj-{XJ(SkG~It0xw1rMYUso^2`rK_SwqA#*}2VX#XV4a+{IFfe#@6>QJ#((Nz#9 zcxtzl)j61J^z}K*m@};AD60w$O8{PKj!c&rUcjQstEuef77C-0MKt>tK0hnR^P$@sd9Tv_u*%Yj- z)698wFU;;gd9L3_0C=T@zsRLsfA31DDvS^x1c8Tnix+CfVRHO5VK3t?SRtvZ4;`+; ztZ$JXO@~f`L@!#ZQI9!La5d1SlWISRc02lAHm(heyMYT9*Ot&}(SRO@?*o31tx7O= z=_wVR~0*N^hAbu zkOzTj!`$OqaphRNl$`!Mm!|??6pi|*2~HU15$DZJcrxMf!g!Qo{Q7=1dm2m_UprF1 zmIh+0nLDWFz5c;wRFcdwidN3V1a-Qu_xzf>{p*6?Bg<;Q0h}Dicla zI@(}(IFE8z9}*_`!{S5y9W~<)4$?4f_UidNEYvV-V*0~q z+ziIU4ZdJ_Z}0=8#EsHRylBAYaD%U)z}MUN6-hmeFW=$Sp_L9ZJtYX7yk5ejvK^67 zu>?jE^?!NO=!VY=Mk3RrV~t+0`XEVz!a)fbD8p8kh|>lMK9<$B1MDETPr7TetrnY) z*r8A42ae&e`S3j-e_w{<4ovuP&RIrF7UoP!1xz2i0dq(uMxVYe1C3K!w5Q_2G2R3l zcdoI`e+Svg47h>st}veDc7(v_LC|$=@J2O(HvT@Q0|ooOP3?cO1Y<&@SkRc;*mJy& za1{9DIc7BiTpqdE%^CR3&PuJfI)3mRC)n3H#x8mf#-EJxkvYK%2i*BXB z)9Wf3$EY2i|1i&d`hNF0GLW7~(Ut#{3eO{;rx2~a%XvSZ@1X9((tO~DW|&juqEf(m z0c6q!bjsL>z$_WooadJlF@NcmmkhkUJn-`a@>-oJ&D*R2jg5iYM zk<2cqnvZ^GAsV4gs>Jpa(@VU521~Kdmww`bwd60wi-jn$`N8uZuP69;8~wm|o!oPt zPVEsfR>y{}K#c>$+$LB5)P+04BJ6l{5XMl783ldx#@-uvp87xd?LhhFH1Ky=__em| z2r>%pU79ox!upkyAzzC!C}a6lTO`kcwJ|K8%w5y_UP?U#`9FSsMTq+`-r)K7@-*u_ zqzU~$`kvHdTZ9Q2$7C;Q zN`q+HA1?DWNuc0tS(p~nB&M%XGAH}H3EyD#=N}* zgEIB{31pYh;NzD2!ApHKV21xr;W2qMq}S-|z$72i$kI)Gj_4a&DqsB7??fJw>#N^# z5wS$$pMN#-kobYdE-kyu&wfXPyASSUo~=R>o=-*$y%s=U&{#4)e)Ap~knl=LA52xKzPqv{4>Gd)+7t_i zp%RTAB%?2c>Ed3g6CTd+S=Fu!H3SnN@orBing1{7l776M_J<`fmb$LMjdXyt%2a?G z6$5!=EmH!xOw3=(KAN-}43}ZSe#_UoVm8p{_*U06hO}UW6xr50?cW&vvisq2Z3Mx$R7UkhgB=Kw0e&9ut@vy z71G{3nAq{F-(*+}7OFOu7*w0X(7|+vlmaJ^wNhdg$h`i&s>c09;gBS>InlneP2wsR z|C#7&WJ35ykozS~j(6@Z$f`6IZyibjk?gdFikJF9K&_kSC1NX(L6kipsWS zWL}ZS?69~SxL|RUQWz%om|ev80l&xW>3mFIuAC6cKp0A$#2)!1woES8KN8nT?n^2^}fT|nyh2&VggWgD2N23Zd^`Nc+i@O~lAr+-(` z!xel!Fz3uymcmpNq?7&5P<{Cl28pq>mD|6CMu)!Y-;xo>^AOJ$4DStoz|_FgZ`p$H zVsp5`R}iH>Fd(IH0p_*s6SZ->0SrK@ z_)v>Pfei2^=uX;SuMSe3c|YTdKfz{c;+sa?`G0uR%UV6m^b>^2b)8f^ECO+SBrKj? zd9YJg^P0tyI{rSU1EFl6=1}|7uLWu4TVL63?x1WLSHj}15*-Vyg%?* zikQf}j0&3{Rf8ePv$@0oq*C+;97+R^yNFXSQES4=udAL8%uyxhj|-oW!TUxqrk{&)`5IY>&tF3JFG11Y{mDbi#83NE&LOY!q9hEw_4 z1yO|yV^}^Z`s297J3+v!$fwgj`rselZ1As(=rCjL7$Yp~8K(4X#9_QG?%}oBSAG|d z3w*KS@U?YUD@<$>=;|Tr1f{ofjrH5!!jPZe)tO!fgMut(MtOo$AYh{C$>&Gfn7&34 zBg7r?ky!naJ4Bg{ljuNRnDj&uNnrK9DNPc6nnP$4O(eG*9V^;&^~_#bcpt6f`k4r# zme9ufr@y$|^^wIF9{U^IPN1{V=9CJC#K`!3dfQ1iTJ-g`dqHO`^wF*yD~&s2sc0i* zUj81R{YVdc!nwc_Gc+xPyKa|%F#m4YmwUIu6IphjdMU0jh>Rz#{d(g>hZK>tTHBRL zqvrEss5@6K$Xpo|wfjvBa_+I3q#%YUcs( z-f?IRtAB~6x`A}DSB9tQIbp^pt|(PI8yNdj=JnHWwlG~S>INXB#QOtYSF{Jk@A83E zQnrXo=rPQH1{Eb8P>9Cv=Tp48c<-#E z4zGU>$^lmn3X`q#>vp=sbetf}DV{}4YL3D%rqivh+QmSsUMit0OdbRjd3#7_{{oR$ z4tZYf5QNE*&eabyUxFO%w{oc!*&wI+TCvaY70^J`V#PFm9jhCU%SK$7&tLzUqih$x z5BNQ{hQ2gm>#4j7B3sA5jw=;n_Q0v<5x;KzfI;G<9SqHZp!QO7gneT$-VUN}UNQN= zw1}P0)3&Xmp-hK4@CVy|5g#mGcJPncK#men5bILT#qkI_01fIZxR-x`jJRZXE%~f?@Yy#HUgq&7WJ)` z;$X=U#{W_nGhUJs^*o>e&^ z{u`=Mgd|12-T{&|zF$A$UJo)HK9YW~IgZZ@;v}#7buur3d?T_Jj{Y++#7>`*EWifj z5G+jd+=_$l?LlT!`ynrY1fzY>?^R-Rk=x1|L{QWP<46ua zCA;aS=LK_h z#IhU3IY1iGV?iUsJ236oi7)JrE`$1SuX}v=`oXAIF?4e|*&zGnJze*mU!i{SPH*o$ zRWRbvSn;z-zRh}o&6x}7RfDRXrXVSbFQKns-S0Jyf@93*@w|b7R0;bM=|5xqi_Wbj zqit)4=urHT)ap45C;U95OaoGIe z^%wsh!RsNuA9!A8jzxV7N!f?hv7sw4ICteHJFy^kKjZS7qm6qzV1o4>oBGoo*n0!x z(bSXK$*-TG<4uCvVl{sdZ+MU6x&i^#uY5zo!&k|OF@HEu4g=nLRxF>y1N-W^j$8z} zL?MHI-S;p%UWXc`uDrCN&I7g#Kg&x(K4N-hI(O0AWalFcmjY#yg3pAV0C$sT@P>~b zD0xTknfU7z%;`zHDXS5(bUebsX6^t668v0odFOtNn^;jM!sr@mZ97v*=o{}Qm9f9GMV!>4Bv zbimrvUwiE%KcKVDA$^_1PZ@OW|Z(~s@$dpBT-EPgw7&(QuF+WdPj;{^2y)P>Qh zRO)ODFt_i0&lh(MSdy{cJ@fe$$b3@kcks|-=rH{K@*ZjnSZDi{<*Jo9$k@vw5k`^- zQd`OFtytz^vd>#eVZ-kLSGv;m^^q9BEtSa1p5ep%r7zkqO&t3Qlk!{_5yN+Q|1dV! zmxbB&F8)4DRvO`2xN!&b6H51*YSG&hyjsdx@Ya0?a?gs_T~NOT(yzsF9Imp4fo`4D z?Uv@CO!DESkNQVI^Ar0ZN&8Wlqn>-jVfH6{`b^2F_P_+pU?3aG^!x@)lr?q}7kvZ% z@wLQg1wBmFefV7D?MHwko#aVUq`=~*8=NHyexn4kf9RXem9>Gy-cwJxBAsD6@3DI{ zi|e0Te7)Ze7x_cGd^`K{b|08u_?2#t7{Pknc-88kQy^1G&iK0acM$&;O!1xi1k(18 z*@=IVg0?lYR!0)LK*C91s_`ZxSml`6Iqlbk#l_>k5f{dJ4DMsTri1liV~?%*PhSNb z{kr~n&&QCjv~d7CpX@ka`9;M97CdP2%wYTh8Y31=#H4;;_UwDV%e8DsF+U;ii$;0S zr!e5n0rvORVbnRAkH>;=+=M*dvS5C_SA4GZ3V*=zDbx`Ne{zrnWK=|iHft+`BATi-m`JJz9ARsrAI{#JxJ};OjD?c2?B?HP#;_VN}+=uzU-7ZO3UV~Qi>RhLP zjDw^HLWFsF^w@kPC?8vlCp5svgOMa1#?BhtAf;Kme}dpENKac;2;12Ja$kg+OYAxc zLQG^H=Fi;4c$4+|bC9M#H^`Pe8J7Cu4verDKkihd3SP(Y5yp0PVE6GlkgMz*aMbY! z%$j-rgFvYddyemWMu1h%o|b5cYh4SSczF<_b8DsUDJMs<~ zAO9&`b>J3=ddnRO7aBmR22po?h$76p@#cqHQVK`|`j0y*PQqBz9UsfR7cl(8A_dH6 z_v(P``m_F*YJ_0YjKXSa3N@bpFeQxE>dTWSAjd@A)OEiv#-sEL#=ovTGRE^ARxb`Y z$=ZZ~q$Z=s`^_GJ%u^{T?Kn(KYKqeVBPemjS7Y^?=#O5bm)u0iNVf3 z;wfyH9k0s`){H%T*B`;`=NvVuJw+IA3!LbTLSkQHx|pH;buW){I`H&6dDEJ3FZ7tY zv_Nl)K#w;vXMSkPz-Z}=XucY5P#|rPUie%SudguOXpOjqdW~GH8?`j>I>ud)#l$RRC}(^U(9596$XQ(o%o zjoXQ!h}&KE)xmO@X5P_#2TsD|Q6gcv+w36iqa^D-W=9y;?P;W}(GD__Y`_V|UwHo@ zOHDasVX6p!AEc>%&b*_hj`ssnQt2ywC;dua5YH{ni1&c3vXd?Xn;% zomb|`5+4ZX`bO@FlZI$X(Cfk9iy%=zSMfy07${&6Jl+^H2vYq_cjT3Lt?T9F(J2`S zx)cjnK<=p80sq1@aHaL@-b)r4pkN<& zilaLd^sJL-EU7*Oy%y!4@i***ahKJWPL1|M`9N3ygtSah!1alt`1u2n_~lb~9objt zcBd?@Iq3<=vgdJmLHq>7iDj`Fe+)ZES&p2X3wCQbw6Ww0mkHz@&u2#z{)v2-LL7AcsodnyZ)}f?>FYptg>v{ zeCQd>|Ec$STNah{;qihgMV%zf0}~*bJn0+F z&R&>dZ%*}vJ`ZH8RR^VohJyg@hA2+)R?Hthud1-R3k5ZKAUT+%e}4l9cpmdU+t6?q z#wSCZiBW}N9L7Vx8K%$_{c1d4V99|`?ym+~LAcA?_6hNGAic2ew_KnFj0#a{IpRJ9 zvldQSP^(N~xa2Lm2mC%n0E=jrN%yOh0z0?6Nsj3(Fq+P8ct`_kEEJkv3gVdGh~>}7+{z}j$Znm28iQcSvEHC2bp#~ zZViDVFfraUTrr6aKOZpXmUgh%l|U?BHRq16sITgE)*8v@$N19Eq zQbB;IW{C=K81@|B_f%89%dsx1Acxgw#Tu=GjlJZMWAQQ=KuS(D@Y1qtG^kO4ig!lw^=nm2XO>hnabboLD4~io%*wEz&Pwa zYa>FoMGx?Emc*L)?f&Oq05^91LENAajFrrY6xkz!@h08P;If=$7{Z2;3X{5(VqyI)JtvxGsGr!y<0n=!q_>u0_tZF=eDS%~hwQ0MGXjh`Qk_jvun z$HVtygV!+g;bF#!ccEAv8@d8x_rFYX+T(`Z&yiTQF8!JV2F&)773uKMgl6(mNiTN}8q8}42$cN=qsC0%eX_f#)&N2)7 zE01AzOqU6G=2h$V0e^=NVz1j8G2WUfi-xOnf5&vuS<6zMGV3EKes(?Ud~Q1^Vi4Dp zK#E}Y4g3DAm|Y;asv(0ocoIg`BXVcv$o}Z7;FsM+nsnU?Qpe%y6V|Ver1UQh z1hovH-Ccsyxg_*3wmjuq{Q?!J*ip&KQEUp*+p{>YjxfOj2m0$Xj|@SA+(E%N!XIJy zE?X(D$nzlQ_+I`uUn^>h&^^|Ez+ElB%Az zjc=gexe(om8X}m^6~RS!=_SaK>nD~I%mjsCu_O$h|By79Pd#07D(L7A0BH}QSI z@3A$%tNwIc!xtDnyi7u0asjg^0JT&01`D82S~Ta)7XnZ|(Ih?4&WqV|Uy(0!O-y0u z({;mRXuf%XCnBBt@1vTrctHi;*O;?iK*0?5ytskWd0<8LEAtT_DKNdA=sh^O$k4y|&7siL1x zbTwJw`2>n=FBXZ4yJ2_~NxDrnChW!d5xz&EZA_qag#C9=!bXnK=*i=;&2F!1=U~ubvLOwnpP;gM~P@SbEbdYKaep_}A zc<4~kouP7p$vuDN32+85|0{Lo{_pVNbL%WH`he+4PK%uZ*FpIuokO{pAQh6%FHh<{S+S$B=frMvmE&f((_Ke<$pH>)1}Qyu0OVh zQAUsTZV4;E_8TH6qqxQ}-W2l*rAaljz;JVs%d1nGz>I<2N=gI;-K#A&Pl~T&_cwF^ zCKv3IxP167h$V2pq+|0IdyemWb;)^67SaM>6Tf@)ZI1l}yAaIv zS>p5`{OyFFb+F%q;Q`k{G%HfAGZBA@}-+KId6k}8E@O1A0!6ir4LFE2_mY}JF!Q=w5zrX2qxY<|*lByK+w z+F+>u+E5_HN0{Q_up{EfBUn0bAO6%v928%BtH!Rg{=J1p`_{BoB24StQBU_W0~Tnw zvhn{u3te=-FlfL2x>*nY%o%QHm#dN=NUS;iRC=c^Oh(Ap_98KO-hc$*qh42NdZPIr|o!7>Ct|9Jfb38#yb9JD(@d7Atl zo&mgA{DrZE1q`0mSJk74~vn!ERv{qhSeAJ-Gy$4hfDen4E> zf^)f`1Wc>_e!Y=N6U!6NKimn5U;2;8LBY)$cQa}f<89&NxCHav?=W4o-#>iSgmDa} zwrl7M3Pga4!({eMWyhcsnRo$9p)$x2kZ42By@a7^7p124qVW0({f@`}>foxtbhSW4 z_*DHa5;A(5tX^KH#@9|3iPJMr_+Dng)2IxznQtG;7Z;1cB@zkID2_P(}br0 zIa7#IGO;+L$I@sfJ&gyD13u;QIpGgbf`0WqPrv`s`O+23fg&lBb?rJp@h zCJz#6v)`0D`+)SZ{x0I6z0lr=`}&IGx=z<$e6}1w0Wzt=9v$Jghp2*#yz1vfkn)Z& zaW>T-IH(Bm6438h|32IcgFT2r;VYSViF|nw(X7f|cE}#@2jn|f-Z_0^7#2lo(O;ux z1zADoKSVZ*0&mJ**F|a{SVJ8DjHd4*C_E6Hm%pe8JfDKw1xMvU9)op-osJ~TVL1p# zCpn-a;4?W+>42MWr* zxS+HGh_ehyIH)lLp7zE4av0wSeP}avvoeK2v}AQKakV^fnpKr4D9wd&6YWU_1vwy< zk|N@A`3son*dzN3X#_#f6&e(Z-@#kTXjw)t2+M-L6T749d#orYd<%B`wS^B{sxk==L95QHnlIg1yi zLKt>^G{B!4r0vg-C_Nwm3dx7q3Vxl2UgUL`XlE8d-W@LWy;M&?^~<-8Ni@nJ_Jt~k z$hn}{cjy$4_yEJ z)%osfw=u|;@w$9Uy%^(fsgeL6*#!XeDFeK!PievAH(w&!m5pJ@>D0IkjU`adGe%@R z@)U%ZSsyj#hyjU(ci7B^nW5uJSJBkIeW2uy(|zx~N+4{4MZd8(6`&S@FIYWZKyMvR z3du2L;9k@Jgus^!7FYGC6pML-bg=U1qVP_Dt2McGgoqp_AAI=2;)yNvIbPyBX8sd| zX&1#^VkLrchj2t4uobxUm6J9)_Q6UDDmM2vD-g!$z|7sKP0=O3Ih*ZDNsc^Kr)$A(&=S$fTPdnI2W`6_1=}6XV}<-AituS zMAI3V`A+<^X|xe=PkMVXvWyH=oJsL+m1_if@sAwDbj4vD()Bax;`;pg1kxy_tARt$ zKRO(|R0IoJ8RN2WvoMXf+@AcSF>nEQ7`{@@!fLfjT%4FP@I89R;7;Z(7&QE|{cNx_ z^vOFN&!Dpg(k+jFZw>zif~-EW2(+z#e(=4*XCg!lB|p5nQ+&M|=E&L$lhhHyQeRKs z-F;QiIoR1lnqv^dzxYPdsG}l3^eE3wh&s{=k{a(wt0efsd<*)$gf)*q#t+D=FyH`- zt9AK$#RXvj=g|)yYF#in=850XcP1G0VAnPA6fxiy(zv3wbNzG8BGrQ)rSl+?>@ua| zo>b@>hIo%U#Q=%G9joJgmcaH>i*RhJIC!G^NwKK+5G+(%`J9t64RZEH-(wGd2CGX_ z+LC6@!1!midk_BZ12wp}`&~__V6N^DTG8-SXeJYoJ8?Z4JQ-snT+aFlZ50+}M;K3o zva8xS_}Ecv9EI7Y)Gj(efO&b^Owph|^9D;(&ws-UmKe;vo3@y~!(c z`LOJ=YK+Dn4VWjUYTdTm6TbKj`gHpuf%XwmE5id~FlU+el+(wzpkTVk^;5wlj8#hi z;`#6-)V3&FHn`jZ8tye2ctpyB@UXR8<@2fV`FO4stNRmBK&5`)Q{D*XW}w5>_RSz^ z<^)6M!2nE8p0%GABne{zIa9{`{s@5nY@E57N)oWLf4rM{^)`s4Otd;g;|BbQc{}4c zKEuRY5?8M{bgc7mAJ@ge01#RpO72Rm0CSDR#lO)=zz_!c)Txhm*8N||(Be7^k|-Ei zw3cqe^rrI5x*uXd=I79E8l{Jznqr5Z^q2h*&6&AH5f}$i$)cCKx_dxY)R$jk?~dW~ z0R^e87aS(MU`nTOMWd7u2t@8BI(oc^ZA5F#_|Yg}Ua}(m>cR>r;E1KEc2H&}ix@kRtMavo790T&jk0`ZHW>3YkYa54hiqEp9-M<}mGl}Spi|<4wrG;7! zl>A0+zEoK3*=>e$y^n|$%&$N>n{Fx1za&MOztv6J{Zc?c`-D54q9#gRC`IEjJ&5Rl zUo6!ZpCf`{c2*1dyP+4gH-StE0m{@aK_zzJ4sx*d+t}?cYjn4CF~e(DZR9}e{A$nM zK`@6hU;ZR@3Xw#3#~<(yBNzE*Wb7+FLEO7s!^6~-gEUi%ebcw+09--piuSALzm{md1xvN`TdQo#<_$Fh~ph zdhDR#E>%!Sm=!wmln{6~8j9~f9)|Y=!cRU~7CfX5i!yULl@~N&P8gp{_p3S(9ntMr zafKSX9nK<&mxnNOFlL06e-1?VE3{Er*~8H8u$Az)su1$szTPJO4S3>gAI6F>fWonZ zllN6pVf>ltdynlN!@OKSXNL|J;G6$lQnUR7OiEK^RT^Ig{@kAOjd#mHo?M-zlS=`R zwsziYAPS1vnj2g$ zn1K*(%?#=eE|6~(dq26LAjBvy_pLm%$e_?Fl4+7vhsDG zUXW=AX(Wr9Z&VRj@OrsbSLG6P-?OLx){G*If1RszrW}P?R=LU28bu(Zj$7?=O$vy! zpDx>FNd`Hu^_o|BpMeoezQY7DtkB6acI`ap07$;hPM)239%k9hi0n5v0uICn-x;)g zfv9%(!Ii_;*WXWW$5Z>VL6|O+)-8V;5b2+NGm)JDI!zn7N?i{HO}Up0d1hjv?Wco3 z8Glp3{F;ciHaa_C8t^+v(^Lp}1ZQ@q7&pSk{okhji3Nb=lC;cm+W}x%Kc`lAh!Z^3 zx|KZt76CbH&$WF?pTSJV{NMRJ^#EQ=J~dmL3LTHh^|MX%f&j4!5lyT0?-6w8l#{}D zKpT_FP);rhfE%MZ?sYB&y42aO`7eD0{`%heE-y1d(Fr!L3rZGH$HSla&eR$#ig|Lg zb)PG2^oy3K>*oVuqDO@;t(L>QYSFy=bjBe4>4#P+U0dk!HK;{bdkQqjQ#|WVBn7d% zJ_f(&DFwJyR#$@D4CFi=>ZLVIfE7)7rak*Rfv-r3r~%J= zJbz*Q^_9~?Ph?Tt$055Ew@e5vRCXxewsTjd*K zP>4ZTh6^lDxawwV@CxRJ+)>Z$W&_@qYW67=3otH@W94a!F31y|X#9OjAC!qs(;Be8 z1hn&Xem`BZAdY_S=@lYV_*C@tw{!DNFq0{yg^PF^hR*m;T27^}e=faA_n|2WSdImb z>6;nBs_O&8R0Moyq@|)ABVoZMrzJHP9hgI!G5b;>7(^N`ROsIF0kI_!C&aUTft!TM zZAHNb2w&uJE4_LT@^F2@3O*b#xZq8|>Au&{&-VZkfhZr$*K+faXPgG^^z3t|nr&hG zd-;!2V*sSa49org)>S=o{*bMlY1tq*=Y+uKXNRm7)6J#|eK3QEM0a*i`g(KXLpiZB!k-3N>Y>OlB z6VvknZp=Zcr!8AS0%Pc|cVBw|>iE9zc{(r5L&H5Zo}L5AH@~f|+%?1Wgk6rWCz?YG z{y{i@4z1ti^lx5M0e5X4Dy5f$|}z_f5VS1LvXM6KG?F1~mJVr6;G z@6T`qd1*Q+@8sq|qF!rz=8s5#^FI1b?ZNVVb@b0L^9wD0xDSBF6s z^B!EB5(P}77-4$Y$Osa;?=`(F{{oYK4in3iwnM4t)!7$&h|ze4)K6Iq6zH72hgXyB z*ERSk0wPU7M5IKd2?(J{hlnUhmo8nT2}tiD1QC(mL69oad+#M6NQWrB_Z}dGmJm{g z=l9O{&aCy#dS}-B31^*s@4e4i*M0V@kRKT?B_@2Dy}EPOVjpINqo)yBlh}v)!R}7Uz27qVb3$dHSvJ@{U=cG!^_hiX@_zW( zmO6tkKrIgkGQ>P{q)vDv)Tv~tAgAGW{*s)bvBN zEW7C?55sL8M)fx!8n_olq9Y$lgH1bH#cb1beRu7!irp^x7n4WO$QjM`b zyy-pJubHD@bi?r>`+HT}Lg6W5%^2$H#cC&?%KV3bhXZVZbDC)dcOJLWC&CSygWq;@ z;O+Eo2S%{Iz)t(;bjz=SO?QMwG#>+BBe_WU#u>kEHKaT`!qO08vdQk~arx?)#_kEQ{Gb|1)}82k zRn=^MBSM~x$9_Ea`J@vQlj{K6yAm)WRHK_>H@P(5^Hb*yGMUZ;dIhx`^ir`3;;+Al zNoSFH8S~Ekh^Q!z_iDQZQReV8+!=oK99+zqT-h>r0KVll23Q;{+>v0F0>)F=Gsa@e~k)WOjvx8IbfSc zFMEOYc?Wzw;+Id)g=fRWBuO-#zUg_H%|>{CK@_MUnEB`|U`Q zud-g+y?eKa@fGQ^B%Ut{b0D(6DV=#Kv?ph8?m8wwAL(h6VUNC7VKP zd++M4xW&CqJ5x&M#49irZz3uioAW>ay0TDr z6@XeB`Y z7RHCQza)qGOG-!WXCeOg!urQ-H&W^V(cemKb5{n*7NHI2PsnvbL`!*d^g@_|5giQP zU+Ic&UdKnicNW)1P&ccl_gi z`dyXwflyB`?Tw^oIcbVVrAU))TaOR6G>1O1WWU7MSxRMW?4L`8wgaX(z;7M8%1b>I zS*Gn97^WEvUL<+aDM_$i3a2Cntz`3igh%b$%=hP!m4%ZgW43G+W7o{Y_Rf&a+w(ilC@j{RMhP_57VwUx4|;Z*N-J+H4yPuvlhj#UMh z2QA0wvs`@K&bb+viq3|t$k5mQ=+?t_rf>*v0(Or%-nDu^TN#HnICDL%*_hh0=|BW~@$d{X0=@*O9~W8!7xc zGOGedL9S=zas#|lZaB;uIPImTGHB~_7bZ+oi3uRvRi1VQ?1&B=`ZIfg7FHfpfI!B zh~Ztx!m8Zu4uPtMt6$vAZ-(Yt}DhMHl01#HaSC|;QsifG+M5w( zlE?>T{h4)@i1fQ}bP_)~y!xHVU^nvk=I2j0?-VMJq=~${+09DlIi2&pw8`*JKeHLm z2Z?tD(Q4~V6>0ft@7~{;+=X#w^iNSac>%*%rF;uj$8KBG>wB6z7j4fB6e_SWwJXdk z$41IagE!;NxJHh%9M&LP;c*13iHnL3)scQe-mkA`Nu20>8a@Mwmo$e`n9bB8)XCo{ zP5XN2gTeDTRD)WE*7d=pQ4aO0D)vc53lk)!E4G0*Np*|L$)F)&AiIqgLA&I@hg@m| z<6CvvPsS*=8HRTPxLx7jCl$1}5H^=Tr+BoFwz)bi<8EH{f+zdk~n(&vDmZhcvNYvR}Y*x&C;!ZKUi z{jnljo~&Vbylq$<<3&nc2a{*(-}5cn{>(V$ZE&aa$_s|~_qhz^#A}gfY0};70fzOv z55>CYY>-k(m*Zoy%GA^l`uB}~mB^`}E*S^&944W$V0RI9noMxR+EziDKf`oPGwjMa zdL|Edsn*$$_$^ytvAV|={vh>n0}CUpMn#7HI>KjRy<%Joza!G1#Q&Od_^JE@cC5JF ziy>5etN>8;NU@)c%4g}GwJC_;Q+J9=blTYN=HSe-PPYJxN{)IJ--9^K)hCZ(6J5{t ztO8u+!@qB-4lSDeQQdEw_0TAWllxyq_ua6FnCALXJ<5`&7C6{i>|<5rG0U)mCS_3f z`J2=vcqlfF9tGDvE5^RA@}9(Z+Fqj$^%;Z+Z?GWvUQ4_efhfgS2 zKejZeofq{Y=guY@(NPHwMXm=i&^gW6x)@rwW=vT7y6#T%Uw_)+HG7M&i*)htAKT~u z3LqzqPa1VKvt*XiHOWk$KEYM{XwC%^H?}({laq`(Go;1dSxNg}>@6+C9*Cx*F;8jt z52np6xNL4-DTq|*{c#m?<{TvpTzS!-*Ubw&-d-`}@zeHA5nGID(i zd@9i66%7{yl~Jk~w^R15i&}_ZB8IwIqTwV;&Svm=o7s}-66abdStgIIZUiw*xjRnX z=W!hs!r2k>#5i_L`ZU&P#iiGtxXr0|cCclfD_@C0o+dSQF3Ij~&(7G4hNc}^R!k3? zj@nM#ppF?QE?MhKgD1eX!+q5j8>l^;UYcRvt)6{IWc}7mR6n0bY}tBBf2dwvTaX-d>*A zEBrCzk+xw1qheJCEC1xTFGnlW@R2ns>B-jJ&&E!kLjvbznc2{z*!IztTIVYv#DGrR zuV;olhrEbW8mK8Gt$;O+XEbmRq-B9CSQ)evayW^W8~z6O0&7UfO9X<6FG*_yE~1?% zW!h6|-%@Erf85l%0Gf#NAE`s7>xtj1iuU>aaFwO#{d&`g$6Uv{YkW>B;GwHTIj5GF z*kP8qr^H=7Fik{rttC*=mD<=m_xJk!ZqLC|8U~qU0dF0P<+IaDMyB*3QvqM&hKu-M z)(?zzJc{u2kAk*a{!hF%jk1Yrq^Z$tb>&(krwTjUg}kw7-oD|G)3$gr-oe^0_Wp%m58`PJwY)zzJBy$MOzWLC%1$>?u8qk=@QT@T zPJqSoOwLf!ap;iWrJj<1LS1cOSXRpbO^c~P)madc3}cyi8u~1mGe=FC<^H*5T;3pM z#&5Iu9tC-)ptVhB`egf)?{O;21NM!;?2aE{FC|*%QY)6u`6=e#VK2Kw>$u(+pE}N& z%xoUVw{Vw9r12n|mM14SEL5a_>_xKl0Qk0;Qxg%}lm{v|_b>JurmRBSw;c}JZae&Q zoVvtgS&bGE=A;SzhK|)Z&~9&woxQ#(s<iWt5Tq*Wm*=)yx$=fWmh;d+OH zUC-YvUcKM?aluXaz}{?_O>cyLGnt`#(m0@Z zMSW*9m1Hs2?5)Dqt?%^dqId%c1N~Jo!P`OCx0a)SpkLO!bAYmHdojm4g-sgb)D`xy zW+wuErK4?vu2g-o0BtX^ET^pEjp{cPY#S_oMp3bm={r1hRcnMEOhS^v6MTKOXOGH2 zvE?g##%IV&lWZt>O@=(f#@};7Yic`FJh5L!Kg+mfgT-ITBkv6&%HpJT)7oD9EP(Th z(Z6-`T-UzV=fzYMN|0=jK*n-$42*GjnN;MD_^+4F07u(V#kEzs)x^~+hD$*?0;;*LJfQTU=Eti2jLfJrFg85Xh9VQ>lZ z0|EsY0GACO4)h=I>$-kSedRIxPk+zhQ~SYIfGSu<%x&&sC99agQ^`qaIE!SlN;;%S zG3q$wRB%%^Kw>ZjJChVQ=C1H}V&1q@e>J>BMK!2IfBvg!LwVHjnxGxXDM&o>*3-Rj z!0Hjl+UX(d)8Moxwdd<0x2p*39v=^vh5EDS5}qe*UkGZk^De}E9m$x7-7_ZgRZT44 zNK(b5`3U;L%%7-4Lnm;{g0j+%U5ZX+iwP$;7U+00KAI{TdbzhW3KEU{b8^usD@Q9S z*7qCZ6grFUx~HdLmPhQq(35e?0&4K2(xtf_R{Cmg zC7^<1I56%$zq485$Lt3p%^Q(wy|5lv+qjuUK9waXS~yHcpcb##|7c}DLv^+FQBnoE zZ@|jl$tx0iMQ;2@7+t!xwUhWGvN6TzOkQ2WIJLgmeM#Bk&xuY)rfX05jCfuF>1q0a zi1#*GeW__Q2G}m8-!kxWdG5fIV6iikn4+X@iEpnJMlD!&=sqvMVQj+XKJ?`;>DH21 zV0sqDT7hC7T}M zzy1=Ny4sUYb>J^maLj`xNAPD>@1{NVoX%{Uja0pP!zK+ zZFpRDb1epH&U~(*VSugf#iS$IGb~oFvrj_I%v2pE| zUDTm+(%8UX0gcMpIPbGB;ajf*Dd+iA{!04ZR6TBq zN-CSDdc?c$%xy>R9b#b`edBdXb@17gOQ-p-oq7sW$kqBt?1uFN(S9G5PI+f6Y`}K zDdqahUaWSYexS(TiIyH&!d+&iLsd2_j%P|4DK&iQ?Xw-4WiAmwUCehTw1Vjet-K;}_$Pje;kP?Ar3eCHZ=o|vv5;~c(v zs@`X8>B!lC2FlKVDw=%jb8Z4+?Q{}QU~s$;ZqdQ>Z15zPbKby3OnDhF4f~ZnEw*JV1{gsa z^Co^YVio2K2wOTrT~~JC_q#{AFx-o%vZLt($o1>y|@ujTNgBQ1w1|Cpp;ye_--}%X{1`1 z%J;7o`~L}&Gf8V zx-H#XY0e<&q~~At_gi*}%CD_slt-*ps3*PUj>|KF3S10^RD4ms+0C@px+~W{Cbpa- z*R&1R1aDAVW)Q~xW$AbK=-rMC;~5!Nl(>{DVJ($>8|x?fV*_I$?Cq(eo_|c86A?un!Y2Pi;uvr%uY&{BnGZtdK|j z>ZQF@eK%OlZi*aB%tMfE240kqQHTZI{nKLi)E(Yr=ej?ChVIQ~GR(T&`g~eDqHZ9l zl*K(z4z*seN_qW>mq&qGtpbRqP$6MrCCwKeG+*s8`%O7kkqz8`n_g;BrnP z->IUbEf88@c;^r-{zuJ<`%QCng@Xk$tkI@$;6qec()E13Jq=n#$=^AYd)lBoIivRX zxj<-&VdPg6zH-t#1GVFT~aRHtGr9G3ff!x3pxHCN|ui z_>}zo3m}jhBt07!{4FQ*S6uFo(1ePpy0DCxHWjeQxa4Kw?57_;(z|~3j4t=}3@*=w z9L*5<_trh)zqL-K-XmZlHz!8`@eYzBDs6#p!wQtW3zL4Z}aOwc-4!fe>>pBWwNW(jEg29 z(NG_MKFoEvJ~p8s=@?)Cxc}?NjAx(oFIz_C6Lz#@hu;2~9&! za#Ohk{}gYY=2~zGc00G9XQpqqxac^_d zm{VBS;q+dc{7GbZ<}`ZMXYQsQ>-s}NUonA@WJ>Cc)n+%K$|r7GU+Jgav-!iQM&v91 z5%=lw9_ZCXP-Foyt&n(6o5AI`}w%2!W!&#N5~L0iuvQ0o+XmEoqzW=qsGs*Ca{zY z#s2og{q=3}`M=7KCNp~FdPA1{>n9%3N|GD-`f!qX$H4AWLPzJ3nGwOIk8kTz>(h!h z7G4zG8vY^~c82YKK{iL+0@C$azG^JBS0iiu4YM`HAbkBC)^har z`t}394_#le$VGKXlXRmklu192xTZ7T?nLi%-?f4A4rPY)IdT%kyqTD?A7cl3Q(~by@i_q9vj1tN?d1v3wHPs^5*j?y%F#fcZ zy;#rynOJK7von$|(x~!m-4?n1ZrJI2k|vms#a_qKteT!_KZ0^P-EY&wi+TH<)5eo_ zeZtJ7_61E6Q8qVoN^0|i*Vs2HVm|gIApbTwX;YN(Im%QztcjtZu@>zfkHDr#7_7>}8 zk@e6gn{&SkvaB59i@hXDVFk`Tzp`g+Rii~b5X~B2G@Gkf#xoMnO5kzkn z{EI@yLM$I$)usb&D(201YTCnZ{*3OlRqkZUWEOHQRv?$0I8Wk87UBceA?17ao@Mdk zw04*>L}X=j6#8Su2(r6=A%`%x;5FfooWD^DS%BH+m#lMI{mHD!+U49mo(v5ee@b3MD^c*dw-xt8bD>%BqOOj{?oTbk(F1g`QBP^< zg{w4>u1J*~_o!VdSK&7JaYuUGAd6EaI~P`(+k>_EFLImuuzZqpA+`YBafyF$|Sak-tnTv4O z$jQi=!+U>}d)2R0_dtPN18@Ne4|p}hU^~QHU?AE-bg2ewU|;X`6itX^&v`<)Hu^F7 zIITTjS}QQiIBMQ(NYde1f$ZV%nRp&zXa0}>&~2HO?)}4OlrqbU18)4oo3n55y8LNn z{fN}P&Zu%Nc1N{-8%r>#hO;Bck;N3KHsf21QGo9UM)LKPFc(nTOF6Tt_4n!z?aaTF z(e@iX9;(xk$BJO{+q@RXV`6u{id?-VaMiZfa($2O=MnB`)Vmr}`E$;Yw1=%p_3AWI z?3cQL$STCe&IJ4kL_I4}XwwOz5$g=&vlz_PhYN^XjIVHxuK=h%&=IyDS80`0PHaBc zNB4^jy^KJM{RNp-v_eEHs%D%)MJv~@^cLbDGt`fh>iWw}z=mF{&={8j&Cl1bvRkn{ zn+_8E6Ie>dE0CUOMG;{ZO8)|BJ^7l6V;8<{6Xs}y=ACzYbXD3JYd#4x;b|FukQ1OZ zm0IcNnCxjw zWDfUT)ugf8%u@xJZzOC;n$H6FGe-KH92?S6-sY#Fb4d>PlKBqXXY(jkb=Jr;^rRcp z+9|_I6{m5BpM@|q9jBdL_Uq`JqXH+ss!Z=7fA>&%iK~fW@_8iKc0HsyI0!fanJhJ0 zMKz8XM znd;W!MzLI?1_J00oP36A@ytBpjT7s81|weY4N{+`J!#5-jGJaCY{sU;>w-C5AaCWB>R{$w=Y~w@3hlF$M-5kI&dUg0}(eK<79T&%p=Q{*V&{tJJVd8Q1!ug zT&n-1Cr0S_@*8p{q1-|a_yD9y;05}IY8-4}OawsG%?U~M?BPCnE{OsP*gS!TL4p00 z2-5|h~=?!J}rFcB&Sd*PJS|DjHu-@UH^ z)vv0KAyf>{H*2dlYmEhcxGW$h`VBrWc3Ste7dr>}zmq>j(^O^KuU9oqTCaqw0$vj< z$z4tO=^OTx62rV=gV-^~X08O)CVXZb>TeSZXO7UMn!VQ6o*Ix(bUYye+7!?!@A=B8 z1Jg+!=Xo!#67+Zi8+I|1qkjgDcfb$E9r8Zg3ygla$E6q~!D|>K1X@0QNA4^YzidBE zTP_t#Wc%kNomWgD)#AvpgHLOkQ9H&%YiFflPjT<|zc7ebB)wv0l8y2n6fgZ`rTrV3 zet}=2^l{AN42JK2Y>|gbx@{pUZxXhzKhye-C3j=^VgMboD}!ouKqX0`wU2XCf|lVT z3fouh-v3QZ@66)c8lrWZsg8F7OCrWMb?nHg>GS&c{gl$PzLWB1Z3g)kn^~*F7{h7E z(K%iwaVlzubfajHxamD~*>0=%tWTh5;uA~6&-3T?EZ#|(-$mQ7F{_889m*|GeSRi4cb9_V;iHOp01 zRiY>=BgZ`_dHARpNyy;E;L$7y7_PmZnc@m+5pRIZiW=L$&$*RlbU0QZKw3Y?Jyj#v zScBeGg%Bs3#w5Nfz8c8AsiggIriz1^(pwX98Mg2BR8LB5>`Y$By(HTy;oIeR4mvMm=CmA`0fJ$m|-CV7gaRQuo>=I=| zjyUaio$CI>%fcpAJJm7|(LJso0y#N~?If<_PdEXZ*xiS5a!0QMUk3KfHebxG>n6IIYZc`M zS!oxZx4}uu;2(bxJD$#S7>Kv%I)7Yv9fOVX$r-EJ9gH&O#AZ`LcjbxF|;g+VDS9#Y4=jU(*D z<(0nFnjIhz68M1OJ8j*7d@=Ka$GnmXuT$iuq68T0_j|Ec1_u~9=NI&$5dh(JF5IWnmDh7IP@Y*s4}TQ8shmDJvseu&E%1exHjG`Yv7 z5y(D!9%y;2(K@)I2$iqf#HW7WsZWJS09uFJw%JL_t*)+!SJ_z`b_fmOc~21MWE+*G2Wb z%2Lpu5(=Y|Ks;5BB_GHP;a5xVmgk_1IkHvWTrTW*l(C<#| zReAW}>2!GdCEH9RiB{#$vK!*Jd+q6BF0w6ZT*Gg{y*kJ;>0<~?(j=XRj zHZ@0#R)MwDYNYAY$&DTGkfyWT&hAu#x;Qk}Y(oJ#xf6!?jThqEB@$IR#>i4$98D4b zZHVa((gX91=>Xq& zUyd)MBd1Dd_RP^L0i8|2w!NT>#PvGKqk`19R-?VYPVqcu+$&Cjccvu$26i4ECu#ll zd?`Nk|En(7JDLt+Zd%;VUPE@Ki;EUsphR0jK#oKe%)?iJuRrARQufkYPVTJ*-9;-C zc+SKHOXv8?*0Gb&>{T2>WPeOTX~BBQt1!nBJ;#r(QiQ!J_}Gy@>8`|J`u%L&&r;>e zJ0?gowf5}~?DLFrEZb$*&Vq%e1KPjmvixY+<@$?B_P;JHg*rMCBxE_N!B9DVdA{7A zRC-EG{qrS9RHciSIQ?^w_#Kz`8wNY*SaaRj0pHK zj&SQaLAYS-Q(#{9Q}juUY*_V{NxWP@@?{9`NB=jHXFK)yudYq7;k6mFYBXTRdb!ig z^C%#HPXi$S_f(@U;8ndIK~~7N#Ue;+P&Tlj@Tzsb$_eMRxYEP%QuliCJV34hvTQpZ zuxNTw+xnkw={626Hyg#)Od|<0HhDZ|1y#-_qnb&VdjsjmdJ=>9p1RNCm!PDedJNjA zLlU=y8R?Cia(BKY!8M>J8vdU@sn(rM-z$wy^cD+f!u~T-1JekR|8_+pRfHMJ-GNtt_#)B9qlSR-2q8PH~6-qd^{#I}v1{w9KPec?F4`+zWEQBm0pRDN19I$HO zec}W^1pQ03Y%b6NVRqSMu}+a_vR5}-Qn-hO?R#OhJ9@S4K3ok)=U-2iSWjFq5&hvhsmnL5Zcf7)n3tsXN79Q=k}*V5v4oXH zlewd%8lHJLJUne2SC<@jC@{QzaHdisb@uK`?K5o~9E{<<@IFse$Vj!acGwNAQxeVQ z_;nd^(gpkS4kq13tH3?iF1MlgDdrE8)Rtvdu04AYXVOs3`T9poD~T7NRpanc4uDTc zRm~@u$icglNYU2!AK9*kYNU6AInN|-E#DNHr44!k@U-f9lH^K;n*}rX(SmOJ<3}-aBSGB}36Ju4-%X_GYl3D297tNO1A+E+(k5JU#jb@$tIn=R<)fyMogc z7-l~X8am15ILt8jEOq8{@~cn4I422vqG>OJs;cU<#nvq%M`bWcmE{XLYfv_zixf`X zOP3wcY-w`Lq-cXg#E{5qncYL+-rzNMgriNV5E78jWza*^NWaX#~vGjGCN*&fu zEfwx8?>8>pe}(q8I@P*&yM#NiE_FUzr{nckwTe9!?$mxTreHd{vl!dvTB??{>pS=Q z@PXwnEulO$^ORvLI7l7KtpW-1#b$vrwKH@k$Cfo*GdZ)mG4lV4V|KoGYvGq%@jMPq}=wn;Pmk|x9wKEdCX?g?JKjUT?S zhmwGJ$gkG4cSEH63Jdk2{Wi>!u4t9~^1w@52u70>=l7V3Y4xw!&;#Ig0M~HxFmk}H489voi8FkC z?8+7lI#Yr?r&*_gaA?=-CDY2?4GU*ljltaAxQV(WKR3s0M!V_x49jM3#=f4*l+i<6 zrM-#=RMVcUT4Gq-qh#Jc#f)aoQD@SDTHfN<+~JsR_Jf|lzsRnsWRwzX*GjQ(aa=J| zH47?6h7?Q}B;xu5=shCbDCh5L|HQAIg(e+TLewpV*28!tD?GPf40rTY{y2 zdD5a6>TsW6!k+=VxwHqH@CIIHL_FO_7=)Z{U%{>WcE_r(K!d!I=pDQ5J$_n*m41nh9yu!E6>?&+-+e?Z=8?2>RFxR?v4{2 zt*N^2H%Xi$E8*T424@ybS9-0ylpuOO%98yn_JAvHZf`-|B~iH}u+hc7;yRuLv6iY_ zBX+O!!Z4yCo~OBZ0IZ(J{C|TMEGqG&hfpFW@%y)i*cxO}>)E4h?T&&u4HuRF0PI5E zz>f>He3lSRq+(=6Jm%_G^ueoorUD$cW)3$XpxszboqiHY;U(Ez1%9Lpw z-fo$YVjm>0UfFF&em9);XDlL*{R6vBHKY(?YJLa8Q3MM2OW9Ol7yik2Uot6saO@4A z-a=DZS0D*#61Mj;{1a0&&Xv-$WF-q)KPPBpXT)*lUW=K;J$cc7~Stc<1I6$ z00YReP0s8CW}&X0BmF+VV2cAPI-wxBo(fw0ym;(2?ktX7CI3I+Zo$jdTu>roSV?Zg z#uEd;SpiN+E)>~ultG*>Iw~WK%y_%e1G2Gj0E1zW-=Y|_n*8cN%vI5i;heMi23XPx;9 zh5JXy&*f3<@nR3Y<|$|u7`;jk-X)Qxi<5JCdp$5i2YmmVz=(JQx1-h+fNdlQ!{YY~ zx*8_sq6BMuHUA)&ASBJ_@h)NEO6V-fhZ>qgFJa%&JC?YE^tIAJ-j=mKZ|1V`&LX+Y zQM20lOv&23H;g{a^c;K_8PzGChO(6(8^1fR#Yk3v_$9Jt=tb`^)X+jLzt_#s!qYoJT{S`miLL zZ^C=zxx>v368xu~xBo?W=A^hKX13!TggpPCRe+NNDnHi(F@L?&`&lUQu(?y8U5l1A zH;4(FLTM)C)(yolVBf{M`*g*Zj2{h-PtQY3fmT}xMeFGUoJa3>1gTkdteE9<7-zLr zubny|B2bbuRWc>p`4ROysLje$H|n$H8gF?|Rk+@92<}^5Nme{4a!^IC206XzH&SH( z(yjAPR2Ijttt+v`zKN>eTg_dr-#!PKlM?>WqHais%N~cy{ngS@IDXFh*y5Aqx>7V{ zufUm&y}g@!Jlw}JR_dHsopp`Lmvz!OTEexG-ne>1O9+M0H5GGu-($v4XZ?o{x629d zK>Pn3bSo?};c3j~d=rnQKEZU5ZvZ;LcarNY-gnL(6W%XB4_YM)`b^`0hbUf8fKPq0Pz`n%0-m|8IQQRQ-SP;Y+!bW-F9Y0!pbfgvJ}t{O*gz@m2mH!5^|+@qfZ-!NI^1$q|?%R*{_f=dB2{TuBS1!q6BqnI?~EK_MorV8~PLe4gjAAn=C(cMSah7I}t2GRDKFOhzY z-qHU1j~$w3Hk>!mUrz$#Ous+4-!2wBL87wj3aSWLDdtBuo2`DvAM>7dvs%AJ8d=LM zLIlh|cH(h-8#f{N^iu@Op|aWLzqW2;gz3);zbX%$RyYOnuUt3%3)C2E=-E3z5+H^V z`b!2EDMbnzK&r!_NyBKn6FTjLYr&?|f`_j{a+ze6z6`^9-An@_w%Oc*S4=1f=1-5q zpy9Lwii^R(%?O{!lVI}oyQrJ(bb$9}k1TR|Wo1JU!H{fllXV*<2hFu?uTlK(vi+D( z6gwLHADlij#&<&}Z_I`SHUqJ2%@;PfyX41``)vE1`=kfTkUXtmi(bnSmfe3%-yaU` zHR8ocmB3s4=fGof3`b5d5El*RIwBz~g#HiQ$?z|C-gug$V;>Lt!vwyj@Ze@Lnraen zQdBhF5Bd+2oa9{{1(s<)G|6Z%5f5;r(By4nym**dn%!v>#8UeK33KS(2&_`^aqb zFEC2Jn382Ry%5P#%K}pYLw!Votg&(NUWQswJ#$E)T9Q#7&-ueGj*6dW@&qG+GW<=KKXzuiRKR*(WpkhG+=ou|D>vbr zFDzHtnCXTtWm53X>n^D;%i_G+8mOitf?>-wUsg`mW1 zxC-V$lkbbQ2Q%S@{!8hH&ch%(dlrszR{dH$Ilq&Vd);oZ&+IYNWdhIM+lL=ZiCo`T zb|@x5b+7e)+Vt4Rot@+)Ju z)r;ytT<5B9EJE0ahi?7TWyXx^weYiD3%{Lr^f z6Z;p!Z!B}Kt2}8Y-km8Pt|s>m{Vmo{4j1@Jnav{iuK@^@{|5MXfU!O^17PUj&w8oX zwFfb|l^iZQ=wPvbVGJ0A+5MGwl@uYYx~@`O`7d>y`e;!!?Sk$}Co>kI6*kXB_2;aX zXU8iqn*?Opj5+}gy(Y$@|9L7BJI;eubS>yTX;|8w&eONfP6Y;NNAZ>GgOit5rQWlz z%<#VmA~e=5<=-AMij!T z=~w0@ZoIbAsYV;C>x%`HP@Q!7WmAJt_-><)h5K+`vSDCcjD*RYig64ol(4!|emWB1 z?|X0kU#nm#r~PSssqpE_irxH|4uKq_b;yvm+M!O9lC7;Pzmd$)zxci<>zlH+LsO9= z@2i>adi{FyX|Ef?1aa99O`=ccx$ikdYe~nsWQ$b-eJ%1z`m3l=#{wzvWJaY5iOH;1 z?V8qnwFtQ!GK;u`mGkATV!hg=ecn{5%Sy@j&|%ewqo!Y%EV({P9d=v0`VMLV2#4i> zoP*eN4&)V$cUlh{-jkoAuBNlAz%D0kL+3Yz-U8oQ>-X0?%WAERuYzIOh||~Fjpqxf z7GvQikn>uWO6Rg*f$+YH14|E!Uw74rHiQlOD)U7Uk7c^y+n6U!8#hzn8+^Pyr0zdq z5a%q{^-#oC|2l=cZrmlp@RaM^lm6OLD3KDg889S(`v5bVDA_Tllb5ojm#qJe<)Lf; zHD3*#L-4B*9at1~3;WZY0um>ca)lt1?UTE-@@d(}AWVfh(95ZJ>u#UF=u%jzwK>Z-82*T_%(FXYpq|4^c6{xh~8yz>t27rUt2A@2{Bm-NeAZwb4)FAnnS ze1B^>Iveq?_YG(1wPfm}y`+(u3|x(-19McdAIPFgmZA^e`qPYrmVX)xKhC#dw8^n3L~idnrSh*x%^e%iCHM z*+IP8GK#$ehIOtzV(4@|%c+4fs??{B>2;3honT$3(I9SfENkId5uI0Mz2Nf1l-aog z{i;Zrha#sfd@3s2xe>27=We#V_?4ZyxxF*P$r3YPP+4m}GuLNPN_2gginA(3>lOrV z;29LW<$oTGVqOFt`$Fz#TVpdJ$W#MjytE`g>pGR*qPYM>vTrUPO0%2 zwGRIa{bXJL2J7|wv7zs=jVYR`O$HM0n>9aN#@x$0wkt^17~wxQNYHzYPL&Lgmpow@Q@YC2;=j5QZl9|+Fo5tQ-S@L(2# zfKGYQB-;-JoIHIx;U%+tgB_TT#gTY*>jYms{($9 z^y&$IrCiUS|7>2f(nom~X@$CVDM2L}@;yLb5X^eLjyiKq?*d>RDLf)~C-02so#z`O zXTKigDUacrvW9HDNdF>efK654lfX{M5%#(D#Mx7zN*pA3&905jsSs z$|*$RNiLVaJm06eX6I9|bJV8OKjk(KH*FQVD&PkfrL!Aif+E>rizbSknL9br-f zGhJ_`m?17nM7{^aeU#)o7N#mmm3H^SHYBeBZJe0wa1=F95f<3Yv(}q|daR+_Hwg4+ z^bw`rwoncwKfr4iktugMW7^${lD24Y_Cj`%ix(FtiR;J56(_YZ^5IHgiIy-_+;*Jl z61sf=tErNE_+$`P=T)tPWp`seSw6atzvWYSZVKt72+54epi6_{D;`({B!6mXW#(H@=%(C*jmSNF`}j^Fqdi)Iu7wBP_x zHml3GBie}E>W>$~CS0c~f{^5QaCv=(!=1gW1QMp##l{Aj185;I!)s*Jva(Yto9gb_ zRsSN0tL}sxl@P*Pr`?cM48jG_e?$4=TBb9Q;#=2cUp=FWC?u$n?*M~;04DUQwgyXf zG6E#tP?hwaFUwM`ZN0T-fyzs*a)L!@QQ%pD3jju{`@Xt?BoEhuY9yr$==uXOJx5G< z#m#`VU-K&pQ6Y2LZZjp-uhrPDOF6MdxwU20Q6xPcI$dD^%n+u~UyK^Kyi5w`kAGOx z&;FhAlOqvMuHEHmh=LguZWuPODU5u|-KfV0QPUk*P!`(@2}{la~ zFA)K88yjGAZWih}4tgAP`E95F(1ATdf`QLFT|k0IRqZ{Sbsmci;NFMj@F{nY^W&MT z{MA=7@)0@(xKCi|-VqJF>BOkf1GPxCRj3U$x3hQ)R5@J^dc|V?eY0EzD8}x(eEjL`krz$p}Jqi;MNEMBMVS>@3v5){-RR+D^W5}bYwN2#3-$o zs8Na%Lrk2v5fL%BAur{uS=QD45McQGsh4nA<86Xx`)b(F;%B5G_5G8_p9ndhNB6y@N8KN2HJ*9v%$28 zbyWW@9_d5Q@zi3c1O7+LK}O;kt0d`{nqFHt{r3O?y z|AJI(QP;|NhG?Wcb{j&rsLhP|sVaGR*KX(fP||_bOvITAGr%`$f4YEooe{`qD6tFP zn$LaRQKYW7r&qUGg$aXvMXNK>*=+AO-^5nmRt*OsUCAR zC7I{wIf|_)&6n%##@}~qrP;h78OY(|rZ+k!6;5(;3W3DW%%5L5>+Z+G(YFN{X?d!( zKK!P~32*2<_dglgP7zu5A5cqn!|#Oo9x!XBMTr5W`Rr;3D(1nM+a}%=H^)kA&O}&V z{BIe{gMX7T{&yK8R>s5#S*J4U^y_`|Q;?*&1tb0m3cVlG$RK_@-UtO6$QT_nV*0@L2OIysMH-7qO-#1`_et7Uk8Yl%JRg0u(nvCClt-wC z91gYa-<5P7XPxbh439jh`m`-5H&L`Fq8@ZTa!Sd6BtKrk`CI9D-Xu1vGS$Wuas&(na?)W zS~@z;r`aBhcL;gKdA}xeW%9Qzj(_04hL2o0KBb% Date: Mon, 16 Jun 2025 16:39:45 -0400 Subject: [PATCH 9/9] best_model --- week4/trips_modeling.Rmd | 180 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 180 insertions(+) create mode 100644 week4/trips_modeling.Rmd diff --git a/week4/trips_modeling.Rmd b/week4/trips_modeling.Rmd new file mode 100644 index 000000000..b9f7a685a --- /dev/null +++ b/week4/trips_modeling.Rmd @@ -0,0 +1,180 @@ +--- +title: "trips_model" +date: '`r Sys.time()`' +output: + html_document: + #code_folding: hide + number_sections: yes + toc: yes + toc_depth: 3 +--- + +```{r setup, include=FALSE} +library(scales) +library(tidyverse) +library(knitr) + +# set plot theme +theme_set(theme_bw()) +``` + +```{r setup, include=FALSE} +set.seed(42) + + +data <- read_tsv("C:/Users/ds3/Documents/coursework/week4/trips_per_day.tsv") +n <- nrow(data) +shuffled_indices <- sample(n) + + +n_test <- floor(0.1 * n) +test_idx <- shuffled_indices[1:n_test] +test_set <- data[test_idx, ] + + +remaining_idx <- shuffled_indices[(n_test + 1):n] +remaining_data <- data[remaining_idx, ] + + +num_folds <- 5 +trips_per_day <- remaining_data %>% + mutate(fold = ((row_number() - 1) %% num_folds) + 1) + +head(trips_per_day) + +``` + +```{r setup, include=FALSE} + +# fit a model for each polynomial degree +K <- 1:8 +avg_validate_err <- c() +se_validate_err <- c() +for (k in K) { + + # do 5-fold cross-validation within each value of k + validate_err <- c() + for (f in 1:num_folds) { + # fit on the training data + trips_per_day_train <- filter(trips_per_day, fold != f) + model <- lm(num_trips ~ poly(tmin, k, raw = T), data=trips_per_day_train) + + # evaluate on the validation data + trips_per_day_validate <- filter(trips_per_day, fold == f) + validate_err[f] <- sqrt(mean((predict(model, trips_per_day_validate) - trips_per_day_validate$num_trips)^2)) + } + + # compute the average validation error across folds + # and the standard error on this estimate + avg_validate_err[k] <- mean(validate_err) + se_validate_err[k] <- sd(validate_err) / sqrt(num_folds) +} + +``` + +```{r setup, include=FALSE} +# plot the validate error, highlighting the value of k with the lowest average error +plot_data <- data.frame(K, avg_validate_err, se_validate_err) +ggplot(plot_data, aes(x=K, y=avg_validate_err)) + + geom_pointrange(aes(ymin=avg_validate_err - se_validate_err, + ymax=avg_validate_err + se_validate_err, + color=avg_validate_err == min(avg_validate_err))) + + geom_line(color = "red") + + scale_x_continuous(breaks=1:12) + + theme(legend.position="none") + + xlab('Polynomial Degree') + + ylab('RMSE on validation data') +``` + + +```{r setup, include=FALSE} +n <- nrow(trips_per_day) +train_size <- floor(0.8 * n) + +# Randomly shuffle row indices +shuffled_indices <- sample(n) + +# Index split +train_indices <- shuffled_indices[1:train_size] +validation_indices <- shuffled_indices[(train_size + 1):n] + +# Create the two datasets +trips_per_day_train <- trips_per_day[train_indices, ] +trips_per_day_validate <- trips_per_day[validation_indices, ] +``` + + + + +```{r setup, include=FALSE} + +# model with features: num_trips prcp snwd snow tmax tmin +model_with_all_features <- trips_per_day_train |> select(-ymd, -date, -fold) + +model_with_all <- lm(num_trips ~ ., data =model_with_all_features) +model_with_all + +validate_err_with_all <- sqrt(mean((predict(model_with_all, trips_per_day_validate) - trips_per_day_validate$num_trips)^2)) + +# High validation error low predictive power +validate_err_with_all +``` + + +```{r setup, include=FALSE} + +# model with features and interaction: num_trips prcp snwd snow tmax tmin + + +model_with_all_interactions <- lm(num_trips ~ .^2, data = model_with_all_features) +model_with_all_interactions + + +validate_err_with_interactions_all <- sqrt(mean((predict(model_with_all_interactions, trips_per_day_validate) - trips_per_day_validate$num_trips)^2)) + +# High validation error low predictive power +validate_err_with_interactions_all +``` + + +```{r setup, include=FALSE} + +# model with features and interaction with tmin: num_trips prcp snwd snow tmax tmin + + +model_with_all_interactions_with_tmin <- model <- lm( + num_trips ~ prcp + snow + tmax + tmin + + prcp:tmin + snwd:tmin + snow:tmin + tmax:tmin, + data = model_with_all_features +) + +model_with_all_interactions_with_tmin + + +validate_err_with_interactions_with_min_all <- sqrt(mean((predict(model_with_all_interactions_with_tmin, trips_per_day_validate) - trips_per_day_validate$num_trips)^2)) + +# High validation error low predictive power +validate_err_with_interactions_with_min_all + +``` + + +```{r setup, include=FALSE} +# Fit the model on the training data with poly(., 4) and interactions with tmin +model_with_poly4_and_linear_interactions <- lm( + num_trips ~ + + data = model_with_all_features +) + +# Predict on the validation set +predictions <- predict(model_with_poly4_and_linear_interactions, trips_per_day_validate) + +# Compute validation RMSE +validate_err_poly4_with_interactions <- sqrt( + mean((predictions - trips_per_day_validate$num_trips)^2) +) + +# Print validation error +validate_err_poly4_with_interactions +``` \ No newline at end of file