maptest.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c254df4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "from geopy.distance import geodesic\n",
    "\n",
    "# Define a function to fetch coordinates using Overpass API\n",
    "def get_coordinates(city_name):\n",
    "    overpass_url = \"http://overpass-api.de/api/interpreter\"\n",
    "    overpass_query = f\"\"\"\n",
    "    [out:json];\n",
    "    node[\"name\"=\"{city_name}\"][\"place\"=\"city\"];\n",
    "    out body;\n",
    "    \"\"\"\n",
    "    response = requests.get(overpass_url, params={'data': overpass_query})\n",
    "    data = response.json()\n",
    "\n",
    "    # Extract the first result's latitude and longitude\n",
    "    if data['elements']:\n",
    "        lat = data['elements'][0]['lat']\n",
    "        lon = data['elements'][0]['lon']\n",
    "        return lat, lon\n",
    "    else:\n",
    "        raise ValueError(f\"Coordinates for {city_name} not found\")\n",
    "\n",
    "# Fetch coordinates for Place1 and Place2\n",
    "place1_coords = get_coordinates(\"Plymouth\")\n",
    "place2_coords = get_coordinates(\"Swansea\")\n",
    "\n",
    "# Calculate the distance between Place1 and Place2\n",
    "distance = geodesic(place1_coords, place2_coords).miles\n",
    "\n",
    "print(f\"The distance between Place1 and Place2 is approximately {distance:.2f} miles.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1ee7269",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "# Example locations (latitude, longitude)\n",
    "locations = [\n",
    "    (50.3755, -4.1427, 'Plymouth'),\n",
    "    (51.5074, -0.1278, 'London'),\n",
    "    (54.9784, -1.6174, 'Newcastle upon Tyne'),\n",
    "    (53.4084, -2.9916, 'Liverpool'),\n",
    "    (51.6214, -3.9436, 'Swansea')\n",
    "]\n",
    "\n",
    "# Create a new figure\n",
    "plt.figure(figsize=(10, 7))\n",
    "\n",
    "# Create a Basemap instance\n",
    "m = Basemap(projection='merc', llcrnrlat=40, urcrnrlat=60,\n",
    "            llcrnrlon=-10, urcrnrlon=5, resolution='i')\n",
    "\n",
    "# Draw coastlines and countries\n",
    "m.drawcoastlines()\n",
    "m.drawcountries()\n",
    "\n",
    "# Draw map boundaries and parallels/meridians\n",
    "m.drawmapboundary(fill_color='aqua')\n",
    "m.fillcontinents(color='lightgreen', lake_color='aqua')\n",
    "\n",
    "# Plot each location\n",
    "for lat, lon, name in locations:\n",
    "    x, y = m(lon, lat)\n",
    "    m.plot(x, y, 'o', markersize=8, label=name)\n",
    "    plt.text(x, y, ' ' + name, fontsize=12, ha='left', va='bottom')\n",
    "\n",
    "# Add a legend\n",
    "plt.legend(loc='best')\n",
    "\n",
    "# Show the plot\n",
    "plt.title('Locations on 2D Map')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5c16636",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "places=[\"Plymouth\",\"Swansea\"]\n",
    "#places=[\"Plymouth\",\"Swansea\",\"Liverpool\",\"Newcastle upon Tyne\",\"London\"]\n",
    "\n",
    "#def circular_pairs(lst):\n",
    "    #return list(zip(lst, lst[1:] + [lst[0]]))\n",
    "    \n",
    "def all_pairs(lst):\n",
    "    pairs = []\n",
    "    for i in range(len(lst)):\n",
    "        for j in range(i + 1, len(lst)):\n",
    "            pairs.append((lst[i], lst[j]))\n",
    "    return pairs\n",
    "\n",
    "result = all_pairs(places)\n",
    "print(result)\n",
    "\n",
    "#def map_dist(pair):\n",
    "actual_dist=[]\n",
    "for place1, place2 in result:\n",
    "    coord1=get_coordinates(place1)\n",
    "    coord2=get_coordinates(place2)\n",
    "    dist=geodesic(coord1, coord2).miles\n",
    "    actual_dist.append(dist)\n",
    "\n",
    "    print(f\"Distance in miles between {place1} & {place2} is {dist:.2f}\")\n",
    "    \n",
    "    #lats.extend([coord1[0], coord2[0]])\n",
    "    #print(lats)\n",
    "    #lons.extend([coord1[1], coord2[1]])\n",
    "    #print(lons)\n",
    "    #pairs.append((coord1, coord2))\n",
    "    #print(pairs)\n",
    "#print(actual_dist[0])\n",
    "#for k in range(10):\n",
    "    #map_dist=actual_dist[k]\n",
    "    #print(map_dist)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fbcb6b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install folium --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "530c2b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47920988",
   "metadata": {},
   "outputs": [],
   "source": [
    "b=folium.Map(location=[52.158485418573186, -1.9805237617708022], zoom_start=6)\n",
    "#b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2446b8a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "def predicted_distance(lat1, lon1, lat2, lon2):\n",
    "    delta_x = (lon2 - lon1) \n",
    "    delta_y = (lat2 - lat1) \n",
    "    return np.sqrt(delta_x**2 + delta_y**2)\n",
    "\n",
    "def partial_deri(lat1_val, lon1_val, lat2_val, lon2_val):\n",
    "    lat1,lon1,lat2,lon2=sp.symbols('lat1 lon1 lat2 lon2')\n",
    "    delta_x = (lon2 - lon1)\n",
    "    delta_y = (lat2 - lat1)\n",
    "    p = sp.sqrt(delta_x**2 + delta_y**2)\n",
    "    \n",
    "    partial_lat1 = sp.diff(p, lat1)\n",
    "    partial_lon1 = sp.diff(p, lon1)\n",
    "    partial_lat2 = sp.diff(p, lat2)\n",
    "    partial_lon2 = sp.diff(p, lon2)\n",
    "    \n",
    "    values = {lat1: lat1_val, lon1: lon1_val, lat2: lat2_val, lon2: lon2_val}\n",
    "\n",
    "    # Substitute the numerical values and compute the partial derivatives\n",
    "    partial_lat1_num = partial_lat1.subs(values).evalf()\n",
    "    partial_lon1_num = partial_lon1.subs(values).evalf()\n",
    "    partial_lat2_num = partial_lat2.subs(values).evalf()\n",
    "    partial_lon2_num = partial_lon2.subs(values).evalf()\n",
    "    \n",
    "    # Return the numerical values of the partial derivatives\n",
    "    return partial_lat1_num, partial_lon1_num, partial_lat2_num, partial_lon2_num\n",
    "    \n",
    "    \n",
    "# Define the cost function\n",
    "def cost_func(locations, coord1, coord2):\n",
    "    m = len(locations)\n",
    "    total_cost = 0\n",
    "    for i in range(m):\n",
    "        for j in range(i + 1, m):\n",
    "            #coord1=tuple(j*np.array(coord1))\n",
    "            #coord2=tuple(j*np.array(coord2))\n",
    "            #print(coord1)\n",
    "            #print(coord2)\n",
    "            #z=map(tuple,np.zeroes((2,),dtype='np.float64'))\n",
    "            #print(z)\n",
    "            \n",
    "            pred_dist = predicted_distance(coord1[0], coord1[1], coord2[0], coord2[1])\n",
    "            partial_deri(coord1[0], coord1[1], coord2[0], coord2[1])\n",
    "            #print(pred_dist)\n",
    "            \n",
    "            coordi1 = get_coordinates(locations[i])\n",
    "            coordi2 = get_coordinates(locations[j])\n",
    "            map_dist = geodesic(coordi1, coordi2).miles\n",
    "            #print(map_dist)\n",
    "            \n",
    "            error = pred_dist - map_dist\n",
    "            #print(error)\n",
    "            \n",
    "            total_cost += error**2\n",
    "            \n",
    "            \n",
    "    return total_cost / (2 * m)\n",
    "\n",
    "#cost=cost_func(places)\n",
    "#print(cost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f85cdf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "\n",
    "def gradient_descent(locations, alpha, num_iterations):\n",
    "    m = len(locations)\n",
    "    coord1=(0.01,0.02)\n",
    "    coord2=(0.03,0.01)\n",
    "    for iteration in range(num_iterations):\n",
    "        # Initialize gradients\n",
    "        d_locations = np.zeros((m, 2))\n",
    "        #print(d_locations)\n",
    "        \n",
    "        for i in range(m):\n",
    "            for j in range(i + 1, m):\n",
    "                coord1=tuple(j*np.array(coord1))\n",
    "                #coord2=tuple(j*np.array(coord2))\n",
    "                \n",
    "                #print(f\"Iteration {iteration+1}, Pair ({i},{j}): Coord1: {coord1}, Coord2: {coord2}\")\n",
    "                \n",
    "                # Check if coord1 and coord2 have two elements\n",
    "                if len(coord1) != 2 or len(coord2) != 2:\n",
    "                    raise ValueError(f\"Invalid coordinates: {coord1}, {coord2}\")\n",
    "                \n",
    "                pred_dist = predicted_distance(coord1[0], coord1[1], coord2[0], coord2[1])\n",
    "                #print(pred_dist)\n",
    "                coordi1 = get_coordinates(locations[i])\n",
    "                coordi2 = get_coordinates(locations[j])\n",
    "                map_dist = geodesic(coordi1, coordi2).miles\n",
    "                #print(map_dist)\n",
    "                error = (pred_dist - map_dist)**2\n",
    "                \n",
    "                \n",
    "            # Optionally print the cost at every iteration\n",
    "            if (iteration + 1) % 1 == 0:\n",
    "                print(coord1,coord2)\n",
    "                cost = cost_func(locations,coord1,coord2)\n",
    "                print(f\"Iteration {iteration+1}: Cost = {cost:.4f}, Location1_lat={d_locations[i][0]}, Location1_lon = {d_locations[i][1]}, Location2_lat={d_locations[j][0]}, Location2_lon = {d_locations[j][1]}\")\n",
    "        #cost = cost_func(locations, theta1, theta2)\n",
    "        #print(f\"Iteration {iteration+1}: Cost = {cost:.4f}, Theta1 = {theta1:.4f}, Theta2 = {theta2:.4f}\")\n",
    "                \n",
    "            # Calculate the partial derivatives\n",
    "                d_locations[i][0],d_locations[i][1],d_locations[j][0],d_locations[j][1] = partial_deri(coord1[0], coord1[1], coord2[0], coord2[1])\n",
    "                \n",
    "                d_locations[i][0] += 2.0*np.sqrt(error)*d_locations[i][0]\n",
    "                d_locations[i][1] += 2.0*np.sqrt(error)*d_locations[i][1]\n",
    "                d_locations[j][0] += 2.0*np.sqrt(error)*d_locations[j][0]\n",
    "                d_locations[j][1] += 2.0*np.sqrt(error)*d_locations[j][1]\n",
    "        \n",
    "        d_locations[i][0] -= alpha * d_locations[i][0] / m\n",
    "        d_locations[i][1] -= alpha * d_locations[i][1] / m\n",
    "        d_locations[j][0] -= alpha * d_locations[j][0] / m\n",
    "        d_locations[j][1] -= alpha * d_locations[j][1] / m\n",
    "        coord1=(d_locations[i][0],d_locations[i][1])\n",
    "        coord2=(d_locations[j][0],d_locations[j][1])\n",
    "        \n",
    "    return locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfb93306",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example usage\n",
    "alpha = 0.01  # Learning rate\n",
    "num_iterations = 3  # Number of iterations\n",
    "#theta1, theta2 = 69.0, 69.0  # Initial values of theta1 and theta2\n",
    "\n",
    "# Assuming 'places' is a list of coordinates and 'get_coordinates' returns the latitude and longitude.\n",
    "updated_locations = gradient_descent(places, alpha, num_iterations)\n",
    "\n",
    "print(\"Updated locations:\", updated_locations)\n",
    "#print(f\"Updated Theta1: {updated_theta1:.4f}, Updated Theta2: {updated_theta2:.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b30f58ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import numpy as np\n",
    "from geopy.distance import geodesic\n",
    "\n",
    "# Define the symbolic variables\n",
    "lat1, lon1, lat2, lon2, map_distance = sp.symbols('lat1 lon1 lat2 lon2 map_distance')\n",
    "\n",
    "# Define the predicted distance function using Euclidean formula\n",
    "predicted_distance = sp.sqrt((lon2 - lon1)**2 + (lat2 - lat1)**2)\n",
    "\n",
    "# Define the cost function (squared error)\n",
    "cost_function = (predicted_distance - map_distance)**2\n",
    "\n",
    "# Compute partial derivatives of the cost function with respect to each variable\n",
    "partial_lat1 = sp.diff(cost_function, lat1)\n",
    "partial_lon1 = sp.diff(cost_function, lon1)\n",
    "partial_lat2 = sp.diff(cost_function, lat2)\n",
    "partial_lon2 = sp.diff(cost_function, lon2)\n",
    "\n",
    "# Convert partial derivatives to lambda functions for efficient numerical computation\n",
    "partial_lat1_fn = sp.lambdify((lat1, lon1, lat2, lon2, map_distance), partial_lat1)\n",
    "partial_lon1_fn = sp.lambdify((lat1, lon1, lat2, lon2, map_distance), partial_lon1)\n",
    "partial_lat2_fn = sp.lambdify((lat1, lon1, lat2, lon2, map_distance), partial_lat2)\n",
    "partial_lon2_fn = sp.lambdify((lat1, lon1, lat2, lon2, map_distance), partial_lon2)\n",
    "\n",
    "# Initialize values\n",
    "lat1_val, lon1_val, lat2_val, lon2_val = 34.0522, -118.2437, 36.7783, -119.4179  # Example: Los Angeles to Fresno\n",
    "learning_rate = 0.01\n",
    "num_iterations = 100\n",
    "\n",
    "# Calculate map distance using geopy\n",
    "map_distance_val = geodesic((lat1_val, lon1_val), (lat2_val, lon2_val)).km\n",
    "\n",
    "# Gradient Descent Algorithm\n",
    "for i in range(num_iterations):\n",
    "    # Compute partial derivatives\n",
    "    grad_lat1 = partial_lat1_fn(lat1_val, lon1_val, lat2_val, lon2_val, map_distance_val)\n",
    "    grad_lon1 = partial_lon1_fn(lat1_val, lon1_val, lat2_val, lon2_val, map_distance_val)\n",
    "    grad_lat2 = partial_lat2_fn(lat1_val, lon1_val, lat2_val, lon2_val, map_distance_val)\n",
    "    grad_lon2 = partial_lon2_fn(lat1_val, lon1_val, lat2_val, lon2_val, map_distance_val)\n",
    "\n",
    "    # Update the variables\n",
    "    lat1_val -= learning_rate * grad_lat1\n",
    "    lon1_val -= learning_rate * grad_lon1\n",
    "    lat2_val -= learning_rate * grad_lat2\n",
    "    lon2_val -= learning_rate * grad_lon2\n",
    "\n",
    "    # Compute the cost to monitor the optimization process\n",
    "    predicted_dist = np.sqrt((lon2_val - lon1_val)**2 + (lat2_val - lat1_val)**2)\n",
    "    cost = (predicted_dist - map_distance_val)**2\n",
    "    print(f\"Iteration {i+1}: Cost = {cost}, lat1 = {lat1_val}, lon1 = {lon1_val}, lat2 = {lat2_val}, lon2 = {lon2_val}\")\n",
    "\n",
    "print(\"Optimization complete.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89fd1201",
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install cartopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24f8d2da",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from geopy.distance import geodesic\n",
    "from sympy import symbols, diff, sqrt, lambdify\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "# Define the predicted distance function using SymPy\n",
    "def predicted_distance_sym(lat1, lon1, lat2, lon2):\n",
    "    delta_x = lon2 - lon1\n",
    "    delta_y = lat2 - lat1\n",
    "    return sqrt(delta_x**2 + delta_y**2)\n",
    "\n",
    "# Actual geographical coordinates of two places\n",
    "actual_lat1, actual_lon1 = 52.5200, 13.4050   # Example: Berlin\n",
    "actual_lat2, actual_lon2 = 48.8566, 2.3522    # Example: Paris\n",
    "\n",
    "# Compute the actual distance between two points using geopy\n",
    "actual_distance = geodesic((actual_lat1, actual_lon1), (actual_lat2, actual_lon2)).kilometers\n",
    "\n",
    "# Initialize the latitudes and longitudes\n",
    "lat1, lon1, lat2, lon2 = 0.01, 0.02, 0.03, 0.04\n",
    "\n",
    "# Define symbols for sympy\n",
    "lat1_sym, lon1_sym, lat2_sym, lon2_sym = symbols('lat1_sym lon1_sym lat2_sym lon2_sym')\n",
    "\n",
    "# Predicted distance with sympy symbols\n",
    "predicted_dist = predicted_distance_sym(lat1_sym, lon1_sym, lat2_sym, lon2_sym)\n",
    "\n",
    "# Cost function as the squared error\n",
    "cost_function = (predicted_dist - actual_distance) ** 2\n",
    "\n",
    "# Compute the partial derivatives of the cost function\n",
    "grad_lat1 = diff(cost_function, lat1_sym)\n",
    "grad_lon1 = diff(cost_function, lon1_sym)\n",
    "grad_lat2 = diff(cost_function, lat2_sym)\n",
    "grad_lon2 = diff(cost_function, lon2_sym)\n",
    "\n",
    "# Convert the gradients to numerical functions using lambdify\n",
    "grad_lat1_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym), grad_lat1, 'numpy')\n",
    "grad_lon1_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym), grad_lon1, 'numpy')\n",
    "grad_lat2_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym), grad_lat2, 'numpy')\n",
    "grad_lon2_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym), grad_lon2, 'numpy')\n",
    "\n",
    "# Gradient Descent parameters\n",
    "learning_rate = 0.0001  # Reduced learning rate to prevent overshooting\n",
    "num_iterations = 10  # Define the number of iterations\n",
    "\n",
    "# Set up the map\n",
    "plt.figure(figsize=(12, 6))\n",
    "map = Basemap(projection='merc', llcrnrlat=-80, urcrnrlat=80, llcrnrlon=-180, urcrnrlon=180, resolution='c')\n",
    "map.drawcoastlines()\n",
    "map.drawcountries()\n",
    "map.drawmapboundary(fill_color='aqua')\n",
    "map.fillcontinents(color='lightgreen', lake_color='aqua')\n",
    "\n",
    "# Plot initial actual locations\n",
    "x_actual1, y_actual1 = map(actual_lon1, actual_lat1)\n",
    "x_actual2, y_actual2 = map(actual_lon2, actual_lat2)\n",
    "map.plot(x_actual1, y_actual1, 'ro', markersize=8, label='Actual Location 1')\n",
    "map.plot(x_actual2, y_actual2, 'bo', markersize=8, label='Actual Location 2')\n",
    "\n",
    "# Initialize plot points for the optimization locations\n",
    "x_lat1, y_lat1 = map(lon1, lat1)\n",
    "x_lat2, y_lat2 = map(lon2, lat2)\n",
    "opt_loc1, = map.plot(x_lat1, y_lat1, 'r^', markersize=8, label='Optimized Location 1')\n",
    "opt_loc2, = map.plot(x_lat2, y_lat2, 'b^', markersize=8, label='Optimized Location 2')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "# Function to update the plot points\n",
    "def update_plot(lat1, lon1, lat2, lon2):\n",
    "    x_lat1, y_lat1 = map(lon1, lat1)\n",
    "    x_lat2, y_lat2 = map(lon2, lat2)\n",
    "    opt_loc1.set_data(x_lat1, y_lat1)\n",
    "    opt_loc2.set_data(x_lat2, y_lat2)\n",
    "    plt.draw()\n",
    "    plt.pause(0.01)  # Adjust pause duration if needed\n",
    "\n",
    "# Perform Gradient Descent optimization with live plot updates\n",
    "for i in range(num_iterations):\n",
    "    # Evaluate gradients numerically\n",
    "    grad_lat1_val = float(grad_lat1_func(lat1, lon1, lat2, lon2))\n",
    "    grad_lon1_val = float(grad_lon1_func(lat1, lon1, lat2, lon2))\n",
    "    grad_lat2_val = float(grad_lat2_func(lat1, lon1, lat2, lon2))\n",
    "    grad_lon2_val = float(grad_lon2_func(lat1, lon1, lat2, lon2))\n",
    "\n",
    "    # Update the variables using the gradients\n",
    "    lat1 -= learning_rate * grad_lat1_val\n",
    "    lon1 -= learning_rate * grad_lon1_val\n",
    "    lat2 -= learning_rate * grad_lat2_val\n",
    "    lon2 -= learning_rate * grad_lon2_val\n",
    "\n",
    "    # Update the plot with new positions\n",
    "    update_plot(lat1, lon1, lat2, lon2)\n",
    "\n",
    "plt.title('Live Optimization of Lat/Lon Coordinates')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "acb4fa79",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "879.6993161290122\n",
      "Location 1: (-377.584830, -458.555659)\n",
      "Location 2: (40.109217, 315.655275)\n",
      "Location 3: (337.565613, 143.020384)\n",
      "879.699316129010\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvxklEQVR4nO3df3BU5aH/8c9uQjYgJoCRhMRgQAG1SIggMSJVajRGJr3etldGucIXfw0KDppqJf4g5bY1apVLqyhXK1CnV0EdoVYowg0CRVMpgVSoAiJgUiABimSTgBvIPt8/SA6sJJiFbB7Ceb9mdpacfc45zz5mzGeeX8djjDECAACwxGu7AgAAwN0IIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMCqDhVGVq1apby8PCUnJ8vj8WjhwoVhX8MYo+eee079+/eXz+dTSkqKfvWrX7V9ZQEAQKtE265AOOrq6pSenq4777xTP/rRj07pGpMnT9bSpUv13HPP6fLLL9f+/fu1f//+Nq4pAABoLU9HfVCex+PRggULdMsttzjHAoGAHn/8cb355ps6cOCABg4cqGeeeUbXXXedJOnzzz/XoEGDtHHjRg0YMMBOxQEAQIgONUzzXSZNmqSSkhLNmzdPn376qf7jP/5DN910k7744gtJ0p/+9Cf17dtX77//vvr06aO0tDTdfffd9IwAAGDRWRNGysvLNWfOHL399tsaMWKELrroIj388MO65pprNGfOHEnStm3b9NVXX+ntt9/W66+/rrlz56q0tFQ/+clPLNceAAD36lBzRk5mw4YNamhoUP/+/UOOBwIBnXfeeZKkYDCoQCCg119/3Sn32muvaciQIdq8eTNDNwAAWHDWhJHa2lpFRUWptLRUUVFRIZ917dpVktSrVy9FR0eHBJZLL71U0tGeFcIIAADt76wJIxkZGWpoaNCePXs0YsSIZssMHz5cR44c0ZdffqmLLrpIkrRlyxZJ0oUXXthudQUAAMd0qNU0tbW12rp1q6Sj4WP69OkaOXKkevTood69e+s///M/9dFHH+n5559XRkaG9u7dq+LiYg0aNEijRo1SMBjUlVdeqa5du2rGjBkKBoOaOHGi4uLitHTpUsvfDgAAd+pQYWTFihUaOXLkCcfHjRunuXPn6vDhw/rlL3+p119/XTt37lRCQoKuuuoqTZs2TZdffrkkadeuXXrggQe0dOlSnXPOOcrNzdXzzz+vHj16tPfXAQAA6mBhBAAAnH3OmqW9AACgYyKMAAAAqzrEappgMKhdu3bp3HPPlcfjsV0dAADQCsYY1dTUKDk5WV5vy/0fHSKM7Nq1S6mpqbarAQAATkFFRYUuuOCCFj/vEGHk3HPPlXT0y8TFxVmuDQAAaA2/36/U1FTn73hLOkQYaRqaiYuLI4wAANDBfNcUCyawAgAAqwgjAADAqrDDyKpVq5SXl6fk5GR5PB4tXLjwO88JBAJ6/PHHdeGFF8rn8yktLU2zZ88+lfoCAICzTNhzRurq6pSenq4777xTP/rRj1p1zq233qqqqiq99tpruvjii7V7924Fg8GwKwsAAM4+YYeR3Nxc5ebmtrr8kiVLtHLlSm3bts15/ktaWlq4twUAAGepiM8Zee+99zR06FA9++yzSklJUf/+/fXwww/r0KFDLZ4TCATk9/tDXgAA4OwU8aW927Zt0+rVqxUbG6sFCxZo3759uv/++/Wvf/1Lc+bMafacoqIiTZs2LdJVAwAAZ4CI94wEg0F5PB797//+r4YNG6abb75Z06dP1+9///sWe0cKCgpUXV3tvCoqKiJdTQAAYEnEe0Z69eqllJQUxcfHO8cuvfRSGWP0z3/+U/369TvhHJ/PJ5/PF+mqAQCAM0DEe0aGDx+uXbt2qba21jm2ZcsWeb3ek+5TDwAA3CHsMFJbW6uysjKVlZVJkrZv366ysjKVl5dLOjrEMnbsWKf87bffrvPOO0/jx4/XZ599plWrVumRRx7RnXfeqc6dO7fNtwAAAB1W2GFk7dq1ysjIUEZGhiQpPz9fGRkZmjp1qiRp9+7dTjCRpK5du2rZsmU6cOCAhg4dqjFjxigvL0+//e1v2+grAACAjsxjjDG2K/Fd/H6/4uPjVV1d3aYPynun9J/auLNaNw1M0lV9z2uz6wIAgNb//Xb1s2lWbtmruR/v0Ge72McEAABbXB1GTv5AYwAA0B5cHUaanPHjVAAAnMVcHUY8jV0jHWDaDAAAZy13hxHbFQAAAC4PI41dI3SMAABgj6vDCAAAsM/VYaRpmMYwhRUAAGtcHUbkTGC1Ww0AANzM1WHE05hGyCIAANjj6jACAADsc3UY8TBMAwCAde4OI43vTGAFAMAed4cRekYAALDO3WGEPVgBALDO1WEEAADY5+owwoPyAACwjzAi5owAAGCTq8MIz+0FAMA+l4eRo+gYAQDAHleHEYZpAACwz91hpPGdTc8AALDH3WGEKSMAAFjn6jDShGEaAADscXUYadqBlSwCAIA97g4jzqQR4ggAALa4O4w0vhNFAACwx9VhBAAA2OfqMOJpHKdhlAYAAHtcHUaasM8IAAD2uDqMsAMrAAD2uTqMAAAA+8IOI6tWrVJeXp6Sk5Pl8Xi0cOHCVp/70UcfKTo6WoMHDw73thHBPiMAANgXdhipq6tTenq6Zs6cGdZ5Bw4c0NixY3X99deHe8uIYZgGAAD7osM9ITc3V7m5uWHfaMKECbr99tsVFRUVVm9KJPGgPAAA7GuXOSNz5szRtm3bVFhY2KrygUBAfr8/5AUAAM5OEQ8jX3zxhaZMmaI//OEPio5uXUdMUVGR4uPjnVdqampE6uZhC1YAAKyLaBhpaGjQ7bffrmnTpql///6tPq+goEDV1dXOq6KiIiL1czY9i8jVAQBAa4Q9ZyQcNTU1Wrt2rdavX69JkyZJkoLBoIwxio6O1tKlS/WDH/zghPN8Pp98Pl8kqybpuDkjzGAFAMCaiIaRuLg4bdiwIeTYSy+9pOXLl+udd95Rnz59Inl7AADQAYQdRmpra7V161bn5+3bt6usrEw9evRQ7969VVBQoJ07d+r111+X1+vVwIEDQ87v2bOnYmNjTzhuBUt7AQCwLuwwsnbtWo0cOdL5OT8/X5I0btw4zZ07V7t371Z5eXnb1TCC2PQMAAD7PKYDTJjw+/2Kj49XdXW14uLi2uy6zyzZpJdXfKk7h/fR1LzL2uy6AACg9X+/eTaN2PQMAACbXB1GPN9dBAAARJi7wwgTWAEAsM7dYYS+EQAArHN1GGnSAebwAgBw1nJ1GPHQMQIAgHXuDiON7/SLAABgj6vDSFPXCKM0AADY4+owcqxnhDQCAIAtrg4jAADAPleHEfYZAQDAPneHER6UBwCAde4OI/SMAABgnavDCAAAsM/VYeTYnmd0jQAAYIu7wwjDNAAAWOfyMMKmZwAA2ObqMAIAAOwjjIgdWAEAsMnVYYQ5IwAA2OfuMMKmZwAAWOfqMNKEnhEAAOxxdRjxeL67DAAAiCx3h5HGdyawAgBgj7vDyLE0AgAALHF1GGlCFgEAwB5XhxGPmDQCAIBt7g4jzj4j9I0AAGCLq8NIE6IIAAD2EEbEPiMAANjk6jDiYaMRAACsc3cYaXynYwQAAHvCDiOrVq1SXl6ekpOT5fF4tHDhwpOWf/fdd3XDDTfo/PPPV1xcnLKysvTBBx+can3bFBNYAQCwL+wwUldXp/T0dM2cObNV5VetWqUbbrhBixcvVmlpqUaOHKm8vDytX78+7MpGClEEAAB7osM9ITc3V7m5ua0uP2PGjJCfn3rqKf3xj3/Un/70J2VkZIR7+zbFjBEAAOwLO4ycrmAwqJqaGvXo0aPFMoFAQIFAwPnZ7/dHpC7OBFa6RgAAsKbdJ7A+99xzqq2t1a233tpimaKiIsXHxzuv1NTUiNTlWBYhjQAAYEu7hpE33nhD06ZN01tvvaWePXu2WK6goEDV1dXOq6KiIqL1Yv4qAAD2tNswzbx583T33Xfr7bffVnZ29knL+nw++Xy+iNfJWdpLGAEAwJp26Rl58803NX78eL355psaNWpUe9yyddj0DAAA68LuGamtrdXWrVudn7dv366ysjL16NFDvXv3VkFBgXbu3KnXX39d0tGhmXHjxuk3v/mNMjMzVVlZKUnq3Lmz4uPj2+hrnJpjm57RNQIAgC1h94ysXbtWGRkZzrLc/Px8ZWRkaOrUqZKk3bt3q7y83Cn/yiuv6MiRI5o4caJ69erlvCZPntxGX+H0MUwDAIA9YfeMXHfddSfdsXTu3LkhP69YsSLcW7QbVvYCAGCfy59Nw5wRAABsc3cYcZ5NY7ceAAC4mbvDiPMv0ggAALa4Oow0oWcEAAB7XB1G2GYEAAD73B1GGgdq6BgBAMAeV4cRORNYiSMAANji7jDSiCgCAIA9rg4jTBkBAMA+d4eRxhmsjNIAAGCPu8NI4ztZBAAAe1wdRpowgRUAAHtcHUbYZwQAAPsIIwAAwCp3hxExgRUAANtcHUaaGKawAgBgjavDiMfZgdVuPQAAcDNXhxEAAGCfq8MIm54BAGCfq8NIE+aMAABgj6vDiLMDK1kEAABr3B1G2GcEAADr3B1GmvYZsVwPAADczNVhxEEaAQDAGleHEWefEdIIAADWuDuM2K4AAABweRhhB1YAAKxzdRhpQhYBAMAel4eRph1YiSMAANji6jDCPiMAANjn7jDS+E6/CAAA9rg6jDRhlAYAAHvCDiOrVq1SXl6ekpOT5fF4tHDhwu88Z8WKFbriiivk8/l08cUXa+7cuadQ1bbnPLXXcj0AAHCzsMNIXV2d0tPTNXPmzFaV3759u0aNGqWRI0eqrKxMDz74oO6++2598MEHYVe2rTlTRugaAQDAmuhwT8jNzVVubm6ry8+aNUt9+vTR888/L0m69NJLtXr1av33f/+3cnJywr19m2ICKwAA9kV8zkhJSYmys7NDjuXk5KikpKTFcwKBgPx+f8grkugXAQDAnoiHkcrKSiUmJoYcS0xMlN/v16FDh5o9p6ioSPHx8c4rNTU1InVjB1YAAOw7I1fTFBQUqLq62nlVVFRE5D6epk3P6BsBAMCasOeMhCspKUlVVVUhx6qqqhQXF6fOnTs3e47P55PP54t01XhSHgAAZ4CI94xkZWWpuLg45NiyZcuUlZUV6Vt/J2fTMzpGAACwJuwwUltbq7KyMpWVlUk6unS3rKxM5eXlko4OsYwdO9YpP2HCBG3btk0/+9nPtGnTJr300kt666239NBDD7XNN2gDhBEAAOwJO4ysXbtWGRkZysjIkCTl5+crIyNDU6dOlSTt3r3bCSaS1KdPHy1atEjLli1Tenq6nn/+ef3ud7+zvqxXYtMzAADOBGHPGbnuuutO+pTb5nZXve6667R+/fpwbxVxTBkBAMC+M3I1TXs5trSXvhEAAGxxdRgBAAD2uTqMOPuM0DECAIA17g4jTBoBAMA6d4eRxnd2YAUAwB5Xh5EmDNMAAGCPu8NI02oau7UAAMDVXB1GPOw0AgCAde4OI+wzAgCAda4OI02IIgAA2OPqMOIM0pBGAACwxt1hhAflAQBgncvDiO0aAAAAV4eRJkxgBQDAHleHkWM7sAIAAFvcHUacpb126wEAgJu5OoyITc8AALDO5WHkKB6UBwCAPa4OIwzTAABgn7vDSOM7YQQAAHvcHUbYaAQAAOtcHUYAAIB9rg4jx4ZpGKcBAMAWd4eRpgmsdqsBAICruTuMNPaN0DECAIA97g4jTs8IaQQAAFsII5KCZBEAAKxxdRjxehimAQDANleHkWPbjJBGAACwxd1hpHECK8M0AADY4+ow4nWeTUMaAQDAllMKIzNnzlRaWppiY2OVmZmpNWvWnLT8jBkzNGDAAHXu3Fmpqal66KGH9M0335xShdsSE1gBALAv7DAyf/585efnq7CwUOvWrVN6erpycnK0Z8+eZsu/8cYbmjJligoLC/X555/rtdde0/z58/XYY4+dduVPl8eZwEoaAQDAlrDDyPTp03XPPfdo/PjxuuyyyzRr1ix16dJFs2fPbrb8xx9/rOHDh+v2229XWlqabrzxRt12223f2ZvSHpzt4K3WAgAAdwsrjNTX16u0tFTZ2dnHLuD1Kjs7WyUlJc2ec/XVV6u0tNQJH9u2bdPixYt18803n0a124aHpb0AAFgXHU7hffv2qaGhQYmJiSHHExMTtWnTpmbPuf3227Vv3z5dc801MsboyJEjmjBhwkmHaQKBgAKBgPOz3+8Pp5qtxgRWAADsi/hqmhUrVuipp57SSy+9pHXr1undd9/VokWL9Itf/KLFc4qKihQfH++8UlNTI1I3lvYCAGBfWD0jCQkJioqKUlVVVcjxqqoqJSUlNXvOk08+qTvuuEN33323JOnyyy9XXV2d7r33Xj3++OPyek/MQwUFBcrPz3d+9vv9EQkkPJsGAAD7wuoZiYmJ0ZAhQ1RcXOwcCwaDKi4uVlZWVrPnHDx48ITAERUVJanl4RGfz6e4uLiQVySwtBcAAPvC6hmRpPz8fI0bN05Dhw7VsGHDNGPGDNXV1Wn8+PGSpLFjxyolJUVFRUWSpLy8PE2fPl0ZGRnKzMzU1q1b9eSTTyovL88JJbZ4jnWNAAAAS8IOI6NHj9bevXs1depUVVZWavDgwVqyZIkzqbW8vDykJ+SJJ56Qx+PRE088oZ07d+r8889XXl6efvWrX7XdtzhFXoZpAACwzmM6wFISv9+v+Ph4VVdXt+mQTWX1N7qqqFhRXo++fMr+UmMAAM4mrf37zbNpxNJeAABscnUYERNYAQCwztVhxNs0gRUAAFjj6jByfBRhqAYAADvcHUaO6xlhqAYAADtcHUa8x3WN0DMCAIAdrg4jHtEzAgCAbe4OI8d9ezY+AwDADneHkeP+zSgNAAB2uDuMHDeBlTACAIAdrg4jIRNYGaYBAMAKV4cRJrACAGCfu8MIS3sBALCOMNKIKAIAgB2uDiPHP5vGBC1WBAAAF3N1GAlZ2kvfCAAAVrg7jPBsGgAArHN1GOHZNAAA2OfqMBKy6ZnFegAA4GauDiPSsRU1QXpGAACwgjDS9A+yCAAAVhBGGrtGmMAKAIAdrg8jTZNYWdoLAIAdrg8jTc+nYcoIAAB2EEaYwAoAgFWEkaZhGrIIAABWuD6MND2fhjACAIAdrg8jTUt7mcAKAIAdhBF6RgAAsIowwgRWAACsIow0vhNFAACww/VhxOttGqYhjgAAYMMphZGZM2cqLS1NsbGxyszM1Jo1a05a/sCBA5o4caJ69eoln8+n/v37a/HixadU4bbm9IyQRQAAsCI63BPmz5+v/Px8zZo1S5mZmZoxY4ZycnK0efNm9ezZ84Ty9fX1uuGGG9SzZ0+98847SklJ0VdffaVu3bq1Rf1PmzOB1XI9AABwq7DDyPTp03XPPfdo/PjxkqRZs2Zp0aJFmj17tqZMmXJC+dmzZ2v//v36+OOP1alTJ0lSWlra6dW6DXmZwAoAgFVhDdPU19ertLRU2dnZxy7g9So7O1slJSXNnvPee+8pKytLEydOVGJiogYOHKinnnpKDQ0Np1fzNsPSXgAAbAqrZ2Tfvn1qaGhQYmJiyPHExERt2rSp2XO2bdum5cuXa8yYMVq8eLG2bt2q+++/X4cPH1ZhYWGz5wQCAQUCAednv98fTjXDQs8IAAB2RXw1TTAYVM+ePfXKK69oyJAhGj16tB5//HHNmjWrxXOKiooUHx/vvFJTUyNWP55NAwCAXWGFkYSEBEVFRamqqirkeFVVlZKSkpo9p1evXurfv7+ioqKcY5deeqkqKytVX1/f7DkFBQWqrq52XhUVFeFUMyweZz0NAACwIawwEhMToyFDhqi4uNg5FgwGVVxcrKysrGbPGT58uLZu3apgMOgc27Jli3r16qWYmJhmz/H5fIqLiwt5RQrDNAAA2BX2ME1+fr5effVV/f73v9fnn3+u++67T3V1dc7qmrFjx6qgoMApf99992n//v2aPHmytmzZokWLFumpp57SxIkT2+5bnAaeTQMAgF1hL+0dPXq09u7dq6lTp6qyslKDBw/WkiVLnEmt5eXl8nqPZZzU1FR98MEHeuihhzRo0CClpKRo8uTJevTRR9vuW5wGnk0DAIBdHtMB9kH3+/2Kj49XdXV1mw/ZjHh2uSr2H9K791+tK3p3b9NrAwDgZq39+82zaRimAQDAKteHkWPPpiGNAABgA2GEZ9MAAGAVYaRpAmuQOAIAgA2EkcZ3oggAAHa4PowwgRUAALtcH0aOPZuGNAIAgA2EETGBFQAAmwgj7MAKAIBVhBHmjAAAYJXrw0jTU3vJIgAA2OH6MMIwDQAAdrk+jHg9dI0AAGCT68NI06Zn9IwAAGCH68OImMAKAIBVrg8jTGAFAMAu14cRhmkAALDL9WGEZ9MAAGCX68MIz6YBAMAuwgjPpgEAwCrCCJueAQBglevDSNOckSBZBAAAK1wfRqK8TRNYSSMAANjg+jDSNEzTQNcIAABWuD6MNPWMEEYAALCDMMI+IwAAWOX6MOJpDCMNpBEAAKxwfRiJamwBhmkAALDD9WHk2HbwhBEAAGwgjDCBFQAAq1wfRqLY9AwAAKtcH0a8bAcPAIBVpxRGZs6cqbS0NMXGxiozM1Nr1qxp1Xnz5s2Tx+PRLbfcciq3jQiGaQAAsCvsMDJ//nzl5+ersLBQ69atU3p6unJycrRnz56Tnrdjxw49/PDDGjFixClXNhJ4Ng0AAHaFHUamT5+ue+65R+PHj9dll12mWbNmqUuXLpo9e3aL5zQ0NGjMmDGaNm2a+vbte1oVbmvH5oyQRgAAsCGsMFJfX6/S0lJlZ2cfu4DXq+zsbJWUlLR43n/913+pZ8+euuuuu1p1n0AgIL/fH/KKFIZpAACwK6wwsm/fPjU0NCgxMTHkeGJioiorK5s9Z/Xq1Xrttdf06quvtvo+RUVFio+Pd16pqanhVDMsTGAFAMCuiK6mqamp0R133KFXX31VCQkJrT6voKBA1dXVzquioiJidWx6UF6QnhEAAKyIDqdwQkKCoqKiVFVVFXK8qqpKSUlJJ5T/8ssvtWPHDuXl5TnHgsHg0RtHR2vz5s266KKLTjjP5/PJ5/OFU7VTxgRWAADsCqtnJCYmRkOGDFFxcbFzLBgMqri4WFlZWSeUv+SSS7RhwwaVlZU5rx/+8IcaOXKkysrKIjr80lpeHpQHAIBVYfWMSFJ+fr7GjRunoUOHatiwYZoxY4bq6uo0fvx4SdLYsWOVkpKioqIixcbGauDAgSHnd+vWTZJOOG5L04PyGKYBAMCOsMPI6NGjtXfvXk2dOlWVlZUaPHiwlixZ4kxqLS8vl9fbcTZ29bK0FwAAq8IOI5I0adIkTZo0qdnPVqxYcdJz586deyq3jJhjS3stVwQAAJfqOF0YEcLSXgAA7HJ9GGEHVgAA7HJ9GGEHVgAA7CKMsM8IAABWuT6MsAMrAAB2uT6MNHaMsOkZAACWuD6MMIEVAAC7CCMM0wAAYJXrw4jHeTaN5YoAAOBSrg8jUWx6BgCAVa4PI16GaQAAsIowwgRWAACscn0YieJBeQAAWOX6MMKD8gAAsIswwjANAABWEUY8PCgPAACbXB9GnE3P6BkBAMAK14eRY0t7LVcEAACXIozwoDwAAKxyfRhpelCeIYwAAGCF68OIhwmsAABY5fow4mx6RhYBAMAK14eRaGcHVmawAgBgA2Gk8bG9R+gaAQDACsKI92gTHGHOCAAAVhBGnJ4RhmkAALCBMNI4Z+QwwzQAAFjh+jDSKappmIaeEQAAbHB9GGECKwAAdhFGGiewHmbOCAAAVrg+jHSKYgdWAABscn0YiW6cM3KYMAIAgBWnFEZmzpyptLQ0xcbGKjMzU2vWrGmx7KuvvqoRI0aoe/fu6t69u7Kzs09avr118rK0FwAAm8IOI/Pnz1d+fr4KCwu1bt06paenKycnR3v27Gm2/IoVK3Tbbbfpww8/VElJiVJTU3XjjTdq586dp135ttD0bJqgkYL0jgAA0O48xpiw/gJnZmbqyiuv1IsvvihJCgaDSk1N1QMPPKApU6Z85/kNDQ3q3r27XnzxRY0dO7ZV9/T7/YqPj1d1dbXi4uLCqe53qj50WOnTlkqSNv/yJvmio9r0+gAAuFVr/36H1TNSX1+v0tJSZWdnH7uA16vs7GyVlJS06hoHDx7U4cOH1aNHjxbLBAIB+f3+kFekNE1glVjeCwCADWGFkX379qmhoUGJiYkhxxMTE1VZWdmqazz66KNKTk4OCTTfVlRUpPj4eOeVmpoaTjXD0rS0V+L5NAAA2NCuq2mefvppzZs3TwsWLFBsbGyL5QoKClRdXe28KioqIlan0J4RJrECANDeosMpnJCQoKioKFVVVYUcr6qqUlJS0knPfe655/T000/r//7v/zRo0KCTlvX5fPL5fOFU7ZR5PB5FeT1qCBp6RgAAsCCsnpGYmBgNGTJExcXFzrFgMKji4mJlZWW1eN6zzz6rX/ziF1qyZImGDh166rWNkGMPy6NnBACA9hZWz4gk5efna9y4cRo6dKiGDRumGTNmqK6uTuPHj5ckjR07VikpKSoqKpIkPfPMM5o6dareeOMNpaWlOXNLunbtqq5du7bhVzl1naK8ChwJMoEVAAALwg4jo0eP1t69ezV16lRVVlZq8ODBWrJkiTOptby8XN7jJoW+/PLLqq+v109+8pOQ6xQWFurnP//56dW+jTTtNcKTewEAaH9hhxFJmjRpkiZNmtTsZytWrAj5eceOHadyi3bVNIn1MD0jAAC0O9c/m0Y6tryXh+UBAND+CCOSoqOYwAoAgC2EER2dwCqx6RkAADYQRiT5oo82Q+AwPSMAALQ3wogkX6ejD8f75nCD5ZoAAOA+hBFJsY09I98cIYwAANDeCCOSYp2eEYZpAABob4QRSbGdGntGGKYBAKDdEUZ0fM8IYQQAgPZGGJEUG300jASOMEwDAEB7I4yIYRoAAGwijIilvQAA2EQY0XFLe1lNAwBAuyOMiJ4RAABsIozouNU0TGAFAKDdEUYkdW4MI4fq6RkBAKC9EUYkxXWOliT5vzlsuSYAALgPYURSXGwnSZL/EGEEAID2RhiRFN/5aBipJowAANDuCCM6FkboGQEAoP0RRnQsjNTVN+hwAytqAABoT4QRSXGNYUSidwQAgPZGGJEU5fXo3NijK2q+PlhvuTYAALgLYaRRcnxnSdLOA99YrgkAAO5CGGmU0r0xjHx9yHJNAABwF8JIo5RuTT0jBy3XBAAAdyGMNOqTcI4kaXNljeWaAADgLoSRRhm9u0mSSr/6WsGgsVsZAABchDDS6HvJ8TrXF62vDx7Wii17bFcHAADXiLZdgTNFTLRXt16ZqtdWb9cDb6zXdQN6qvs5nRTtPZrXPB7JI0/je/M8LXzgaeaDlq7R0geeZj5o8X4tXbqZD5q7btjXbqFwe9cjnPZvScv3OzPavzldYqLVvUsndT8nRolxserdo4uivGFcAAAsI4wcJ/+G/iqrOKDSr77Wog27bVcHOCWxnby6JClO1/Y/X7mXJ+mSpDjbVQKAk/IYY874CRJ+v1/x8fGqrq5WXFxk/8faEDT667Z/aVNljaoPHVYwaGRkZIxkJLXUWkYtftCaQ0ePt3Dx5g63fI2WqnHiB+H+l2+ufpGsR3OHW65zGG0X5n/DSLV/GL8yjdduvv3rAg06cLBe+w/Wa9eBQ/rmcOgjDa7q20P3fr+vRg7oGVYvEQCcrtb+/T6lMDJz5kz9+te/VmVlpdLT0/XCCy9o2LBhLZZ/++239eSTT2rHjh3q16+fnnnmGd18882tvl97hhGgI2sIGpXvP6i/7divZZ9VafmmPWponJD9veQ4PfCDfrrxskR5GcYB0A5a+/c77Ams8+fPV35+vgoLC7Vu3Tqlp6crJydHe/Y0P+nz448/1m233aa77rpL69ev1y233KJbbrlFGzduDPfWAL5DlNejPgnn6NahqXp17FD95Wcjde/3+6pLTJT+scuvCX8o1c2//Yve/3SXE1IAwLawe0YyMzN15ZVX6sUXX5QkBYNBpaam6oEHHtCUKVNOKD969GjV1dXp/fffd45dddVVGjx4sGbNmtWqe9IzApye/XX1em31Nv3+469UGzgiSbro/HP0/4b30dUXnae+CecwhAOgzbX273dYE1jr6+tVWlqqgoIC55jX61V2drZKSkqaPaekpET5+fkhx3JycrRw4cIW7xMIBBQIBJyf/X5/ONUE8C09zonRIzmX6N4RF2n2R9s156Pt+nJvnZ5ceLSH0hft1QXdO+v8c306JyZa5/ii5Yv2yuvxyOs9uiLJ27iizOsJb4USzkz8J8S3/fiKCzQwJd7KvcMKI/v27VNDQ4MSExNDjicmJmrTpk3NnlNZWdls+crKyhbvU1RUpGnTpoVTNQCtEN+lkx66ob/uGtFHb3xSrg837dH6igMKHAnqy711+nJvne0qArAko3f3jhFG2ktBQUFIb4rf71dqaqrFGgFnl7jYTppw7UWacO1FOtwQ1K4Dh1Sx/5C+Plivg/VHVBtoUOBIw9FVZMYoaKRg47ua3tFhtbj6D67Wr2dXa/cOK4wkJCQoKipKVVVVIcerqqqUlJTU7DlJSUlhlZckn88nn88XTtUAnKJOUV5deN45uvC8c2xXBYBLhbWaJiYmRkOGDFFxcbFzLBgMqri4WFlZWc2ek5WVFVJekpYtW9ZieQAA4C5hD9Pk5+dr3LhxGjp0qIYNG6YZM2aorq5O48ePlySNHTtWKSkpKioqkiRNnjxZ1157rZ5//nmNGjVK8+bN09q1a/XKK6+07TcBAAAdUthhZPTo0dq7d6+mTp2qyspKDR48WEuWLHEmqZaXl8vrPdbhcvXVV+uNN97QE088occee0z9+vXTwoULNXDgwLb7FgAAoMNiO3gAABAREduBFQAAoC0RRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWhb0dvA1Nm8T6/X7LNQEAAK3V9Hf7uzZ77xBhpKamRpKUmppquSYAACBcNTU1io+Pb/HzDvFsmmAwqF27duncc8+Vx+Nps+v6/X6lpqaqoqKCZ95EGG3dPmjn9kE7tw/auX1Esp2NMaqpqVFycnLIQ3S/rUP0jHi9Xl1wwQURu35cXBy/6O2Etm4ftHP7oJ3bB+3cPiLVzifrEWnCBFYAAGAVYQQAAFjl6jDi8/lUWFgon89nuypnPdq6fdDO7YN2bh+0c/s4E9q5Q0xgBQAAZy9X94wAAAD7CCMAAMAqwggAALCKMAIAAKxydRiZOXOm0tLSFBsbq8zMTK1Zs8Z2lc5oq1atUl5enpKTk+XxeLRw4cKQz40xmjp1qnr16qXOnTsrOztbX3zxRUiZ/fv3a8yYMYqLi1O3bt101113qba2NqTMp59+qhEjRig2Nlapqal69tlnI/3VzhhFRUW68sorde6556pnz5665ZZbtHnz5pAy33zzjSZOnKjzzjtPXbt21Y9//GNVVVWFlCkvL9eoUaPUpUsX9ezZU4888oiOHDkSUmbFihW64oor5PP5dPHFF2vu3LmR/npnjJdfflmDBg1yNnnKysrSn//8Z+dz2jgynn76aXk8Hj344IPOMdq6bfz85z+Xx+MJeV1yySXO52d8OxuXmjdvnomJiTGzZ882//jHP8w999xjunXrZqqqqmxX7Yy1ePFi8/jjj5t3333XSDILFiwI+fzpp5828fHxZuHChebvf/+7+eEPf2j69OljDh065JS56aabTHp6uvnrX/9q/vKXv5iLL77Y3Hbbbc7n1dXVJjEx0YwZM8Zs3LjRvPnmm6Zz587mf/7nf9rra1qVk5Nj5syZYzZu3GjKysrMzTffbHr37m1qa2udMhMmTDCpqammuLjYrF271lx11VXm6quvdj4/cuSIGThwoMnOzjbr1683ixcvNgkJCaagoMAps23bNtOlSxeTn59vPvvsM/PCCy+YqKgos2TJknb9vra89957ZtGiRWbLli1m8+bN5rHHHjOdOnUyGzduNMbQxpGwZs0ak5aWZgYNGmQmT57sHKet20ZhYaH53ve+Z3bv3u289u7d63x+preza8PIsGHDzMSJE52fGxoaTHJysikqKrJYq47j22EkGAyapKQk8+tf/9o5duDAAePz+cybb75pjDHms88+M5LM3/72N6fMn//8Z+PxeMzOnTuNMca89NJLpnv37iYQCDhlHn30UTNgwIAIf6Mz0549e4wks3LlSmPM0Tbt1KmTefvtt50yn3/+uZFkSkpKjDFHQ6PX6zWVlZVOmZdfftnExcU57fqzn/3MfO973wu51+jRo01OTk6kv9IZq3v37uZ3v/sdbRwBNTU1pl+/fmbZsmXm2muvdcIIbd12CgsLTXp6erOfdYR2duUwTX19vUpLS5Wdne0c83q9ys7OVklJicWadVzbt29XZWVlSJvGx8crMzPTadOSkhJ169ZNQ4cOdcpkZ2fL6/Xqk08+ccp8//vfV0xMjFMmJydHmzdv1tdff91O3+bMUV1dLUnq0aOHJKm0tFSHDx8OaedLLrlEvXv3Dmnnyy+/XImJiU6ZnJwc+f1+/eMf/3DKHH+NpjJu/P1vaGjQvHnzVFdXp6ysLNo4AiZOnKhRo0ad0B60ddv64osvlJycrL59+2rMmDEqLy+X1DHa2ZVhZN++fWpoaAhpdElKTExUZWWlpVp1bE3tdrI2raysVM+ePUM+j46OVo8ePULKNHeN4+/hFsFgUA8++KCGDx+ugQMHSjraBjExMerWrVtI2W+383e1YUtl/H6/Dh06FImvc8bZsGGDunbtKp/PpwkTJmjBggW67LLLaOM2Nm/ePK1bt05FRUUnfEZbt53MzEzNnTtXS5Ys0csvv6zt27drxIgRqqmp6RDt3CGe2gu40cSJE7Vx40atXr3adlXOSgMGDFBZWZmqq6v1zjvvaNy4cVq5cqXtap1VKioqNHnyZC1btkyxsbG2q3NWy83Ndf49aNAgZWZm6sILL9Rbb72lzp07W6xZ67iyZyQhIUFRUVEnzCSuqqpSUlKSpVp1bE3tdrI2TUpK0p49e0I+P3LkiPbv3x9SprlrHH8PN5g0aZLef/99ffjhh7rggguc40lJSaqvr9eBAwdCyn+7nb+rDVsqExcX1yH+x9UWYmJidPHFF2vIkCEqKipSenq6fvOb39DGbai0tFR79uzRFVdcoejoaEVHR2vlypX67W9/q+joaCUmJtLWEdKtWzf1799fW7du7RC/064MIzExMRoyZIiKi4udY8FgUMXFxcrKyrJYs46rT58+SkpKCmlTv9+vTz75xGnTrKwsHThwQKWlpU6Z5cuXKxgMKjMz0ymzatUqHT582CmzbNkyDRgwQN27d2+nb2OPMUaTJk3SggULtHz5cvXp0yfk8yFDhqhTp04h7bx582aVl5eHtPOGDRtCgt+yZcsUFxenyy67zClz/DWayrj59z8YDCoQCNDGbej666/Xhg0bVFZW5ryGDh2qMWPGOP+mrSOjtrZWX375pXr16tUxfqdPewpsBzVv3jzj8/nM3LlzzWeffWbuvfde061bt5CZxAhVU1Nj1q9fb9avX28kmenTp5v169ebr776yhhzdGlvt27dzB//+Efz6aefmn/7t39rdmlvRkaG+eSTT8zq1atNv379Qpb2HjhwwCQmJpo77rjDbNy40cybN8906dLFNUt777vvPhMfH29WrFgRskTv4MGDTpkJEyaY3r17m+XLl5u1a9earKwsk5WV5XzetETvxhtvNGVlZWbJkiXm/PPPb3aJ3iOPPGI+//xzM3PmTFcthZwyZYpZuXKl2b59u/n000/NlClTjMfjMUuXLjXG0MaRdPxqGmNo67by05/+1KxYscJs377dfPTRRyY7O9skJCSYPXv2GGPO/HZ2bRgxxpgXXnjB9O7d28TExJhhw4aZv/71r7ardEb78MMPjaQTXuPGjTPGHF3e++STT5rExETj8/nM9ddfbzZv3hxyjX/961/mtttuM127djVxcXFm/PjxpqamJqTM3//+d3PNNdcYn89nUlJSzNNPP91eX9G65tpXkpkzZ45T5tChQ+b+++833bt3N126dDH//u//bnbv3h1ynR07dpjc3FzTuXNnk5CQYH7605+aw4cPh5T58MMPzeDBg01MTIzp27dvyD3Odnfeeae58MILTUxMjDn//PPN9ddf7wQRY2jjSPp2GKGt28bo0aNNr169TExMjElJSTGjR482W7dudT4/09vZY4wxp9+/AgAAcGpcOWcEAACcOQgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArPr/CJWBjrliJw4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from geopy.distance import geodesic\n",
    "from sympy import symbols, diff, sqrt, lambdify\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define the predicted distance function using SymPy\n",
    "def predicted_distance_sym(lat1, lon1, lat2, lon2):\n",
    "    delta_x = lon2 - lon1\n",
    "    delta_y = lat2 - lat1\n",
    "    return sqrt(delta_x**2 + delta_y**2)\n",
    "\n",
    "# Normalization functions for latitude and longitude\n",
    "def normalize_latitude(lat):\n",
    "    # Bring latitude within -360 to +360\n",
    "    lat = lat % 360\n",
    "    # Adjust to fit within -90 to +90\n",
    "    if lat > 90 and lat <= 270:\n",
    "        lat = 180 - lat\n",
    "    elif lat > 270:\n",
    "        lat = lat - 360\n",
    "    return lat\n",
    "\n",
    "def normalize_longitude(lon):\n",
    "    # Bring longitude within -180 to +180\n",
    "    lon = ((lon + 180) % 360) - 180\n",
    "    return lon\n",
    "\n",
    "# Function to clip gradients to prevent oscillations\n",
    "def clip_gradient(grad, max_value=1.0):\n",
    "    return np.clip(grad, -max_value, max_value)\n",
    "\n",
    "# Actual geographical coordinates of three places\n",
    "actual_lat1, actual_lon1 = 52.5200, 13.4050   # Berlin\n",
    "actual_lat2, actual_lon2 = 48.8566, 2.3522    # Paris\n",
    "actual_lat3, actual_lon3 = 51.5074, -0.1278   # London\n",
    "\n",
    "# Compute the actual distances between each pair of points\n",
    "actual_distance_12 = geodesic((actual_lat1, actual_lon1), (actual_lat2, actual_lon2)).kilometers\n",
    "actual_distance_13 = geodesic((actual_lat1, actual_lon1), (actual_lat3, actual_lon3)).kilometers\n",
    "actual_distance_23 = geodesic((actual_lat2, actual_lon2), (actual_lat3, actual_lon3)).kilometers\n",
    "print(actual_distance_12 )\n",
    "\n",
    "# Initialize the latitudes and longitudes for three locations\n",
    "lat1, lon1 = 0.01, 0.02\n",
    "lat2, lon2 = 0.03, 0.04\n",
    "lat3, lon3 = 0.05, 0.06\n",
    "\n",
    "# Define symbols for sympy\n",
    "lat1_sym, lon1_sym, lat2_sym, lon2_sym, lat3_sym, lon3_sym = symbols(\n",
    "    'lat1_sym lon1_sym lat2_sym lon2_sym lat3_sym lon3_sym')\n",
    "\n",
    "# Predicted distances with sympy symbols for each pair\n",
    "predicted_dist_12 = predicted_distance_sym(lat1_sym, lon1_sym, lat2_sym, lon2_sym)\n",
    "predicted_dist_13 = predicted_distance_sym(lat1_sym, lon1_sym, lat3_sym, lon3_sym)\n",
    "predicted_dist_23 = predicted_distance_sym(lat2_sym, lon2_sym, lat3_sym, lon3_sym)\n",
    "\n",
    "# Cost function as the sum of squared errors for all pairs\n",
    "cost_function = ((predicted_dist_12 - actual_distance_12) ** 2 +\n",
    "                 (predicted_dist_13 - actual_distance_13) ** 2 +\n",
    "                 (predicted_dist_23 - actual_distance_23) ** 2)\n",
    "\n",
    "# Compute the partial derivatives of the cost function\n",
    "grad_lat1 = diff(cost_function, lat1_sym)\n",
    "grad_lon1 = diff(cost_function, lon1_sym)\n",
    "grad_lat2 = diff(cost_function, lat2_sym)\n",
    "grad_lon2 = diff(cost_function, lon2_sym)\n",
    "grad_lat3 = diff(cost_function, lat3_sym)\n",
    "grad_lon3 = diff(cost_function, lon3_sym)\n",
    "\n",
    "# Convert the gradients to numerical functions using lambdify\n",
    "grad_lat1_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym, lat3_sym, lon3_sym), grad_lat1, 'numpy')\n",
    "grad_lon1_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym, lat3_sym, lon3_sym), grad_lon1, 'numpy')\n",
    "grad_lat2_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym, lat3_sym, lon3_sym), grad_lat2, 'numpy')\n",
    "grad_lon2_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym, lat3_sym, lon3_sym), grad_lon2, 'numpy')\n",
    "grad_lat3_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym, lat3_sym, lon3_sym), grad_lat3, 'numpy')\n",
    "grad_lon3_func = lambdify((lat1_sym, lon1_sym, lat2_sym, lon2_sym, lat3_sym, lon3_sym), grad_lon3, 'numpy')\n",
    "\n",
    "# Gradient Descent parameters\n",
    "learning_rate = 0.01\n",
    "num_iterations = 5000\n",
    "\n",
    "cost=[]\n",
    "x=[]\n",
    "# Gradient Descent loop\n",
    "for i in range(num_iterations):\n",
    "    # Calculate gradients\n",
    "    x.append(i+1)\n",
    "    #grad_lat1_val = clip_gradient(grad_lat1_func(lat1, lon1, lat2, lon2, lat3, lon3))\n",
    "    #grad_lon1_val = clip_gradient(grad_lon1_func(lat1, lon1, lat2, lon2, lat3, lon3))\n",
    "    #grad_lat2_val = clip_gradient(grad_lat2_func(lat1, lon1, lat2, lon2, lat3, lon3))\n",
    "    #grad_lon2_val = clip_gradient(grad_lon2_func(lat1, lon1, lat2, lon2, lat3, lon3))\n",
    "    #grad_lat3_val = clip_gradient(grad_lat3_func(lat1, lon1, lat2, lon2, lat3, lon3))\n",
    "    #grad_lon3_val = clip_gradient(grad_lon3_func(lat1, lon1, lat2, lon2, lat3, lon3))\n",
    "    grad_lat1_val = grad_lat1_func(lat1, lon1, lat2, lon2, lat3, lon3)\n",
    "    grad_lon1_val = grad_lon1_func(lat1, lon1, lat2, lon2, lat3, lon3)\n",
    "    grad_lat2_val = grad_lat2_func(lat1, lon1, lat2, lon2, lat3, lon3)\n",
    "    grad_lon2_val = grad_lon2_func(lat1, lon1, lat2, lon2, lat3, lon3)\n",
    "    grad_lat3_val = grad_lat3_func(lat1, lon1, lat2, lon2, lat3, lon3)\n",
    "    grad_lon3_val = grad_lon3_func(lat1, lon1, lat2, lon2, lat3, lon3)\n",
    "\n",
    "    # Update locations\n",
    "    lat1 -= learning_rate * grad_lat1_val\n",
    "    lon1 -= learning_rate * grad_lon1_val\n",
    "    lat2 -= learning_rate * grad_lat2_val\n",
    "    lon2 -= learning_rate * grad_lon2_val\n",
    "    lat3 -= learning_rate * grad_lat3_val\n",
    "    lon3 -= learning_rate * grad_lon3_val\n",
    "\n",
    "    # Normalize the latitude and longitude values\n",
    "    #lat1 = normalize_latitude(lat1)\n",
    "    #lon1 = normalize_longitude(lon1)\n",
    "    #lat2 = normalize_latitude(lat2)\n",
    "    #lon2 = normalize_longitude(lon2)\n",
    "    #lat3 = normalize_latitude(lat3)\n",
    "    #lon3 = normalize_longitude(lon3)\n",
    "\n",
    "    # Calculate updated cost\n",
    "    updated_cost = ((predicted_distance_sym(lat1, lon1, lat2, lon2) - actual_distance_12) ** 2 +\n",
    "                    (predicted_distance_sym(lat1, lon1, lat3, lon3) - actual_distance_13) ** 2 +\n",
    "                    (predicted_distance_sym(lat2, lon2, lat3, lon3) - actual_distance_23) ** 2)\n",
    "    cost.append(updated_cost)\n",
    "    # Print updated locations and cost\n",
    "    #print(f\"Iteration {i+1}:\")\n",
    "    #print(f\"Location 1: ({lat1:.6f}, {lon1:.6f})\")\n",
    "    #print(f\"Location 2: ({lat2:.6f}, {lon2:.6f})\")\n",
    "    #print(f\"Location 3: ({lat3:.6f}, {lon3:.6f})\")\n",
    "    #print(f\"Cost: {updated_cost}\\n\")\n",
    "print(f\"Location 1: ({lat1:.6f}, {lon1:.6f})\")\n",
    "print(f\"Location 2: ({lat2:.6f}, {lon2:.6f})\")\n",
    "print(f\"Location 3: ({lat3:.6f}, {lon3:.6f})\")\n",
    "#dist=geodesic((lat1,lon1),(lat2,lon2)).kilometers\n",
    "dist=predicted_distance_sym(lat1,lon1,lat2,lon2)\n",
    "print(dist)\n",
    "plt.plot(x, cost)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "68d4c1bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{(0, 1): 179.33440487789025, (0, 2): 245.2131093544739, (0, 3): 281.40151393040526, (0, 4): 184.86820889185262, (0, 5): 278.4203752471717, (0, 6): 201.67608828560105, (0, 7): 287.18577273764487, (0, 8): 129.21272236838507, (0, 9): 136.21271091363224, (0, 10): 263.929700863042, (0, 11): 102.56681350250993, (0, 12): 292.38694884909773, (1, 2): 74.32729111111344, (1, 3): 107.75756069774037, (1, 4): 49.64663220290803, (1, 5): 149.7484816832519, (1, 6): 184.18856216989573, (1, 7): 351.9364370615073, (1, 8): 61.73356848014629, (1, 9): 93.09378160749522, (1, 10): 229.4222591407563, (1, 11): 79.24018510077295, (1, 12): 231.6647822379754, (2, 3): 36.2899060618998, (2, 4): 64.53266318105935, (2, 5): 95.29764800574429, (2, 6): 186.19474203690518, (2, 7): 363.40319469219304, (2, 8): 116.72474749440087, (2, 9): 165.908622742257, (2, 10): 285.50827715806537, (2, 11): 150.7797332027635, (2, 12): 201.77374941509623, (3, 4): 100.19719061069813, (3, 5): 98.57067740461687, (3, 6): 211.74666526784688, (3, 7): 388.85918601393314, (3, 8): 153.01464750658798, (3, 9): 196.279030167369, (3, 10): 303.33898718129, (3, 11): 185.85440655334594, (3, 12): 213.17613884646408, (4, 5): 107.7459106989666, (4, 6): 139.05013925146707, (4, 7): 312.1015128721118, (4, 8): 55.91058230220265, (4, 9): 133.64535905861703, (4, 10): 278.1729860457778, (4, 11): 100.09925138355548, (4, 12): 182.03275112322288, (5, 6): 139.1204745622709, (5, 7): 308.169510992298, (5, 8): 157.74422598498097, (5, 9): 240.45333784687833, (5, 10): 376.18516502262236, (5, 11): 206.57549175316058, (5, 12): 116.0336453633846, (6, 7): 177.361773764673, (6, 8): 140.69278491503488, (6, 9): 239.54838451138178, (6, 10): 396.58857274689046, (6, 11): 184.97702648652432, (6, 12): 93.19730234418839, (7, 8): 297.0864210932367, (7, 9): 383.0866149188769, (7, 10): 536.8734001178462, (7, 11): 326.65755696732236, (7, 12): 211.08483688228927, (8, 9): 100.08756645218632, (8, 10): 256.1618921868582, (8, 11): 51.21853470855889, (8, 12): 208.11056790649857, (9, 10): 157.42831995183784, (9, 11): 56.77118052993751, (9, 12): 307.18610475814734, (10, 11): 213.70391916892197, (10, 12): 459.04453603663666, (11, 12): 258.2133598908325}\n",
      "0\n",
      "100\n",
      "200\n",
      "300\n",
      "400\n",
      "500\n",
      "600\n",
      "700\n",
      "800\n",
      "900\n",
      "1000\n",
      "1100\n",
      "1200\n",
      "1300\n",
      "1400\n",
      "1500\n",
      "1600\n",
      "1700\n",
      "1800\n",
      "1900\n",
      "2000\n",
      "2100\n",
      "2200\n",
      "2300\n",
      "2400\n",
      "2500\n",
      "2600\n",
      "2700\n",
      "2800\n",
      "2900\n",
      "3000\n",
      "3100\n",
      "3200\n",
      "3300\n",
      "3400\n",
      "3500\n",
      "3600\n",
      "3700\n",
      "3800\n",
      "3900\n",
      "4000\n",
      "4100\n",
      "4200\n",
      "4300\n",
      "4400\n",
      "4500\n",
      "4600\n",
      "4700\n",
      "4800\n",
      "4900\n",
      "5000\n",
      "5100\n",
      "5200\n",
      "5300\n",
      "5400\n",
      "5500\n",
      "5600\n",
      "5700\n",
      "5800\n",
      "5900\n",
      "6000\n",
      "6100\n",
      "6200\n",
      "6300\n",
      "6400\n",
      "6500\n",
      "6600\n",
      "6700\n",
      "6800\n",
      "6900\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJS0lEQVR4nO3dd3RUdf7/8dekTRLSqAklEJoU6aCQ4ApKNCL6A3dF5LBSFGyIsFjW6Cqi68ay+IVVRNEVUBcjoICyCCIIrhKVLkV6CSJJQEgFkpD5/P6IGRkTMIGZuSnPxzlz1rnzuXfe94aTvPZT7rUZY4wAAACqCR+rCwAAAHAnwg0AAKhWCDcAAKBaIdwAAIBqhXADAACqFcINAACoVgg3AACgWiHcAACAaoVwAwAAqhXCDQBUMn379lXfvn2tLgOosgg3QCW0b98+3XPPPWrRooUCAwMVFham3r17a9q0aTp9+rTbv+/UqVN6+umntXr1arcf22oxMTG66aabnO8ry7nu2LFDTz/9tA4ePGhpHUB15Gd1AQBc/fe//9XgwYNlt9s1fPhwdejQQQUFBfrqq6/0yCOPaPv27Zo5c6Zbv/PUqVOaPHmyJFX7HoPKcq47duzQ5MmT1bdvX8XExLh89tlnn1lTFFBNEG6ASuTAgQO6/fbb1axZM61atUoNGzZ0fjZ27Fjt3btX//3vfy2ssPI5e/asHA6HAgICLK0jLy9PtWrVcsuxrD4XoKpjWAqoRF588UXl5ubq3//+t0uwKdGqVSuNHz/e+f7s2bN69tln1bJlS9ntdsXExOjxxx9Xfn6+y37r169XQkKC6tWrp6CgIDVv3lx33nmnJOngwYOqX7++JGny5Mmy2Wyy2Wx6+umnL1jr/v37NXjwYNWpU0fBwcHq1auXS/BKT0+Xn5+fs5fkXLt27ZLNZtOrr77q3JaZmakJEyYoOjpadrtdrVq10gsvvCCHw+Fsc/DgQdlsNv3zn//U1KlTnee9Y8eOC9Z67v6/d647d+7Urbfeqjp16igwMFA9evTQxx9/7HKc2bNny2azac2aNbr//vvVoEEDNWnSRJJ06NAh3X///WrTpo2CgoJUt25dDR482GX4afbs2Ro8eLAk6ZprrnHWUTJUVtacm4yMDN11112KjIxUYGCgOnfurDlz5pQ6v5LrM3PmTOf1ueKKK7Ru3TqXtmlpaRo1apSaNGkiu92uhg0bauDAgQyToVqg5waoRD755BO1aNFCcXFx5Wo/evRozZkzR7feeqseeughffvtt0pKStIPP/yghQsXSir+o3j99derfv36euyxxxQREaGDBw/qo48+kiTVr19fM2bM0H333adbbrlFf/zjHyVJnTp1Ou/3pqenKy4uTqdOndKDDz6ounXras6cOfp//+//acGCBbrlllsUGRmpPn36aN68eZo0aZLL/h988IF8fX2df+BPnTqlPn366MiRI7rnnnvUtGlTrV27VomJiTp69KimTp3qsv+sWbN05swZ3X333bLb7apTp065rtfvnev27dvVu3dvNW7cWI899phq1aqlefPmadCgQfrwww91yy23uBzv/vvvV/369fXUU08pLy9PkrRu3TqtXbtWt99+u5o0aaKDBw9qxowZ6tu3r3bs2KHg4GBdffXVevDBB/Wvf/1Ljz/+uNq1aydJzv/9rdOnT6tv377au3evHnjgATVv3lzz58/XyJEjlZmZ6RJ4JWnu3LnKycnRPffcI5vNphdffFF//OMftX//fvn7+0uS/vSnP2n79u0aN26cYmJilJGRoRUrVig1NbXUMBlQ5RgAlUJWVpaRZAYOHFiu9ps3bzaSzOjRo122P/zww0aSWbVqlTHGmIULFxpJZt26dec91rFjx4wkM2nSpHJ994QJE4wk87///c+5LScnxzRv3tzExMSYoqIiY4wxb7zxhpFktm7d6rJ/+/btzbXXXut8/+yzz5patWqZ3bt3u7R77LHHjK+vr0lNTTXGGHPgwAEjyYSFhZmMjIxy1dqsWTMzYMCAcp1rv379TMeOHc2ZM2ec2xwOh4mLizOtW7d2bps1a5aRZK666ipz9uxZl2OcOnWq1HFTUlKMJPPOO+84t82fP99IMl988UWp9n369DF9+vRxvp86daqRZN577z3ntoKCAhMbG2tCQkJMdna2MebX61O3bl1z4sQJZ9vFixcbSeaTTz4xxhhz8uRJI8m89NJLpb4bqA4YlgIqiezsbElSaGhoudovXbpUkjRx4kSX7Q899JAkOYeIIiIiJElLlixRYWGhO0rV0qVLdeWVV+qqq65ybgsJCdHdd9+tgwcPOoeJ/vjHP8rPz08ffPCBs922bdu0Y8cODRkyxLlt/vz5+sMf/qDatWvr+PHjzld8fLyKior05Zdfunz/n/70J+fwkrucOHFCq1at0m233aacnBxnDT///LMSEhK0Z88eHTlyxGWfMWPGyNfX12VbUFCQ878LCwv1888/q1WrVoqIiNDGjRsvqralS5cqKipKQ4cOdW7z9/fXgw8+qNzcXK1Zs8al/ZAhQ1S7dm3n+z/84Q+SiocSS2oMCAjQ6tWrdfLkyYuqCajManS4+fLLL3XzzTerUaNGstlsWrRoUYWPYYzRP//5T1122WWy2+1q3LixnnvuOfcXi2ovLCxMkpSTk1Ou9ocOHZKPj49atWrlsj0qKkoRERE6dOiQJKlPnz7605/+pMmTJ6tevXoaOHCgZs2aVWpeTkUcOnRIbdq0KbW9ZFil5Lvr1aunfv36ad68ec42H3zwgfz8/JxDQpK0Z88eLVu2TPXr13d5xcfHSyoeWjtX8+bNL7r289m7d6+MMXryySdL1VEyrFaeOk6fPq2nnnrKOXeoXr16ql+/vjIzM5WVlXVRtR06dEitW7eWj4/rr+zfXu8STZs2dXlfEnRKgozdbtcLL7ygTz/9VJGRkbr66qv14osvKi0t7aLqAyqbGj3nJi8vT507d9add97p8ou2IsaPH6/PPvtM//znP9WxY0edOHFCJ06ccHOlqAnCwsLUqFEjbdu2rUL72Wy23/18wYIF+uabb/TJJ59o+fLluvPOOzVlyhR98803CgkJuZSyf9ftt9+uUaNGafPmzerSpYvmzZunfv36qV69es42DodD1113nR599NEyj3HZZZe5vD+3d8RdSiYuP/zww0pISCizzW+DZFl1jBs3TrNmzdKECRMUGxur8PBw2Ww23X777S6Toz3pt71JJYwxzv+eMGGCbr75Zi1atEjLly/Xk08+qaSkJK1atUpdu3b1Sp2Ap9TocNO/f3/179//vJ/n5+friSee0Pvvv6/MzEx16NBBL7zwgnMVww8//KAZM2Zo27Ztzv8X64n/R4ma46abbtLMmTOVkpKi2NjYC7Zt1qyZHA6H9uzZ4zIRNT09XZmZmWrWrJlL+169eqlXr1567rnnNHfuXA0bNkzJyckaPXr07waksr57165dpbbv3LnT+XmJQYMG6Z577nEOTe3evVuJiYku+7Vs2VK5ubnOnhpPOt+5tmjRQlLxcM+l1LFgwQKNGDFCU6ZMcW47c+aMMjMzy1VHWZo1a6bvv/9eDofDpfemrOtdES1bttRDDz2khx56SHv27FGXLl00ZcoUvffeexd1PKCyqNHDUr/ngQceUEpKipKTk/X9999r8ODBuuGGG7Rnzx5Jv65sWbJkiZo3b66YmBiNHj2anhtctEcffVS1atXS6NGjlZ6eXurzffv2adq0aZKkG2+8UZJKrSR6+eWXJUkDBgyQVDwUce7/Y5ekLl26SJJzaCo4OFiSSv0BPp8bb7xR3333nVJSUpzb8vLyNHPmTMXExKh9+/bO7REREUpISNC8efOUnJysgIAADRo0yOV4t912m1JSUrR8+fJS35WZmamzZ8+Wq67yON+5NmjQQH379tUbb7yho0ePltrv2LFj5Tq+r69vqev9yiuvqKioyGVbyT1xynPNb7zxRqWlpbnMXTp79qxeeeUVhYSEqE+fPuWqrcSpU6d05swZl20tW7ZUaGjoJQ1XApVFje65uZDU1FTNmjVLqampatSokaTi7uply5Zp1qxZ+sc//qH9+/fr0KFDmj9/vt555x0VFRXpL3/5i2699VatWrXK4jNAVdSyZUvNnTtXQ4YMUbt27VzuULx27Vrn8l9J6ty5s0aMGKGZM2cqMzNTffr00Xfffac5c+Zo0KBBuuaaayRJc+bM0WuvvaZbbrlFLVu2VE5Ojt58802FhYU5A1JQUJDat2+vDz74QJdddpnq1KmjDh06qEOHDmXW+dhjj+n9999X//799eCDD6pOnTqaM2eODhw4oA8//LDU3JAhQ4boz3/+s1577TUlJCQ4JzmXeOSRR/Txxx/rpptu0siRI9W9e3fl5eVp69atWrBggQ4ePOgyjHUpLnSu06dP11VXXaWOHTtqzJgxatGihdLT05WSkqIff/xRW7Zs+d3j33TTTXr33XcVHh6u9u3bKyUlRZ9//rnq1q3r0q5Lly7y9fXVCy+8oKysLNntdl177bVq0KBBqWPefffdeuONNzRy5Eht2LBBMTExWrBggb7++mtNnTq13JPQS+zevVv9+vXTbbfdpvbt28vPz08LFy5Uenq6br/99godC6iULF2rVYlIMgsXLnS+X7JkiZFkatWq5fLy8/Mzt912mzHGmDFjxhhJZteuXc79NmzYYCSZnTt3evsUUI3s3r3bjBkzxsTExJiAgAATGhpqevfubV555RWXZcqFhYVm8uTJpnnz5sbf399ER0ebxMRElzYbN240Q4cONU2bNjV2u900aNDA3HTTTWb9+vUu37l27VrTvXt3ExAQUK5l4fv27TO33nqriYiIMIGBgebKK680S5YsKbNtdna2CQoKKrWc+Vw5OTkmMTHRtGrVygQEBJh69eqZuLg4889//tMUFBQYY35d6lyRJcy/XQr+e+e6b98+M3z4cBMVFWX8/f1N48aNzU033WQWLFjgbFOyFLys5fUnT540o0aNMvXq1TMhISEmISHB7Ny50zRr1syMGDHCpe2bb75pWrRoYXx9fV2Whf92KbgxxqSnpzuPGxAQYDp27GhmzZrl0uZC1+fc8zx+/LgZO3asadu2ralVq5YJDw83PXv2NPPmzbvwxQSqCJsxv+k/raFsNpsWLlzo7C7/4IMPNGzYMG3fvr3U5LyQkBBFRUVp0qRJ+sc//uGyvPb06dMKDg7WZ599puuuu86bpwAAAMSw1Hl17dpVRUVFysjIcN4j4rd69+6ts2fPat++fWrZsqWk4u5e6eIn+AEAgEtTo3tucnNztXfvXknFYebll1/WNddcozp16qhp06b685//rK+//lpTpkxR165ddezYMa1cuVKdOnXSgAED5HA4dMUVVygkJERTp06Vw+HQ2LFjFRYWxlN9AQCwSI0ON6tXr3ZOujzXiBEjNHv2bBUWFurvf/+73nnnHR05ckT16tVTr169NHnyZHXs2FGS9NNPP2ncuHH67LPPVKtWLfXv319Tpkwp97NuAACAe9XocAMAAKof7nMDAACqFcINAACoVmrcaimHw6GffvpJoaGhFb7lPAAAsIYxRjk5OWrUqFGpG4X+Vo0LNz/99JOio6OtLgMAAFyEw4cPq0mTJhdsU+PCTcltyg8fPqywsDCLqwEAAOWRnZ2t6Ojocj1upMaFm5KhqLCwMMINAABVTHmmlDChGAAAVCuEGwAAUK0QbgAAQLVSacLN888/L5vNpgkTJlyw3fz589W2bVsFBgaqY8eOWrp0qXcKBAAAVUKlCDfr1q3TG2+8oU6dOl2w3dq1azV06FDddddd2rRpkwYNGqRBgwZp27ZtXqoUAABUdpaHm9zcXA0bNkxvvvmmateufcG206ZN0w033KBHHnlE7dq107PPPqtu3brp1Vdf9VK1AACgsrM83IwdO1YDBgxQfHz877ZNSUkp1S4hIUEpKSnn3Sc/P1/Z2dkuLwAAUH1Zep+b5ORkbdy4UevWrStX+7S0NEVGRrpsi4yMVFpa2nn3SUpK0uTJky+pTgAAUHVY1nNz+PBhjR8/Xv/5z38UGBjose9JTExUVlaW83X48GGPfRcAALCeZT03GzZsUEZGhrp16+bcVlRUpC+//FKvvvqq8vPz5evr67JPVFSU0tPTXbalp6crKirqvN9jt9tlt9vdWzwAAKi0LOu56devn7Zu3arNmzc7Xz169NCwYcO0efPmUsFGkmJjY7Vy5UqXbStWrFBsbKy3ygYAAJWcZT03oaGh6tChg8u2WrVqqW7dus7tw4cPV+PGjZWUlCRJGj9+vPr06aMpU6ZowIABSk5O1vr16zVz5kyv1w8AACony1dLXUhqaqqOHj3qfB8XF6e5c+dq5syZ6ty5sxYsWKBFixaVCklWKDjr0E+Zp3Uk87TVpQAAUKPZjDHG6iK8KTs7W+Hh4crKynLrU8HXHTyhwa+nKKZusFY/co3bjgsAACr297tS99xUJXa/4ktZcNZhcSUAANRshBs3CSgJN0WEGwAArES4cZMA3+JLmU/PDQAAliLcuIndv3jpOuEGAABrEW7cpKTnpuCsQzVsjjYAAJUK4cZNSubcSFJhEeEGAACrEG7cxH5OuMk/W2RhJQAA1GyEGzcpGZaSWA4OAICVCDdu4uNjk7+vTRLLwQEAsBLhxo3OnVQMAACsQbhxI5aDAwBgPcKNG9FzAwCA9Qg3blSyHJyeGwAArEO4caNfww1LwQEAsArhxo14MjgAANYj3LhRAOEGAADLEW7ciCeDAwBgPcKNG5UsBafnBgAA6xBu3Mi5FJw7FAMAYBnCjRsxoRgAAOsRbtzIzlJwAAAsR7hxI1ZLAQBgPcKNGxFuAACwHuHGjVgKDgCA9Qg3bmT3J9wAAGA1wo0bBfj+cp8bloIDAGAZwo0bMecGAADrEW7c6Nel4IQbAACsQrhxo197brjPDQAAViHcuBHDUgAAWI9w40YMSwEAYD3CjRvxbCkAAKxnabiZMWOGOnXqpLCwMIWFhSk2NlaffvrpedvPnj1bNpvN5RUYGOjFii/MOSzFUnAAACzjZ+WXN2nSRM8//7xat24tY4zmzJmjgQMHatOmTbr88svL3CcsLEy7du1yvrfZbN4q93c573NDzw0AAJaxNNzcfPPNLu+fe+45zZgxQ9988815w43NZlNUVJQ3yqsw7lAMAID1Ks2cm6KiIiUnJysvL0+xsbHnbZebm6tmzZopOjpaAwcO1Pbt2y943Pz8fGVnZ7u8PKXk2VL03AAAYB3Lw83WrVsVEhIiu92ue++9VwsXLlT79u3LbNumTRu9/fbbWrx4sd577z05HA7FxcXpxx9/PO/xk5KSFB4e7nxFR0d76lScc27ouQEAwDo2Y4yxsoCCggKlpqYqKytLCxYs0FtvvaU1a9acN+Ccq7CwUO3atdPQoUP17LPPltkmPz9f+fn5zvfZ2dmKjo5WVlaWwsLC3HYekrTvWK76TVmj0EA/bX06wa3HBgCgJsvOzlZ4eHi5/n5bOudGkgICAtSqVStJUvfu3bVu3TpNmzZNb7zxxu/u6+/vr65du2rv3r3nbWO322W3291W74WwFBwAAOtZPiz1Ww6Hw6Wn5UKKioq0detWNWzY0MNVlc+5S8Et7hADAKDGsrTnJjExUf3791fTpk2Vk5OjuXPnavXq1Vq+fLkkafjw4WrcuLGSkpIkSc8884x69eqlVq1aKTMzUy+99JIOHTqk0aNHW3kaTvZfloIbI511GPn7Vp5l6gAA1BSWhpuMjAwNHz5cR48eVXh4uDp16qTly5fruuuukySlpqbKx+fXzqWTJ09qzJgxSktLU+3atdW9e3etXbu2XPNzvKFkKbhUPKnY37fSdYwBAFDtWT6h2NsqMiGpohwOoxaPL5UkbXzyOtWpFeDW4wMAUFNV5O83XQtu5ONjk59P8VAUk4oBALAG4cbNfr3XTZHFlQAAUDMRbtyM5eAAAFiLcONm3KUYAABrEW7c7Nx73QAAAO8j3LiZ3a/4Xjf5hYQbAACsQLhxM+eTwem5AQDAEoQbNwtgQjEAAJYi3LgZS8EBALAW4cbNWAoOAIC1CDduRrgBAMBahBs34z43AABYi3DjZiVLwem5AQDAGoQbNytZCs6EYgAArEG4cTO7P8NSAABYiXDjZoH+v9yhmHADAIAlCDduFvjLhOIzhQxLAQBgBcKNm9n9ebYUAABWIty4Wcl9bs4woRgAAEsQbtyMnhsAAKxFuHGzQHpuAACwFOHGzei5AQDAWoQbN6PnBgAAaxFu3Kyk5+YMPTcAAFiCcONmgX48fgEAACsRbtyMOTcAAFiLcONmgf703AAAYCXCjZsF+jHnBgAAKxFu3MxOzw0AAJYi3LhZSc9NYZFRkcNYXA0AADUP4cbNSnpuJHpvAACwgqXhZsaMGerUqZPCwsIUFham2NhYffrppxfcZ/78+Wrbtq0CAwPVsWNHLV261EvVlo/9l54biXk3AABYwdJw06RJEz3//PPasGGD1q9fr2uvvVYDBw7U9u3by2y/du1aDR06VHfddZc2bdqkQYMGadCgQdq2bZuXKz8/Xx+b/H1tkui5AQDACjZjTKWaGFKnTh299NJLuuuuu0p9NmTIEOXl5WnJkiXObb169VKXLl30+uuvl+v42dnZCg8PV1ZWlsLCwtxW97k6TlqunPyz+uLhvmper5ZHvgMAgJqkIn+/K82cm6KiIiUnJysvL0+xsbFltklJSVF8fLzLtoSEBKWkpHijxHIrmXdzppCeGwAAvM3P6gK2bt2q2NhYnTlzRiEhIVq4cKHat29fZtu0tDRFRka6bIuMjFRaWtp5j5+fn6/8/Hzn++zsbPcUfgEl827yzzLnBgAAb7O856ZNmzbavHmzvv32W913330aMWKEduzY4bbjJyUlKTw83PmKjo5227HPh54bAACsY3m4CQgIUKtWrdS9e3clJSWpc+fOmjZtWplto6KilJ6e7rItPT1dUVFR5z1+YmKisrKynK/Dhw+7tf6yBNJzAwCAZSwPN7/lcDhchpHOFRsbq5UrV7psW7FixXnn6EiS3W53LjUveXlaID03AABYxtI5N4mJierfv7+aNm2qnJwczZ07V6tXr9by5cslScOHD1fjxo2VlJQkSRo/frz69OmjKVOmaMCAAUpOTtb69es1c+ZMK0+jFObcAABgHUvDTUZGhoYPH66jR48qPDxcnTp10vLly3XddddJklJTU+Xj82vnUlxcnObOnau//e1vevzxx9W6dWstWrRIHTp0sOoUykTPDQAA1rE03Pz73/++4OerV68utW3w4MEaPHiwhypyD3puAACwTqWbc1MdlPTc5NNzAwCA1xFuPICeGwAArEO48QDm3AAAYB3CjQfY/Yt7bgg3AAB4H+HGAwL9fplzw7AUAABeR7jxAHpuAACwDuHGA+z03AAAYBnCjQcE0nMDAIBlCDceQM8NAADWIdx4AD03AABYh3DjAfTcAABgHcKNB/zac0O4AQDA2wg3HvBrzw3DUgAAeBvhxgNKem7y6bkBAMDrCDcewIRiAACsQ7jxgOCA4nBzmnADAIDXEW48oKTn5nRhkYwxFlcDAEDNQrjxgKBfem6MYTk4AADeRrjxgJKngkvS6QKGpgAA8CbCjQf4+foowLf40jLvBgAA7yLceEigP+EGAAArEG48pGTeDcNSAAB4F+HGQ4K41w0AAJYg3HjIucvBAQCA9xBuPIRhKQAArEG48ZAgem4AALAE4cZDmHMDAIA1CDceEsiwFAAAliDceMivw1I8fgEAAG8i3HgIc24AALAG4cZDSlZLMecGAADvItx4iPM+N8y5AQDAqywNN0lJSbriiisUGhqqBg0aaNCgQdq1a9cF95k9e7ZsNpvLKzAw0EsVlx/DUgAAWMPScLNmzRqNHTtW33zzjVasWKHCwkJdf/31ysvLu+B+YWFhOnr0qPN16NAhL1VcfkE8OBMAAEv4Wfnly5Ytc3k/e/ZsNWjQQBs2bNDVV1993v1sNpuioqI8Xd4lcc65YVgKAACvqlRzbrKysiRJderUuWC73NxcNWvWTNHR0Ro4cKC2b9/ujfIqJCigODfScwMAgHdVmnDjcDg0YcIE9e7dWx06dDhvuzZt2ujtt9/W4sWL9d5778nhcCguLk4//vhjme3z8/OVnZ3t8vIG5twAAGANS4elzjV27Fht27ZNX3311QXbxcbGKjY21vk+Li5O7dq10xtvvKFnn322VPukpCRNnjzZ7fX+niBWSwEAYIlK0XPzwAMPaMmSJfriiy/UpEmTCu3r7++vrl27au/evWV+npiYqKysLOfr8OHD7ij5dwUFMKEYAAArWNpzY4zRuHHjtHDhQq1evVrNmzev8DGKioq0detW3XjjjWV+brfbZbfbL7XUCuM+NwAAWMPScDN27FjNnTtXixcvVmhoqNLS0iRJ4eHhCgoKkiQNHz5cjRs3VlJSkiTpmWeeUa9evdSqVStlZmbqpZde0qFDhzR69GjLzqMszLkBAMAaloabGTNmSJL69u3rsn3WrFkaOXKkJCk1NVU+Pr+Onp08eVJjxoxRWlqaateure7du2vt2rVq3769t8ouFx6/AACANWzGGGN1Ed6UnZ2t8PBwZWVlKSwszGPfk3mqQF2eWSFJ2vNcf/n7VorpTQAAVEkV+fvNX1wPKZlzI9F7AwCANxFuPMTu5yObrfi/mXcDAID3EG48xGazOScVnylwWFwNAAA1B+HGg1gxBQCA9xFuPCiQcAMAgNcRbjyoZDk4N/IDAMB7CDce5JxzQ88NAABeQ7jxIObcAADgfYQbDwpkWAoAAK8j3HhQkD9PBgcAwNsINx7EnBsAALyPcONBJaulTjEsBQCA1xBuPCjIv/ih6wxLAQDgPYQbDwou6bnJP2txJQAA1ByEGw8KtjMsBQCAtxFuPCjYn3ADAIC3EW48KDigeM7NqQKGpQAA8BbCjQcxLAUAgPcRbjwomKXgAAB4HeHGgxiWAgDA+wg3HhTMs6UAAPA6wo0HlYSbPMINAABeQ7jxoJJhKXpuAADwHsKNB5X03BQUOVRY5LC4GgAAagbCjQeVPDhTYsUUAADeQrjxoABfH/n52CQxNAUAgLcQbjzIZrM5e2/yWA4OAIBXXFS4eeaZZ3Tq1KlS20+fPq1nnnnmkouqTmoxqRgAAK+6qHAzefJk5ebmltp+6tQpTZ48+ZKLqk64SzEAAN51UeHGGCObzVZq+5YtW1SnTp1LLqo6YVgKAADv8qtI49q1a8tms8lms+myyy5zCThFRUXKzc3Vvffe6/YiqzKGpQAA8K4KhZupU6fKGKM777xTkydPVnh4uPOzgIAAxcTEKDY21u1FVmVBDEsBAOBVFQo3I0aMkCQ1b95cvXv3lp9fhXYvJSkpSR999JF27typoKAgxcXF6YUXXlCbNm0uuN/8+fP15JNP6uDBg2rdurVeeOEF3XjjjZdUi6fUspeEG4alAADwhouacxMaGqoffvjB+X7x4sUaNGiQHn/8cRUUFJT7OGvWrNHYsWP1zTffaMWKFSosLNT111+vvLy88+6zdu1aDR06VHfddZc2bdqkQYMGadCgQdq2bdvFnIrHBfmXPBmcnhsAALzhosLNPffco927d0uS9u/fryFDhig4OFjz58/Xo48+Wu7jLFu2TCNHjtTll1+uzp07a/bs2UpNTdWGDRvOu8+0adN0ww036JFHHlG7du307LPPqlu3bnr11Vcv5lQ8jtVSAAB410WFm927d6tLly6SioeI+vTpo7lz52r27Nn68MMPL7qYrKwsSbrgiquUlBTFx8e7bEtISFBKSkqZ7fPz85Wdne3y8qbgkmGpfIalAADwhoteCu5wFD8I8vPPP3fOd4mOjtbx48cvqhCHw6EJEyaod+/e6tChw3nbpaWlKTIy0mVbZGSk0tLSymyflJSk8PBw5ys6Ovqi6rtYwSXDUoX03AAA4A0XFW569Oihv//973r33Xe1Zs0aDRgwQJJ04MCBUsGjvMaOHatt27YpOTn5ovY/n8TERGVlZTlfhw8fduvxf49zWIqeGwAAvOKiljtNnTpVw4YN06JFi/TEE0+oVatWkqQFCxYoLi6uwsd74IEHtGTJEn355Zdq0qTJBdtGRUUpPT3dZVt6erqioqLKbG+322W32ytck7s4h6WYcwMAgFdcVLjp1KmTtm7dWmr7Sy+9JF9f33IfxxijcePGaeHChVq9erWaN2/+u/vExsZq5cqVmjBhgnPbihUrKu39dUp6bk4zLAUAgFdc0o1qNmzY4FwS3r59e3Xr1q1C+48dO1Zz587V4sWLFRoa6pw3Ex4erqCgIEnS8OHD1bhxYyUlJUmSxo8frz59+mjKlCkaMGCAkpOTtX79es2cOfNSTsVjgn+5Q3Eew1IAAHjFRYWbjIwMDRkyRGvWrFFERIQkKTMzU9dcc42Sk5NVv379ch1nxowZkqS+ffu6bJ81a5ZGjhwpSUpNTZWPz69Tg+Li4jR37lz97W9/0+OPP67WrVtr0aJFF5yEbCWWggMA4F0XFW7GjRun3Nxcbd++Xe3atZMk7dixQyNGjNCDDz6o999/v1zHMcb8bpvVq1eX2jZ48GANHjy4QjVbhWEpAAC866LCzbJly/T55587g41UPCw1ffp0XX/99W4rrjr4dViKcAMAgDdc1FJwh8Mhf3//Utv9/f2d979BMWfPDc+WAgDAKy4q3Fx77bUaP368fvrpJ+e2I0eO6C9/+Yv69evntuKqA+dTwQuLyjUMBwAALs1FhZtXX31V2dnZiomJUcuWLdWyZUs1b95c2dnZeuWVV9xdY5VW65dhKWOkM4X0agEA4GkXNecmOjpaGzdu1Oeff66dO3dKktq1a1fqmU+Qgvx/ve/PqYKzzp4cAADgGRXquVm1apXat2+v7Oxs2Ww2XXfddRo3bpzGjRunK664Qpdffrn+97//earWKsnHx+YMOCwHBwDA8yoUbqZOnaoxY8YoLCys1Gfh4eG655579PLLL7utuOqCe90AAOA9FQo3W7Zs0Q033HDez6+//npt2LDhkouqbpyTilkxBQCAx1Uo3KSnp5e5BLyEn5+fjh07dslFVTclk4rpuQEAwPMqFG4aN26sbdu2nffz77//Xg0bNrzkoqqbIIalAADwmgqFmxtvvFFPPvmkzpw5U+qz06dPa9KkSbrpppvcVlx1UcvOsBQAAN5SoaXgf/vb3/TRRx/psssu0wMPPKA2bdpIknbu3Knp06erqKhITzzxhEcKrcqC/BmWAgDAWyoUbiIjI7V27Vrdd999SkxMdN5x12azKSEhQdOnT1dkZKRHCq3KSnpu8vLpuQEAwNMqfBO/Zs2aaenSpTp58qT27t0rY4xat26t2rVre6K+aqGWnYdnAgDgLRd1h2JJql27tq644gp31lJthZSEG+bcAADgcRf1bClUTEm4yTlDuAEAwNMIN17w67AU4QYAAE8j3HhB6C/hJpdwAwCAxxFuvKCk5yaXYSkAADyOcOMFIYH03AAA4C2EGy8I+eU+N4QbAAA8j3DjBSH24oeNMqEYAADPI9x4QcmwVA7hBgAAjyPceEFIQHG4KTjrUMFZh8XVAABQvRFuvKDk2VISQ1MAAHga4cYL/Hx9FOhffKmZVAwAgGcRbrykZFIx4QYAAM8i3HgJy8EBAPAOwo2XcCM/AAC8g3DjJbUCeAQDAADeQLjxktBAngwOAIA3EG68pBZPBgcAwCssDTdffvmlbr75ZjVq1Eg2m02LFi26YPvVq1fLZrOVeqWlpXmn4EsQ8ku4yWFYCgAAj7I03OTl5alz586aPn16hfbbtWuXjh496nw1aNDAQxW6T0m4YVgKAADP8rPyy/v376/+/ftXeL8GDRooIiLC/QV5UAjDUgAAeEWVnHPTpUsXNWzYUNddd52+/vrrC7bNz89Xdna2y8sKzLkBAMA7qlS4adiwoV5//XV9+OGH+vDDDxUdHa2+fftq48aN590nKSlJ4eHhzld0dLQXK/4V97kBAMA7LB2Wqqg2bdqoTZs2zvdxcXHat2+f/u///k/vvvtumfskJiZq4sSJzvfZ2dmWBBznsBQTigEA8KgqFW7KcuWVV+qrr7467+d2u112u92LFZUtlJ4bAAC8okoNS5Vl8+bNatiwodVl/K6wwOIHZ2adLrS4EgAAqjdLe25yc3O1d+9e5/sDBw5o8+bNqlOnjpo2barExEQdOXJE77zzjiRp6tSpat68uS6//HKdOXNGb731llatWqXPPvvMqlMot/Cg4nCTTbgBAMCjLA0369ev1zXXXON8XzI3ZsSIEZo9e7aOHj2q1NRU5+cFBQV66KGHdOTIEQUHB6tTp076/PPPXY5RWYX9Em7yCop0tsghP98q32kGAEClZDPGGKuL8Kbs7GyFh4crKytLYWFhXvvewiKHWj/xqSRp45PXqU6tAK99NwAAVV1F/n7TfeAl/r4+qhXgK4mhKQAAPIlw40UlQ1PZZwg3AAB4CuHGi0pWTGWfZjk4AACeQrjxopIVUywHBwDAcwg3XhQWVLw4jWEpAAA8h3DjRb8OSxFuAADwFMKNF4UxLAUAgMcRbryI1VIAAHge4caLwn55eCarpQAA8BzCjRexWgoAAM8j3HgRw1IAAHge4caLWC0FAIDnEW686NdhKebcAADgKYQbL+ImfgAAeB7hxotK5twUnHXoTGGRxdUAAFA9EW68KCTAT74+NklS5il6bwAA8ATCjRf5+NhUO7i49+bkqQKLqwEAoHoi3HhZ7eAASdLJPMINAACeQLjxstq1isPNCXpuAADwCMKNl9Wh5wYAAI8i3HiZs+cmjwnFAAB4AuHGy+rUYkIxAACeRLjxspIJxScYlgIAwCMIN15W55dhKXpuAADwDMKNl/0654ZwAwCAJxBuvIzVUgAAeBbhxsvqcJ8bAAA8inDjZSXDUmcKHTpdwMMzAQBwN8KNl9UK8FWAb/Flp/cGAAD3I9x4mc1mcw5NHc/Jt7gaAACqH8KNBRqE2SVJxwg3AAC4HeHGAg1Ci8NNBuEGAAC3szTcfPnll7r55pvVqFEj2Ww2LVq06Hf3Wb16tbp16ya73a5WrVpp9uzZHq/T3eqHBkqSMnLOWFwJAADVj6XhJi8vT507d9b06dPL1f7AgQMaMGCArrnmGm3evFkTJkzQ6NGjtXz5cg9X6l703AAA4Dl+Vn55//791b9//3K3f/3119W8eXNNmTJFktSuXTt99dVX+r//+z8lJCR4qky3K5lzk5FNuAEAwN2q1JyblJQUxcfHu2xLSEhQSkrKeffJz89Xdna2y8tq9UNKJhQzLAUAgLtVqXCTlpamyMhIl22RkZHKzs7W6dOny9wnKSlJ4eHhzld0dLQ3Sr2gBmElc27ouQEAwN2qVLi5GImJicrKynK+Dh8+bHVJzjk3x3Pz5XAYi6sBAKB6sXTOTUVFRUUpPT3dZVt6errCwsIUFBRU5j52u112u90b5ZVbvV+GpQqLjDJPFzpv6gcAAC5dleq5iY2N1cqVK122rVixQrGxsRZVdHEC/HycgYbl4AAAuJel4SY3N1ebN2/W5s2bJRUv9d68ebNSU1MlFQ8pDR8+3Nn+3nvv1f79+/Xoo49q586deu211zRv3jz95S9/saL8S+JcDs6KKQAA3MrScLN+/Xp17dpVXbt2lSRNnDhRXbt21VNPPSVJOnr0qDPoSFLz5s313//+VytWrFDnzp01ZcoUvfXWW1VqGXiJ+tzrBgAAj7B0zk3fvn1lzPkn1JZ19+G+fftq06ZNHqzKOyJ/WTGVns2wFAAA7lSl5txUJ1GEGwAAPIJwY5HI8OJwk5ZFuAEAwJ0INxaJ/GXODT03AAC4F+HGIlElPTeEGwAA3IpwY5GSOTfHcvJ1tshhcTUAAFQfhBuL1A2xy9fHJoeRfs4rsLocAACqDcKNRXx9bM4b+TGpGAAA9yHcWKjk6eDMuwEAwH0INxaKCmPFFAAA7ka4sVDJpGKGpQAAcB/CjYUiWQ4OAIDbEW4sVNJzw5PBAQBwH8KNhSKZUAwAgNsRbizkfDI4c24AAHAbwo2FSh7BkJN/Vnn5Zy2uBgCA6oFwY6EQu59C7H6SGJoCAMBdCDcWiyy51w1DUwAAuAXhxmI8HRwAAPci3FgsKixIknSUnhsAANyCcGOxRhHFPTdHs05bXAkAANUD4cZijSKKe25+yqTnBgAAdyDcWKzhL3Nufsqk5wYAAHcg3Fjs154bwg0AAO5AuLFYSc9N9hlu5AcAgDsQbiwWGuiv0MDiG/kxqRgAgEtHuKkEGoUXD00dYVIxAACXjHBTCTiXgzPvBgCAS0a4qQQalkwq5kZ+AABcMsJNJdCYFVMAALgN4aYS4F43AAC4D+GmEoiuEyxJSj1xyuJKAACo+gg3lUCzX8LNT5mnVXDWYXE1AABUbZUi3EyfPl0xMTEKDAxUz5499d1335237ezZs2Wz2VxegYGBXqzW/eqH2hXk7yuHkY4wNAUAwCWxPNx88MEHmjhxoiZNmqSNGzeqc+fOSkhIUEZGxnn3CQsL09GjR52vQ4cOebFi97PZbGr6S+/NoZ/zLK4GAICqzfJw8/LLL2vMmDEaNWqU2rdvr9dff13BwcF6++23z7uPzWZTVFSU8xUZGenFij2jaV3m3QAA4A6WhpuCggJt2LBB8fHxzm0+Pj6Kj49XSkrKeffLzc1Vs2bNFB0drYEDB2r79u3nbZufn6/s7GyXV2UUU7ek54ZwAwDApbA03Bw/flxFRUWlel4iIyOVlpZW5j5t2rTR22+/rcWLF+u9996Tw+FQXFycfvzxxzLbJyUlKTw83PmKjo52+3m4Q9O6tSQRbgAAuFSWD0tVVGxsrIYPH64uXbqoT58++uijj1S/fn298cYbZbZPTExUVlaW83X48GEvV1w+zZzLwZlzAwDApfCz8svr1asnX19fpaenu2xPT09XVFRUuY7h7++vrl27au/evWV+brfbZbfbL7lWT2t2zpwbh8PIx8dmcUUAAFRNlvbcBAQEqHv37lq5cqVzm8Ph0MqVKxUbG1uuYxQVFWnr1q1q2LChp8r0ikYRQfL1selMoUMZOflWlwMAQJVl+bDUxIkT9eabb2rOnDn64YcfdN999ykvL0+jRo2SJA0fPlyJiYnO9s8884w+++wz7d+/Xxs3btSf//xnHTp0SKNHj7bqFNzC39fHuRx8/7Fci6sBAKDqsnRYSpKGDBmiY8eO6amnnlJaWpq6dOmiZcuWOScZp6amysfn1wx28uRJjRkzRmlpaapdu7a6d++utWvXqn379ladgtu0ahCiA8fztCcjV3Gt6lldDgAAVZLNGGOsLsKbsrOzFR4erqysLIWFhVldjosXl+3Ua6v36c+9murvgzpaXQ4AAJVGRf5+Wz4shV+1jgyRJO1JZ1gKAICLRbipRFo3CJUk7c0g3AAAcLEIN5VIy/ohstmkn/MK9HMuK6YAALgYhJtKJCjAV01qB0mi9wYAgItFuKlkSoamdqblWFwJAABVE+GmkunQOFyS9P2PWRZXAgBA1US4qWQ6NykON1t+zLS2EAAAqijCTSXTqUmEJGnfsVzl5p+1thgAAKogwk0lUz/UrsYRQTJG2srQFAAAFUa4qYQ6RxcPTa0/eMLiSgAAqHoIN5VQbMvi50r9b+9xiysBAKDqIdxUQle3Lg43Gw+dZN4NAAAVRLiphJrVraWmdYJ11mGUsu9nq8sBAKBKIdxUUte2bSBJWrr1qMWVAABQtRBuKqmbOzeSJH22PU2nC4osrgYAgKqDcFNJdWsaoSa1g5RXUKTPdqRZXQ4AAFUG4aaSstls+lO3JpKkOWsPWlsMAABVCOGmEhvWq6kCfH20MTVTm1JPWl0OAABVAuGmEmsQGuice/Pqqr0WVwMAQNVAuKnk7uvbUr4+Nq3cmaGvuakfAAC/i3BTybVqEKI7ejWTJD25eJvOFLJyCgCACyHcVAF/ib9MDULt2n8sT1M+22V1OQAAVGqEmyogPNhfSX/sKEl666sDPFATAIALINxUEf3aRerW7k1kjDQ+ebNO5hVYXRIAAJUS4aYKefKm9mpaJ1hHMk/rweRNKnIYq0sCAKDSIdxUIeFB/nrjju4K9PfR//Yc19Mfb5cxBBwAAM5FuKli2jUM00u3dpbNJr37zSFN/mQHPTgAAJyDcFMF3dy5kZJuKZ5gPHvtQY2c9Z1+PHnK4qoAAKgcbKaGjWtkZ2crPDxcWVlZCgsLs7qcS/Lxlp/06IItOlPokN3PRwO7NNIfWtdXu4Zhal6vlnx9bFaXCACAW1Tk7zfhporbnZ6jJxdt07cHXJeHB/r7qE1kqNpGhaldw1B1aVpbnZuEy2Yj8AAAqh7CzQVUt3AjScYYrT90Up9s+Unf/5ilnWnZOlPoKNWuRb1aejihjW7s2NCCKgEAuHiEmwuojuHmt4ocRod+ztMPR3O0My1bPxzN1tp9P+tUQfGjG+6+uoUS+7elFwcAUGVU5O93pZhQPH36dMXExCgwMFA9e/bUd999d8H28+fPV9u2bRUYGKiOHTtq6dKlXqq0avD1salF/RAN6NRQD13fRm+NuELrnojX2GtaSpJmfrlfD83booKzpXt3AACo6iwPNx988IEmTpyoSZMmaePGjercubMSEhKUkZFRZvu1a9dq6NChuuuuu7Rp0yYNGjRIgwYN0rZt27xcedVSy+6nRxLa6qVbO8nXx6aPNh3RiLe/U3r2GatLAwDArSwflurZs6euuOIKvfrqq5Ikh8Oh6OhojRs3To899lip9kOGDFFeXp6WLFni3NarVy916dJFr7/++u9+X00Ylvo9q3dlaOx/NiqvoEihgX4a1rOZrr6sniLDAhXg6yNGqwAAlyLAz0cNQgPdesyK/P32c+s3V1BBQYE2bNigxMRE5zYfHx/Fx8crJSWlzH1SUlI0ceJEl20JCQlatGhRme3z8/OVn5/vfJ+dnX3phVdxfds00MKxvfXI/C3a8mOWXl+zT6+v2Wd1WQCAaqJb0wh9dH9vy77f0nBz/PhxFRUVKTIy0mV7ZGSkdu7cWeY+aWlpZbZPS0srs31SUpImT57snoKrkcsiQ/XR/b21Yke6Ptnyk7b8mKms04XMwwEAXDJ/X2tnvVgabrwhMTHRpacnOztb0dHRFlZUefj62HRDhyjd0CHK6lIAAHAbS8NNvXr15Ovrq/T0dJft6enpiooq+w9uVFRUhdrb7XbZ7Xb3FAwAACo9S/uNAgIC1L17d61cudK5zeFwaOXKlYqNjS1zn9jYWJf2krRixYrztgcAADWL5cNSEydO1IgRI9SjRw9deeWVmjp1qvLy8jRq1ChJ0vDhw9W4cWMlJSVJksaPH68+ffpoypQpGjBggJKTk7V+/XrNnDnTytMAAACVhOXhZsiQITp27JieeuoppaWlqUuXLlq2bJlz0nBqaqp8fH7tYIqLi9PcuXP1t7/9TY8//rhat26tRYsWqUOHDladAgAAqEQsv8+Nt3GfGwAAqp4q9/gFAAAAdyHcAACAaoVwAwAAqhXCDQAAqFYINwAAoFoh3AAAgGqFcAMAAKoVwg0AAKhWCDcAAKBasfzxC95WckPm7OxsiysBAADlVfJ3uzwPVqhx4SYnJ0eSFB0dbXElAACgonJychQeHn7BNjXu2VIOh0M//fSTQkNDZbPZ3Hrs7OxsRUdH6/DhwzXyuVU1/fwlrkFNP3+Ja8D51+zzlzx3DYwxysnJUaNGjVweqF2WGtdz4+PjoyZNmnj0O8LCwmrsP2qJ85e4BjX9/CWuAedfs89f8sw1+L0emxJMKAYAANUK4QYAAFQrhBs3stvtmjRpkux2u9WlWKKmn7/ENajp5y9xDTj/mn3+UuW4BjVuQjEAAKje6LkBAADVCuEGAABUK4QbAABQrRBuAABAtUK4cZPp06crJiZGgYGB6tmzp7777jurS7ooX375pW6++WY1atRINptNixYtcvncGKOnnnpKDRs2VFBQkOLj47Vnzx6XNidOnNCwYcMUFhamiIgI3XXXXcrNzXVp8/333+sPf/iDAgMDFR0drRdffNHTp1YuSUlJuuKKKxQaGqoGDRpo0KBB2rVrl0ubM2fOaOzYsapbt65CQkL0pz/9Senp6S5tUlNTNWDAAAUHB6tBgwZ65JFHdPbsWZc2q1evVrdu3WS329WqVSvNnj3b06dXLjNmzFCnTp2cN+CKjY3Vp59+6vy8up//bz3//POy2WyaMGGCc1t1vwZPP/20bDaby6tt27bOz6v7+UvSkSNH9Oc//1l169ZVUFCQOnbsqPXr1zs/r86/C2NiYkr9/G02m8aOHSupivz8DS5ZcnKyCQgIMG+//bbZvn27GTNmjImIiDDp6elWl1ZhS5cuNU888YT56KOPjCSzcOFCl8+ff/55Ex4ebhYtWmS2bNli/t//+3+mefPm5vTp0842N9xwg+ncubP55ptvzP/+9z/TqlUrM3ToUOfnWVlZJjIy0gwbNsxs27bNvP/++yYoKMi88cYb3jrN80pISDCzZs0y27ZtM5s3bzY33nijadq0qcnNzXW2uffee010dLRZuXKlWb9+venVq5eJi4tzfn727FnToUMHEx8fbzZt2mSWLl1q6tWrZxITE51t9u/fb4KDg83EiRPNjh07zCuvvGJ8fX3NsmXLvHq+Zfn444/Nf//7X7N7926za9cu8/jjjxt/f3+zbds2Y0z1P/9zfffddyYmJsZ06tTJjB8/3rm9ul+DSZMmmcsvv9wcPXrU+Tp27Jjz8+p+/idOnDDNmjUzI0eONN9++63Zv3+/Wb58udm7d6+zTXX+XZiRkeHys1+xYoWRZL744gtjTNX4+RNu3ODKK680Y8eOdb4vKioyjRo1MklJSRZWdel+G24cDoeJiooyL730knNbZmamsdvt5v333zfGGLNjxw4jyaxbt87Z5tNPPzU2m80cOXLEGGPMa6+9ZmrXrm3y8/Odbf7617+aNm3aePiMKi4jI8NIMmvWrDHGFJ+vv7+/mT9/vrPNDz/8YCSZlJQUY0xxQPTx8TFpaWnONjNmzDBhYWHOc3700UfN5Zdf7vJdQ4YMMQkJCZ4+pYtSu3Zt89Zbb9Wo88/JyTGtW7c2K1asMH369HGGm5pwDSZNmmQ6d+5c5mc14fz/+te/mquuuuq8n9e034Xjx483LVu2NA6Ho8r8/BmWukQFBQXasGGD4uPjndt8fHwUHx+vlJQUCytzvwMHDigtLc3lXMPDw9WzZ0/nuaakpCgiIkI9evRwtomPj5ePj4++/fZbZ5urr75aAQEBzjYJCQnatWuXTp486aWzKZ+srCxJUp06dSRJGzZsUGFhocs1aNu2rZo2bepyDTp27KjIyEhnm4SEBGVnZ2v79u3ONuceo6RNZfs3U1RUpOTkZOXl5Sk2NrZGnf/YsWM1YMCAUnXWlGuwZ88eNWrUSC1atNCwYcOUmpoqqWac/8cff6wePXpo8ODBatCggbp27ao333zT+XlN+l1YUFCg9957T3feeadsNluV+fkTbi7R8ePHVVRU5PJDlKTIyEilpaVZVJVnlJzPhc41LS1NDRo0cPncz89PderUcWlT1jHO/Y7KwOFwaMKECerdu7c6dOggqbi+gIAARUREuLT97TX4vfM7X5vs7GydPn3aE6dTIVu3blVISIjsdrvuvfdeLVy4UO3bt68x55+cnKyNGzcqKSmp1Gc14Rr07NlTs2fP1rJlyzRjxgwdOHBAf/jDH5STk1Mjzn///v2aMWOGWrdureXLl+u+++7Tgw8+qDlz5kiqWb8LFy1apMzMTI0cOVJS1fn3X+OeCg6U19ixY7Vt2zZ99dVXVpfidW3atNHmzZuVlZWlBQsWaMSIEVqzZo3VZXnF4cOHNX78eK1YsUKBgYFWl2OJ/v37O/+7U6dO6tmzp5o1a6Z58+YpKCjIwsq8w+FwqEePHvrHP/4hSeratau2bdum119/XSNGjLC4Ou/697//rf79+6tRo0ZWl1Ih9Nxconr16snX17fUTPH09HRFRUVZVJVnlJzPhc41KipKGRkZLp+fPXtWJ06ccGlT1jHO/Q6rPfDAA1qyZIm++OILNWnSxLk9KipKBQUFyszMdGn/22vwe+d3vjZhYWGV4o9HQECAWrVqpe7duyspKUmdO3fWtGnTasT5b9iwQRkZGerWrZv8/Pzk5+enNWvW6F//+pf8/PwUGRlZ7a/Bb0VEROiyyy7T3r17a8S/gYYNG6p9+/Yu29q1a+ccmqspvwsPHTqkzz//XKNHj3Zuqyo/f8LNJQoICFD37t21cuVK5zaHw6GVK1cqNjbWwsrcr3nz5oqKinI51+zsbH377bfOc42NjVVmZqY2bNjgbLNq1So5HA717NnT2ebLL79UYWGhs82KFSvUpk0b1a5d20tnUzZjjB544AEtXLhQq1atUvPmzV0+7969u/z9/V2uwa5du5SamupyDbZu3eryi23FihUKCwtz/sKMjY11OUZJm8r6b8bhcCg/P79GnH+/fv20detWbd682fnq0aOHhg0b5vzv6n4Nfis3N1f79u1Tw4YNa8S/gd69e5e6BcTu3bvVrFkzSTXjd6EkzZo1Sw0aNNCAAQOc26rMz98t05JruOTkZGO3283s2bPNjh07zN13320iIiJcZopXFTk5OWbTpk1m06ZNRpJ5+eWXzaZNm8yhQ4eMMcXLHyMiIszixYvN999/bwYOHFjm8seuXbuab7/91nz11VemdevWLssfMzMzTWRkpLnjjjvMtm3bTHJysgkODrZ8+aMxxtx3330mPDzcrF692mUp5KlTp5xt7r33XtO0aVOzatUqs379ehMbG2tiY2Odn5csg7z++uvN5s2bzbJly0z9+vXLXAb5yCOPmB9++MFMnz690iyDfeyxx8yaNWvMgQMHzPfff28ee+wxY7PZzGeffWaMqf7nX5ZzV0sZU/2vwUMPPWRWr15tDhw4YL7++msTHx9v6tWrZzIyMowx1f/8v/vuO+Pn52eee+45s2fPHvOf//zHBAcHm/fee8/Zprr/LiwqKjJNmzY1f/3rX0t9VhV+/oQbN3nllVdM06ZNTUBAgLnyyivNN998Y3VJF+WLL74wkkq9RowYYYwpXgL55JNPmsjISGO3202/fv3Mrl27XI7x888/m6FDh5qQkBATFhZmRo0aZXJyclzabNmyxVx11VXGbrebxo0bm+eff95bp3hBZZ27JDNr1ixnm9OnT5v777/f1K5d2wQHB5tbbrnFHD161OU4Bw8eNP379zdBQUGmXr165qGHHjKFhYUubb744gvTpUsXExAQYFq0aOHyHVa68847TbNmzUxAQICpX7++6devnzPYGFP9z78svw031f0aDBkyxDRs2NAEBASYxo0bmyFDhrjc46W6n78xxnzyySemQ4cOxm63m7Zt25qZM2e6fF7dfxcuX77cSCp1TsZUjZ+/zRhj3NMHBAAAYD3m3AAAgGqFcAMAAKoVwg0AAKhWCDcAAKBaIdwAAIBqhXADAACqFcINAACoVgg3AGqcmJgYTZ061eoyAHgI4QaAR40cOVKDBg2SJPXt21cTJkzw2nfPnj1bERERpbavW7dOd999t9fqAOBdflYXAAAVVVBQoICAgIvev379+m6sBkBlQ88NAK8YOXKk1qxZo2nTpslms8lms+ngwYOSpG3btql///4KCQlRZGSk7rjjDh0/fty5b9++ffXAAw9owoQJqlevnhISEiRJL7/8sjp27KhatWopOjpa999/v3JzcyVJq1ev1qhRo5SVleX8vqefflpS6WGp1NRUDRw4UCEhIQoLC9Ntt92m9PR05+dPP/20unTponfffVcxMTEKDw/X7bffrpycHM9eNAAXhXADwCumTZum2NhYjRkzRkePHtXRo0cVHR2tzMxMXXvtteratavWr1+vZcuWKT09XbfddpvL/nPmzFFAQIC+/vprvf7665IkHx8f/etf/9L27ds1Z84crVq1So8++qgkKS4uTlOnTlVYWJjz+x5++OFSdTkcDg0cOFAnTpzQmjVrtGLFCu3fv19Dhgxxabdv3z4tWrRIS5Ys0ZIlS7RmzRo9//zzHrpaAC4Fw1IAvCI8PFwBAQEKDg5WVFSUc/urr76qrl276h//+Idz29tvv63o6Gjt3r1bl112mSSpdevWevHFF12Oee78nZiYGP3973/Xvffeq9dee00BAQEKDw+XzWZz+b7fWrlypbZu3aoDBw4oOjpakvTOO+/o8ssv17p163TFFVdIKg5Bs2fPVmhoqCTpjjvu0MqVK/Xcc89d2oUB4Hb03ACw1JYtW/TFF18oJCTE+Wrbtq2k4t6SEt27dy+17+eff65+/fqpcePGCg0N1R133KGff/5Zp06dKvf3//DDD4qOjnYGG0lq3769IiIi9MMPPzi3xcTEOIONJDVs2FAZGRkVOlcA3kHPDQBL5ebm6uabb9YLL7xQ6rOGDRs6/7tWrVounx08eFA33XST7rvvPj333HOqU6eOvvrqK911110qKChQcHCwW+v09/d3eW+z2eRwONz6HQDcg3ADwGsCAgJUVFTksq1bt2768MMPFRMTIz+/8v9K2rBhgxwOh6ZMmSIfn+JO6Hnz5v3u9/1Wu3btdPjwYR0+fNjZe7Njxw5lZmaqffv25a4HQOXBsBQAr4mJidG3336rgwcP6vjx43I4HBo7dqxOnDihoUOHat26ddq3b5+WL1+uUaNGXTCYtGrVSoWFhXrllVe0f/9+vfvuu86Jxud+X25urlauXKnjx4+XOVwVHx+vjh07atiwYdq4caO+++47DR8+XH369FGPHj3cfg0AeB7hBoDXPPzww/L19VX79u1Vv359paamqlGjRvr6669VVFSk66+/Xh07dtSECRMUERHh7JEpS+fOnfXyyy/rhRdeUIcOHfSf//xHSUlJLm3i4uJ07733asiQIapfv36pCclS8fDS4sWLVbt2bV199dWKj49XixYt9MEHH7j9/AF4h80YY6wuAgAAwF3ouQEAANUK4QYAAFQrhBsAAFCtEG4AAEC1QrgBAADVCuEGAABUK4QbAABQrRBuAABAtUK4AQAA1QrhBgAAVCuEGwAAUK0QbgAAQLXy/wGnbNf5jgJVBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from geopy.distance import geodesic\n",
    "from sympy import symbols, diff, sqrt, lambdify\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define the predicted distance function using SymPy\n",
    "def predicted_distance_sym(lat1, lon1, lat2, lon2):\n",
    "    delta_x = lon2 - lon1\n",
    "    delta_y = lat2 - lat1\n",
    "    return sqrt(delta_x**2 + delta_y**2)\n",
    "\n",
    "# Normalization functions for latitude and longitude\n",
    "def normalize_latitude(lat):\n",
    "    lat = lat % 360\n",
    "    if lat > 90 and lat <= 270:\n",
    "        lat = 180 - lat\n",
    "    elif lat > 270:\n",
    "        lat = lat - 360\n",
    "    return lat\n",
    "\n",
    "def normalize_longitude(lon):\n",
    "    lon = ((lon + 180) % 360) - 180\n",
    "    return lon\n",
    "\n",
    "# Function to clip gradients to prevent oscillations\n",
    "def clip_gradient(grad, max_value=1.0):\n",
    "    return np.clip(grad, -max_value, max_value)\n",
    "\n",
    "# Define actual geographical coordinates for 10 cities (randomly chosen for illustration)\n",
    "actual_coords = [\n",
    "    (52.631175943871, 1.4438059026964032),  # Norwich\n",
    "    (50.910307779600494, -1.685596568246761), #Southhampton\n",
    "    (50.70943994894834, -3.352854870350599), #Exeter\n",
    "    (50.37549715914587, -3.9886187080727304), #Plymouth\n",
    "    (51.44713035442605, -2.444620848400215), #Bristol\n",
    "    (51.67537358796746, -4.918247021684257), #Pembroke\n",
    "    (53.40128004270814, -3.221778430465314), #Liverpool\n",
    "    (55.954637922470745, -2.815774829240603), #Edinburgh\n",
    "    (51.76, -1.247), #Oxford\n",
    "    (50.76455315081721, 0.4288200714415762), #Eastbourne\n",
    "    (48.8566, 2.3522),   # Paris\n",
    "    (51.5074, -0.1278),  # London\n",
    "    (53.319328516885, -5.470531808221821), #Dublin\n",
    "    \n",
    "]\n",
    "\n",
    "# Compute actual distances between each pair of points\n",
    "num_cities = len(actual_coords)\n",
    "actual_distances = {}\n",
    "for i in range(num_cities):\n",
    "    for j in range(i + 1, num_cities):\n",
    "        actual_distances[(i, j)] = geodesic(actual_coords[i], actual_coords[j]).miles\n",
    "print(actual_distances)\n",
    "\n",
    "# Initialize the latitudes and longitudes randomly for 10 locations\n",
    "#latitudes = np.random.uniform(-90, 90, num_cities)\n",
    "#longitudes = np.random.uniform(-180, 180, num_cities)\n",
    "#latitudes = np.random.uniform(-1, 1, num_cities)\n",
    "#longitudes = np.random.uniform(-1, 1, num_cities)\n",
    "latitudes = np.random.random(num_cities)\n",
    "longitudes = np.random.random(num_cities)\n",
    "\n",
    "# Define symbols for sympy\n",
    "lat_symbols = symbols(f'lat0:{num_cities}') #####Mangesh : indexed beginnign with 0\n",
    "lon_symbols = symbols(f'lon0:{num_cities}')\n",
    "\n",
    "# Predicted distances with sympy symbols for each pair\n",
    "predicted_distances = {}\n",
    "for i in range(num_cities):\n",
    "    for j in range(i + 1, num_cities):\n",
    "        predicted_distances[(i, j)] = predicted_distance_sym(\n",
    "            lat_symbols[i], lon_symbols[i], lat_symbols[j], lon_symbols[j])\n",
    "\n",
    "# Cost function as the sum of squared errors for all pairs\n",
    "cost_function = 0\n",
    "for i in range(num_cities):\n",
    "    for j in range(i + 1, num_cities):\n",
    "        cost_function += (predicted_distances[(i, j)] - actual_distances[(i, j)]) ** 2\n",
    "\n",
    "#cost_function = sum([(predicted_distances[(i, j)] - actual_distances[(i, j)]) ** 2\n",
    "#                     for i in range(num_cities) for j in range(i + 1, num_cities)])\n",
    "\n",
    "#cost_function = sum((predicted_distances[(i, j)] - actual_distances[(i, j)]) ** 2\n",
    "#                    for i in range(num_cities) for j in range(i + 1, num_cities))\n",
    "\n",
    "# Compute the partial derivatives of the cost function\n",
    "grad_lat_funcs = []\n",
    "grad_lon_funcs = []\n",
    "for i in range(num_cities):\n",
    "    grad_lat = diff(cost_function, lat_symbols[i])\n",
    "    grad_lon = diff(cost_function, lon_symbols[i])\n",
    "    grad_lat_funcs.append(lambdify(lat_symbols + lon_symbols, grad_lat, 'numpy'))\n",
    "    grad_lon_funcs.append(lambdify(lat_symbols + lon_symbols, grad_lon, 'numpy'))\n",
    "\n",
    "# Gradient Descent parameters\n",
    "learning_rate = 0.5\n",
    "num_iterations = 7000\n",
    "cost_history = []  # To store cost values over iterations\n",
    "\n",
    "# Gradient Descent loop\n",
    "for iteration in range(num_iterations):\n",
    "    # Collect all current latitudes and longitudes\n",
    "    current_lat_lon = np.concatenate((latitudes, longitudes))\n",
    "\n",
    "    # Calculate gradients for all cities\n",
    "    grad_lats = np.array([clip_gradient(grad_lat_func(*current_lat_lon)) for grad_lat_func in grad_lat_funcs])\n",
    "    grad_lons = np.array([clip_gradient(grad_lon_func(*current_lat_lon)) for grad_lon_func in grad_lon_funcs])\n",
    "\n",
    "    # Update all locations\n",
    "    latitudes -= learning_rate * grad_lats\n",
    "    longitudes -= learning_rate * grad_lons\n",
    "\n",
    "    # Normalize the latitude and longitude values\n",
    "    #latitudes = np.array([normalize_latitude(lat) for lat in latitudes])\n",
    "    #longitudes = np.array([normalize_longitude(lon) for lon in longitudes])\n",
    "\n",
    "    # Calculate updated cost\n",
    "    updated_cost = 0\n",
    "    for i in range(num_cities):\n",
    "        for j in range(i + 1, num_cities):\n",
    "            updated_cost += (predicted_distance_sym(latitudes[i], longitudes[i], latitudes[j], longitudes[j]) - \n",
    "                        actual_distances[(i, j)]) ** 2\n",
    "    #updated_cost = sum((predicted_distance_sym(latitudes[i], longitudes[i], latitudes[j], longitudes[j]) - \n",
    "       #                 actual_distances[(i, j)]) ** 2\n",
    "        #               for i in range(num_cities) for j in range(i + 1, num_cities))\n",
    "\n",
    "    # Append cost to history for plotting\n",
    "    cost_history.append(float(updated_cost))\n",
    "    \n",
    "    if iteration %100 == 0:\n",
    "        print(iteration)\n",
    "\n",
    "    # Print updated locations and cost\n",
    "    #print(f\"Iteration {iteration + 1}:\")\n",
    "    #for i in range(num_cities):\n",
    "        #print(f\"Location {i+1}: ({latitudes[i]:.6f}, {longitudes[i]:.6f})\")\n",
    "        #for j in range(i + 1, num_cities):\n",
    "            #print(predicted_distance_sym(latitudes[i], longitudes[i], latitudes[j], longitudes[j]))\n",
    "    #print(f\"Cost: {updated_cost}\\n\")\n",
    "\n",
    "\n",
    "\n",
    "# Plotting the cost history\n",
    "plt.plot(cost_history)\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('Cost')\n",
    "plt.title('Cost over Iterations')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1744f27c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "179.587332661204\n",
      "245.192617621716\n",
      "281.830582593153\n",
      "185.494798269730\n",
      "278.946082743110\n",
      "201.566641468071\n",
      "286.445843753566\n",
      "129.343265297859\n",
      "137.104537464585\n",
      "263.783845234383\n",
      "102.741628280412\n",
      "238.805762927901\n",
      "74.4427038113642\n",
      "108.128883088339\n",
      "50.2615647321685\n",
      "150.591464244847\n",
      "184.019047895465\n",
      "351.855108127317\n",
      "62.0925045318956\n",
      "93.0137451831755\n",
      "229.386490126473\n",
      "79.2339618219690\n",
      "134.961104488118\n",
      "36.7577916446513\n",
      "63.8672797391620\n",
      "95.5028511036602\n",
      "185.245584654360\n",
      "362.754627885411\n",
      "116.558458300435\n",
      "166.006004563562\n",
      "286.212285775712\n",
      "150.650162222406\n",
      "102.615936382379\n",
      "100.023712341019\n",
      "98.9388995815437\n",
      "211.310711518658\n",
      "388.762336756360\n",
      "153.315645494595\n",
      "196.554911015743\n",
      "304.189409380568\n",
      "186.116314371938\n",
      "120.597198057987\n",
      "107.802172356417\n",
      "138.488613645151\n",
      "311.831489763861\n",
      "56.4034480085696\n",
      "134.388276210993\n",
      "278.765617074358\n",
      "100.538245024900\n",
      "85.6697306471936\n",
      "138.726501930733\n",
      "308.045886804284\n",
      "158.171480138751\n",
      "241.254001546564\n",
      "377.137623169917\n",
      "207.012813329269\n",
      "45.5181120458499\n",
      "177.682722822733\n",
      "140.085862177132\n",
      "239.714211797995\n",
      "396.210462588935\n",
      "184.485393047907\n",
      "96.2407100205532\n",
      "296.464921176130\n",
      "383.243690183447\n",
      "536.153725135654\n",
      "326.030537939756\n",
      "270.676600461942\n",
      "100.766297769081\n",
      "256.395028442466\n",
      "51.2387339383761\n",
      "125.779522185171\n",
      "156.815490094682\n",
      "57.5423758101161\n",
      "218.554634636263\n",
      "213.774552735842\n",
      "364.237755712790\n",
      "176.731282966691\n"
     ]
    }
   ],
   "source": [
    "for i in range(num_cities):\n",
    "        #print(f\"Location {i+1}: ({latitudes[i]:.6f}, {longitudes[i]:.6f})\")\n",
    "        for j in range(i + 1, num_cities):\n",
    "            print(predicted_distance_sym(latitudes[i], longitudes[i], latitudes[j], longitudes[j]))\n",
    "        plt.plot(latitudes[i],longitudes[i],'.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "fe9e8320",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Location 0: (50.341535, 89.563229)\n",
      "Location 1: (-58.959504, -52.932009)\n",
      "Location 2: (-129.293209, -77.322296)\n",
      "Location 3: (-153.958036, -104.576315)\n",
      "Location 4: (-98.205510, -21.531755)\n",
      "Location 5: (-206.002069, -20.431651)\n",
      "Location 6: (-150.494806, 106.706015)\n",
      "Location 7: (-159.874336, 284.141001)\n",
      "Location 8: (-50.514657, 8.583548)\n",
      "Location 9: (33.834623, -46.543993)\n",
      "Location 10: (145.793894, -156.344810)\n",
      "Location 11: (0.000000, 0.000000)\n",
      "Location 12: (-176.179392, 13.955940)\n",
      "Not Reflected\n"
     ]
    }
   ],
   "source": [
    "%matplotlib qt\n",
    "\n",
    "def quadrant(theta2):\n",
    "    if theta2>0 and theta2<=np.pi/2:\n",
    "        return 1\n",
    "    elif theta2>np.pi/2 and theta2<=np.pi:\n",
    "        return 2\n",
    "    elif theta2>np.pi and theta2<=np.pi*3/2:\n",
    "        return 3\n",
    "    elif theta2>np.pi*3/2 and theta2<=2*np.pi:\n",
    "        return 4\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "pp = ['Norwich','Southhampton', 'Exeter','Plymouth','Bristol','Pembroke','Liverpool','Edinbrugh','Oxford','Eastbourne','Paris', 'London','Dublin']\n",
    "for i in range(num_cities):\n",
    "    print(f\"Location {i}: ({latitudes[i]:.6f}, {longitudes[i]:.6f})\")\n",
    "    #plt.plot(latitudes[i],longitudes[i],'.')\n",
    "    #plt.text(latitudes[i],longitudes[i],pp[i]+str(i),fontsize = 10)\n",
    "    #for j in range(i + 1, num_cities):\n",
    "    #    print(predicted_distance_sym(latitudes[i], longitudes[i], latitudes[j], longitudes[j]))\n",
    "    \n",
    "#Set London as the center\n",
    "idx_London = pp.index('London')\n",
    "idx_Paris = pp.index('Paris')\n",
    "idx_Norwich = pp.index('Norwich')\n",
    "for i in range(num_cities):\n",
    "    latitudes[i]  -= latitudes[idx_London]\n",
    "    longitudes[i] -= longitudes[idx_London]\n",
    "    \n",
    "#Check \n",
    "theta_Paris = np.arctan2(longitudes[idx_Paris], latitudes[idx_Paris])\n",
    "theta_Norwich = np.arctan2(longitudes[idx_Norwich], latitudes[idx_Norwich])\n",
    "\n",
    "#If paris rotated by +90 deg takes it to the same quadrant as Norwich then no reflection necessary\n",
    "quadrant_paris = quadrant(np.arctan2(latitudes[idx_Paris],-longitudes[idx_Paris]))\n",
    "quadrant_Norwich = quadrant(theta_Norwich)\n",
    "r = 1\n",
    "if quadrant_Norwich == quadrant_paris:\n",
    "    print(\"Not Reflected\")\n",
    "else:\n",
    "    print(\"Reflected\")\n",
    "    r = -1\n",
    "for i in range(num_cities):\n",
    "    latitudes[i] *= r\n",
    "\n",
    "#Fix the angle between London and Paris\n",
    "theta_actual = 313*np.pi/180\n",
    "theta_sim = np.arctan2(longitudes[idx_Paris], latitudes[idx_Paris])\n",
    "theta_rot = theta_actual - theta_sim\n",
    "\n",
    "#rotate the coordinates of all the cities by theta_rot\n",
    "\n",
    "for i in range(num_cities):\n",
    "    latitudes[i], longitudes[i] = np.cos(theta_rot)*latitudes[i]-np.sin(theta_rot)*longitudes[i], np.sin(theta_rot)*latitudes[i]+np.cos(theta_rot)*longitudes[i]\n",
    "plt.axis('equal')\n",
    "for i in range(num_cities):\n",
    "    \n",
    "    plt.plot(latitudes[i],longitudes[i],'.')\n",
    "    plt.text(latitudes[i],longitudes[i],pp[i],fontsize = 10)   \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d06376bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb3929159d0>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.plot(latitudes[i],longitudes[i],'.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "id": "700a712e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: matplotlib in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (3.5.2)\n",
      "Requirement already satisfied: basemap in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (1.4.1)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (1.4.2)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (9.2.0)\n",
      "Requirement already satisfied: cycler>=0.10 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (0.11.0)\n",
      "Requirement already satisfied: python-dateutil>=2.7 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (2.8.2)\n",
      "Requirement already satisfied: numpy>=1.17 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (1.24.3)\n",
      "Requirement already satisfied: fonttools>=4.22.0 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (4.25.0)\n",
      "Requirement already satisfied: packaging>=20.0 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (21.3)\n",
      "Requirement already satisfied: pyparsing>=2.2.1 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from matplotlib) (3.0.9)\n",
      "Requirement already satisfied: basemap-data<1.4,>=1.3.2 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from basemap) (1.3.2)\n",
      "Requirement already satisfied: pyproj<3.7.0,>=1.9.3 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from basemap) (3.6.1)\n",
      "Requirement already satisfied: pyshp<2.4,>=1.2 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from basemap) (2.3.1)\n",
      "Requirement already satisfied: certifi in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from pyproj<3.7.0,>=1.9.3->basemap) (2022.9.24)\n",
      "Requirement already satisfied: six>=1.5 in /Users/mashi/opt/anaconda3/lib/python3.9/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "pip install matplotlib basemap\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b2613c08",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mashi/opt/anaconda3/lib/python3.9/site-packages/mpl_toolkits/mplot3d/proj3d.py:113: RuntimeWarning: invalid value encountered in divide\n",
      "  txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "# Create a figure and a 3D axis\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# Set up the map projection using Basemap\n",
    "m = Basemap(projection='ortho', lat_0=0, lon_0=0)\n",
    "\n",
    "# Get a grid of latitudes and longitudes\n",
    "lons = np.linspace(-180, 180, 100)\n",
    "lats = np.linspace(-90, 90, 100)\n",
    "lon_grid, lat_grid = np.meshgrid(lons, lats)\n",
    "\n",
    "# Convert to map projection coordinates\n",
    "x, y = m(lon_grid, lat_grid)\n",
    "\n",
    "# Generate a 3D surface plot\n",
    "z = np.sin(np.radians(lat_grid))  # Simple Z function\n",
    "\n",
    "ax.plot_surface(x, y, z, cmap='viridis')\n",
    "\n",
    "# Plot lat/lon points if needed\n",
    "lon_points, lat_points = np.array([10, 20, 30]), np.array([10, 15, 25])\n",
    "x_points, y_points = m(lon_points, lat_points)\n",
    "z_points = np.sin(np.radians(lat_points))\n",
    "ax.scatter(x_points, y_points, z_points, color='r', s=50)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6b702eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sum(a, b):\n",
    "    \"\"\"\n",
    "    >>> sum(4, 3)\n",
    "    7\n",
    "\n",
    "    >>> sum(-4, 5)\n",
    "    1\n",
    "    \"\"\"\n",
    "    return a + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d74de889",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<enumerate at 0x7fb392a70e00>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "college_years = ['Freshman', 'Sophomore', 'Junior', 'Senior']\n",
    "enumerate(college_years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "a42e2791",
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "invalid index to scalar variable.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/x7/39dnrb7j2qgck_5lhpy2d04w0000gn/T/ipykernel_45126/1062576365.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     59\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maligned_est\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maligned_est\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'red'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Estimated Coordinates (Aligned)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     60\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlon\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maligned_est\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 61\u001b[0;31m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlon\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlon\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34mf\"Est {i+1}\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m12\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'red'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     62\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     63\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Longitude'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: invalid index to scalar variable."
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.spatial import procrustes\n",
    "\n",
    "# Function to compute the distance between two points (lat, lon)\n",
    "def haversine(lat1, lon1, lat2, lon2):\n",
    "    # Radius of the Earth in kilometers (mean radius)\n",
    "    R = 6371.0\n",
    "    lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])\n",
    "\n",
    "    dlat = lat2 - lat1\n",
    "    dlon = lon2 - lon1\n",
    "\n",
    "    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2\n",
    "    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))\n",
    "\n",
    "    distance = R * c\n",
    "    return distance\n",
    "\n",
    "# Example latitudes and longitudes (original coordinates)\n",
    "latitudes = [10, 20, 30]\n",
    "longitudes = [50, 60, 70]\n",
    "\n",
    "num_cities = len(latitudes)\n",
    "\n",
    "# Simulated estimated latitudes and longitudes (not original but preserving distances)\n",
    "estimated_latitudes = [12, 22, 33]\n",
    "estimated_longitudes = [52, 61, 71]\n",
    "\n",
    "# Step 1: Ensure distances are preserved\n",
    "# Compute original distances\n",
    "original_distances = np.zeros((num_cities, num_cities))\n",
    "for i in range(num_cities):\n",
    "    for j in range(i + 1, num_cities):\n",
    "        original_distances[i, j] = haversine(latitudes[i], longitudes[i], latitudes[j], longitudes[j])\n",
    "\n",
    "# Compute estimated distances\n",
    "estimated_distances = np.zeros((num_cities, num_cities))\n",
    "for i in range(num_cities):\n",
    "    for j in range(i + 1, num_cities):\n",
    "        estimated_distances[i, j] = haversine(estimated_latitudes[i], estimated_longitudes[i], estimated_latitudes[j], estimated_longitudes[j])\n",
    "\n",
    "# Step 2: Align the estimated points with the original ones using Procrustes analysis\n",
    "original_coords = np.array([latitudes, longitudes]).T\n",
    "estimated_coords = np.array([estimated_latitudes, estimated_longitudes]).T\n",
    "\n",
    "# Procrustes analysis for optimal rotation and translation in 2D\n",
    "_, aligned_est, _ = procrustes(original_coords, estimated_coords)\n",
    "\n",
    "# Step 3: Plot the original and estimated (aligned) points\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "# Plot original points\n",
    "#plt.scatter(original_coords[:, 1], original_coords[:, 0], color='blue', label='Original Coordinates')\n",
    "#for i, (lat, lon) in enumerate(zip(latitudes, longitudes)):\n",
    "#    plt.text(lon, lat, f\"City {i+1}\", fontsize=12, color='blue')\n",
    "\n",
    "# Plot estimated points (after alignment)\n",
    "plt.scatter(aligned_est[:, 1], aligned_est[:, 0], color='red', label='Estimated Coordinates (Aligned)')\n",
    "for i, (lat, lon) in enumerate(aligned_est):\n",
    "    plt.text(lon[1], lon[0], f\"Est {i+1}\", fontsize=12, color='red')\n",
    "\n",
    "plt.xlabel('Longitude')\n",
    "plt.ylabel('Latitude')\n",
    "plt.title('Original vs Estimated (Aligned) Coordinates')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9434196e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}