{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Bifurcation Diagram\n", "\n", "In this notebook, we generate the bifurcation diagram." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib qt\n", "import numpy as np\n", "from TriangleFunctions import *\n", "import math\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start by looking at triangle code, and determining method to go about it. I will need to come up with a way to keep iterating until get a point that matches, and then graph all the previous points until the last one that matches. \n", "\n", "That is, take all the previous points up until the last one that matches, and compute over them to determine their y value on the graph. The given set of points will then be paired with an alpha value. Can geometrically find their y value useing the method in my notebook. \n", "\n", "Will also need to start each one on both sides to ensure get all possible points.\n", "\n", "## Data Setup" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def return_iso_points(angle1, right_point):\n", " \"\"\"anlge1 is top angle, angle2 is right angle.\n", " The left point is on origin, the right point is defined by you, the top point is found computationally\"\"\"\n", " \n", " #find left angle\n", " angle2 = ((180 - angle1) / 2)\n", " angle3 = angle2\n", " \n", " #define left line, slope and intersection. Use tan to define slope\n", " m1 = math.tan(math.radians(angle3))\n", " b1 = 0\n", " mb1 = [m1, b1]\n", " \n", " #define right line\n", " m2 = -(math.tan(math.radians(angle2)))\n", " b2 = newB(m2, right_point)\n", " mb2 = [m2, b2]\n", " \n", " # find out when the line defined by left and right side intersect to get your final point.\n", " pointT = inter(mb2, mb1)\n", " \n", " return (0, 0), pointT, right_point" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "## want to edit this function. Add into method the while loop, and iteration step value of the angle.\n", "# need to use return point and feed as input into triangle so as to get initial conditions of the triangle\n", "def orbit_points(startPoint, points):\n", " \"\"\"\n", " startPoint: the point to start the process from within the triangle\n", " points: the set of points that define the outer triangle\n", " \"\"\"\n", " \n", " # current point to draw from\n", " current_point = startPoint\n", " # track visited points for graphing\n", " visited_points = []\n", " indx = 0\n", " \n", " # define lines of triangle. Of form [m, b]\n", " # bottom line is on x axis so don't have to draw it.\n", " l1 = findFunc(points[0], points[1])\n", " l2 = findFunc(points[1], points[2])\n", " \n", " # define perpindicular functions. Returns m'\n", " lp1 = pFunc(l1[0])\n", " lp2 = pFunc(l2[0])\n", " \n", " # iterate until point is in loop. When point is in set, get its index\n", " while True:\n", " \n", " if(current_point in visited_points):\n", " # get indx\n", " indx = visited_points.index(current_point)\n", " break\n", " \n", " \n", " visited_points.append(current_point)\n", " \n", " # for each perpendicular line, define the exact function that passes through the current point to find next step\n", " lOne = [lp1, newB(lp1, current_point)]\n", " inter1 = inter(l1, lOne)\n", " \n", " lTwo = [lp2, newB(lp2, current_point)]\n", " inter2 = inter(l2, lTwo)\n", " \n", " # this corresponds to the botton line which lies along the x axis so need a different function rule for it. \n", " inter3 = [current_point[0], 0]\n", " \n", " dist1 = dist(current_point, inter1)\n", " dist2 = dist(current_point, inter2)\n", " #dist3 = dist(current_point, inter3)\n", " ## NEW CODE\n", " dist3 = current_point[1]\n", " \n", " # use largest to pick the longest distance to travel within triangle\n", " largest = heapq.nlargest(3, [dist1, dist2, dist3])\n", " cur_largest = largest[0]\n", " next_indx = largest.index(cur_largest) + 1\n", " next_largest = largest[next_indx]\n", " \n", " \"\"\"### NEW CODE UPDATE, AFTER BIF 2, OLD IN RED COMMENTS\n", " if(cur_largest == dist1):\n", " current_point = inter1\n", " \n", " if(cur_largest == dist2):\n", " current_point = inter2\n", " \n", " if(cur_largest == dist3):\n", " current_point = inter3\"\"\"\n", " \n", " \n", " # prevent program from drawing outside of triangle\n", " \n", " while(True):\n", " if((cur_largest == dist1) and (inter1[1] >= points[1][1])):\n", " cur_largest = largest[next_indx]\n", " if((cur_largest == dist2) and (inter2[1] >= points[1][1])):\n", " cur_largest = largest[next_indx]\n", " if((cur_largest == dist3) and (inter3[1] >= points[1][1])):\n", " cur_largest = largest[next_indx]\n", " else:\n", " break\n", " \n", " # if two distances are the same pick a random one to follow\n", " if(cur_largest == next_largest):\n", " # make a random choice between the two\n", " largest_pick = random.choice([cur_largest, next_largest])\n", "\n", " if(largest_pick == dist1):\n", " current_point = inter1\n", " #visited_points.append(current_point)\n", "\n", " if(largest_pick == dist2):\n", " current_point = inter2\n", " #visited_points.append(current_point)\n", "\n", " if(largest_pick == dist3):\n", " current_point = inter3\n", " #visited_points.append(current_point)\n", " \n", " # if all different distances and within rules, move to the new point\n", " if(cur_largest != next_largest): \n", " \n", " if(cur_largest == dist1):\n", " current_point = inter1\n", " #visited_points.append(current_point)\n", " if(cur_largest == dist2):\n", " current_point = inter2\n", " #visited_points.append(current_point)\n", " if(cur_largest == dist3):\n", " current_point = inter3\n", " #visited_points.append(current_point)\n", " \n", " \n", " \n", " return visited_points, indx" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "# test\n", "triangle = return_iso_points(2, [1, 0])\n", "pnts, indx = orbit_points([0, 0], triangle)\n", "#print(pnts), print(indx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate Data" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.07427972123566462" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.random.uniform(low=0, high=(1/3))" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# setup for loop to get the set\n", "def gen_data(angle_start, angle_step_size):\n", " points = []\n", " indices = []\n", " x_angle = []\n", " \n", " num_steps = round((angle_start/angle_step_size))\n", " \n", " curr_angle = angle_start\n", " \n", " for i in range(0, num_steps):\n", " \n", " x_angle.append(curr_angle)\n", " \n", " # ADD RANDOM POINT\n", " rand_pnt = np.random.uniform(low=0, high=(1/3))\n", " \n", " triangle = return_iso_points(curr_angle, [(1/3),0])\n", " \n", " # START IN MIDDLE TO MAKE UPPER LEVELS WORK\n", " #pnts, indx = orbit_points([((1/6)-.001), 0], triangle)\n", " \n", " ## NEW RANDOM POINT TEST\n", " pnts, indx = orbit_points([rand_pnt, 0], triangle)\n", " \n", " points.append(pnts)\n", " indices.append(indx)\n", " \n", " curr_angle = curr_angle - angle_step_size \n", " \n", " return points, indices, x_angle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Clean Data\n", "Begin by recovering the wanted data points from each set" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# get points we care about\n", "def gen_pure_data(points, indices):\n", " pure_points = []\n", " for i in range(0, len(indices)):\n", " curr_indx = indices[i]\n", " curr_point = points[i]\n", " \n", " pp = curr_point[curr_indx:]\n", " pure_points.append(pp)\n", " return pure_points" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# method get_point_dist will give the edge value of the given point.\n", "def get_y_dist(angle, point):\n", " #math.sin(x) takes radians, angle comes in degrees\n", " # math.radians(x) converts from degrees to radians\n", " \n", " # construct angle off of alpha, assume angle is alpha\n", " angle = (180 - angle) / 2\n", " \n", " # first get y value of point, as that is what matters:\n", " h = point[1]\n", " angle = math.radians(angle)\n", " dist = h / (math.sin(angle))\n", " return dist" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "# for each pure point get distances, will have to simultaneously determine side and add length\n", "# KEEP IN MIND THE INTERVAL LENGTH WE ARE DOING IS 1, SO AXIS IS 3\n", "def pure_dists(pure_points, angle_set):\n", " scaled_dists = []\n", " for i in range(0, len(pure_points)):\n", " curr_pnt = pure_points[i]\n", " sub_set = []\n", " for j in range(0, len(curr_pnt)):\n", " curr_sub_pnt = curr_pnt[j]\n", " # determine if x or not\n", " if(curr_sub_pnt[1] == 0):\n", " pnt = (1/3) + curr_sub_pnt[0]\n", " sub_set.append(pnt)\n", " continue\n", " \n", " # determine side of axis of point\n", " # left axis\n", " # WAS .5\n", " if(curr_sub_pnt[0] < (1/6)):\n", " dist = get_y_dist(angle_set[i], curr_sub_pnt)\n", " pnt = (1/3) - dist\n", " sub_set.append(pnt)\n", " # if on right axis \n", " if(curr_sub_pnt[0] > (1/6)):\n", " dist = get_y_dist(angle_set[i], curr_sub_pnt)\n", " pnt = (2/3) + dist\n", " sub_set.append(pnt)\n", " \n", " scaled_dists.append(sub_set)\n", " \n", " return scaled_dists" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# reflect points across axis. Axis is 1.5.\n", "# make into set at end to ensure unique points\n", "def mirror_points(scaled_dists):\n", " \n", " axis = .5\n", " mirrored_dists = []\n", " \n", " for i in scaled_dists:\n", " sub_mirrors = []\n", " for j in i:\n", " # first append point, then look at opposite\n", " sub_mirrors.append(j)\n", " \n", " d = abs((axis - j))\n", " if (j >= axis):\n", " nd = axis - d\n", " sub_mirrors.append(nd)\n", " \"\"\"if(nd in sub_mirrors):\n", " continue\n", " else:\n", " sub_mirrors.append(nd)\"\"\"\n", " if (j < axis):\n", " nd = axis + d\n", " sub_mirrors.append(nd)\n", " \"\"\"if(nd in sub_mirrors):\n", " continue\n", " else:\n", " sub_mirrors.append(nd)\"\"\"\n", " \n", " mirrored_dists.append(sub_mirrors) \n", " \n", " return mirrored_dists" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Data Demonstration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we make everything" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1min 11s, sys: 73.9 ms, total: 1min 12s\n", "Wall time: 1min 12s\n" ] } ], "source": [ "%%time\n", "# enter in angle to start at, and step size\n", "points, indices, x_angle = gen_data(150, 0.05)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "pure_points = gen_pure_data(points, indices)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "dists = pure_dists(pure_points, x_angle)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "final_point_set = mirror_points(dists)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For both possible orientations, mirrored on one another. To get 180 mark at end, append empy point to final set, and x_angle" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "final_point_set.append([])\n", "x_angle.append(180)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "# angles is x, final_point_set is y\n", "for xe, ye in zip(x_angle, final_point_set):\n", " plt.xticks(np.arange(min(x_angle)-.02, max(x_angle)+1, 30))\n", " plt.yticks(np.arange(0, 1.001, .5))\n", " plt.scatter([xe] * len(ye), ye, s=.04, color='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both orientations, no mirroring" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "for xe, ye in zip(x_angle, dists):\n", " plt.xticks(np.arange(min(x_angle)-.01, max(x_angle)+1, 90))\n", " plt.yticks(np.arange(0, 1.001, .5))\n", " plt.scatter([xe] * len(ye), ye, s=.09, color='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that copied bars are forming at 0 1, 2, 3.\n", "\n", "The thing giving the jumps is when the orientation is flipped. \n", "\n", "# flipping\n", "Mirror image just gives you flipped view, not actual diagram though\n", "All kinks worked out, now acts as expected although the same pretty mirror imagery stuff is not going on as in the one above :)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "def gen_data2(angle_start, angle_step_size):\n", " points = []\n", " indices = []\n", " points2 = []\n", " indices2 = []\n", " x_angle = []\n", " \n", " num_steps = round((angle_start/angle_step_size))\n", " \n", " curr_angle = angle_start\n", " \n", " for i in range(0, num_steps):\n", " \n", " x_angle.append(curr_angle)\n", " \n", " triangle = return_iso_points(curr_angle, [4,0])\n", " \n", " pnts, indx = orbit_points([(1/12), 0], triangle)\n", " pnts2, indx2 = orbit_points([(3/12), 0], triangle)\n", " \n", " points.append(pnts)\n", " indices.append(indx)\n", " points2.append(pnts)\n", " indices2.append(indx)\n", " \n", " curr_angle = curr_angle - angle_step_size \n", " \n", " return points, indices, points2, indices2, x_angle" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "#COPY\n", "def gen_pure_data2(points, indices, points2, indices2):\n", " pure_points = []\n", " for i in range(0, len(indices)):\n", " curr_indx = indices[i]\n", " curr_point = points[i]\n", " \n", " pp = curr_point[curr_indx:]\n", " \n", " curr_indx2 = indices2[i]\n", " curr_point2 = points2[i]\n", " \n", " pp2 = curr_point[curr_indx2:]\n", " \n", " ppp = pp + pp2\n", " \n", " pure_points.append(ppp)\n", " return pure_points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For just one orientation throughout" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "ename": "KeyboardInterrupt", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n", "\u001b[0;32m\u001b[0m in \u001b[0;36mgen_data2\u001b[0;34m(angle_start, angle_step_size)\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 18\u001b[0m \u001b[0mpnts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0morbit_points\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;36m12\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[0mtriangle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m \u001b[0mpnts2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindx2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0morbit_points\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m12\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[0mtriangle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 21\u001b[0m \u001b[0mpoints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpnts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m\u001b[0m in \u001b[0;36morbit_points\u001b[0;34m(startPoint, points)\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 26\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 27\u001b[0;31m \u001b[0;32mif\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcurrent_point\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mvisited_points\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 28\u001b[0m \u001b[0;31m# get indx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 29\u001b[0m \u001b[0mindx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvisited_points\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcurrent_point\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mKeyboardInterrupt\u001b[0m: " ] } ], "source": [ "%%time\n", "# enter in angle to start at, and step size\n", "points1, indices1, points2, indices2, x_angle = gen_data2(150, 0.1)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'points1' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mpure_points2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgen_pure_data2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoints1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindices1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpoints2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindices2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNameError\u001b[0m: name 'points1' is not defined" ] } ], "source": [ "pure_points2 = gen_pure_data2(points1, indices1, points2, indices2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dists2 = pure_dists(pure_points2, x_angle)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "for xe, ye in zip(x_angle, dists2):\n", " plt.xticks(np.arange(min(x_angle)-.1, max(x_angle)+1, 90))\n", " plt.yticks(np.arange(0, 1.001, .5))\n", " plt.scatter([xe] * len(ye), ye, s=.1, color='black')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }