1 line
No EOL
8 KiB
Text
1 line
No EOL
8 KiB
Text
{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"fractaldimension.ipynb","provenance":[],"authorship_tag":"ABX9TyPDGTZCMhlplMpdJFPM0/7i"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"code","metadata":{"id":"uo4kxbshsfPU","executionInfo":{"status":"ok","timestamp":1622567114654,"user_tz":240,"elapsed":118,"user":{"displayName":"Dylan Torrance","photoUrl":"","userId":"06723129446585093176"}}},"source":["import random\n","import matplotlib.pyplot as plt\n","from scipy.optimize import curve_fit\n","import numpy as np\n","import math"],"execution_count":27,"outputs":[]},{"cell_type":"code","metadata":{"id":"G01gxmY2QI64"},"source":["####This is poor generator stuff, you can ignore \n","\n","shrink = np.array([[1,0,0,0,0,0,0,0],\n"," [0,0,1,0,0,0,0,0],\n"," [0,0,0,0,0,1,0,0],\n"," [0,0,0,0,0,0,0,1]])\n","q = np.linalg.inv(np.array([[1,-3,-3,-3],\n"," [-3,1,-3,-3],\n"," [-3,-3,1,-3],\n"," [-3,-3,-3,1]]))\n","\n","#take transpose of row reduced gram matrix to get linear relations \n","#pivot columns work \n","\n","linrel = 1/2 * np.array([[2,0,0,0],\n"," [1,1,-1,1],\n"," [0,2,0,0],\n"," [1,1,1,-1],\n"," [-1,1,1,1],\n"," [0,0,2,0],\n"," [1,-1,1,1],\n"," [0,0,0,2]])\n","alpha = np.array([2*math.sqrt(2),2*math.sqrt(2),2*math.sqrt(2),2*math.sqrt(2),0,0,0,0])\n","\n","wtpw = np.matmul(linrel,np.matmul(q,shrink))\n","i = np.identity(8)\n","\n","curv = np.array([-1,2,7,4,11,8,3,6])\n","\n","print(i - 2* np.matmul(np.outer(alpha,alpha),wtpw)/abs(np.matmul(alpha,np.matmul(wtpw,alpha))))\n","\n","print(np.matmul((i - 2*np.matmul(np.outer(alpha,alpha),wtpw)/abs(np.matmul(alpha,np.matmul(wtpw,alpha)))), curv))\n","\n"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"xY_oNcRDl797"},"source":["##This is just a demonstration of how it might walk over the words, can ignore unless you want to see that, not used in the actual calculation.\n","\n","word = [1]\n","while word:\n"," print(word)\n"," if random.uniform(0, len(word))<1:\n"," word = over(word)\n"," else:\n"," word = down(word)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"5GPpjyI9ivpl","executionInfo":{"status":"ok","timestamp":1622567124989,"user_tz":240,"elapsed":490,"user":{"displayName":"Dylan Torrance","photoUrl":"","userId":"06723129446585093176"}}},"source":["##These are the main functions used for \n","\n","\n","#function that returns the new curvature for the nth generator, just for checking whether we should go deeper into the word tree or start going down \n","def transval(n,root):\n"," return 2*sum(root)-3*root[n-1]\n","\n","#function for when the word would have given a curvature that is too big \n","def down(word): \n"," #if not word:\n"," #return word #this is older stuff you can probably ignore\n","\n"," if word[-1] == 4: #if the word ends in the last generator, drop that and repeat this whole function with that new word\n"," word.pop()\n"," word = down(word)\n"," elif len(word)>1 and word[-1]+1==word[-2]: #if the word has length 2 and the last two generators are n+1, n\n"," if word[-1]==3: #if the word ended in \"43\", then drop those two (since that is as far as we could go in this part of the tree), and if it is nonempty, then call this whole function again \n"," word.pop()\n"," word.pop()\n"," if word: \n"," word = down(word)\n"," else: #if it ended in somehting like \"21\", then that becomes \"23\", like it skips over \"22,\" since that wouldn't be reduced.\n"," word[-1]=word[-1]+2\n"," else: #if it wasn't one of the weird above cases then just change the last generator to the one after that \n"," word[-1]=word[-1]+1\n"," return word\n","\n","#function for when the word gives a curvature that's still under the bound\n","#just adds a 1 or 2 at the end of the word, which is basically like the left-most branch above the word in the tree\n","def over(word):\n"," if word[-1] == 1:\n"," word.append(2)\n"," else:\n"," word.append(1)\n"," return word"],"execution_count":30,"outputs":[]},{"cell_type":"code","metadata":{"id":"-_JTSiM3cJqf","executionInfo":{"status":"ok","timestamp":1622567127105,"user_tz":240,"elapsed":3,"user":{"displayName":"Dylan Torrance","photoUrl":"","userId":"06723129446585093176"}}},"source":["#this is the function we want to fit the curve to\n","def func(x, a, b):\n"," return a * x**b\n","\n","#this function fits the curve and doesn't actually return anything but prints the fractal dimension\n","#takes in a number of how many times we want to divide , the bound from below, and the list of curvatures under that bound\n","def curvefit(n, bound, curvatures):\n"," x = []\n"," y = []\n"," for j in range(n): #n is just how many sub-bounds we want to have, so if our bound was like curvatures under 30,000, and we had n = 30\n"," #then we would be counting like curvatures under 1000, 2000, 3000, ... 30,0000. \n"," x.append(int((j+1)*bound/n))\n"," y.append(sum(i < (j+1)*bound/n for i in circlelist)) #this is just counting the number of curvatures under each of our sub-bounds \n"," popt,pcov=curve_fit(func,x,y) #curve fit to our list of like [1000,2000,...,30000] matched with the count of circles under each of those bounds\n"," print(\"fractal dimension: \" + str(popt[1])) #popt is like the [a,b] that is the best curve fit for ax^b, so we just take b "],"execution_count":31,"outputs":[]},{"cell_type":"code","metadata":{"id":"9_nsjk3fklRW","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1622567131441,"user_tz":240,"elapsed":2442,"user":{"displayName":"Dylan Torrance","photoUrl":"","userId":"06723129446585093176"}},"outputId":"2dd3d1c8-10f9-4965-aaa5-8bf49513eaa2"},"source":["root = [-10, 18, 23, 27] #some base tuple, it will change as we go through the words since we are kinda only paying attention to the end of the word\n","circlelist = list(root) #list of the curvatures, it starts with the root stuff and we will add to it \n","word = [1] #word that we are currently considering, it will change as we go through these, the first word we start with is just the word that consists of the first generator\n","recentroot = [] #this will be an array of arrays. the recentroot[i] is the most recent tuple that wasn't too big for word of length i so that when we go back down we know which tuple to operate on\n","recentroot.append(root) #the root thing will be our 0th element of this array\n","bound = 100000 #bound for what we want our curvatures to be under\n","\n","while word: #this runs until you get back to the empty word\n"," #print(\"root : \" + str(root))\n"," #print(\"word : \" + str(word))\n"," n = word[-1] #we only really pay attention to the last generator in our word\n"," l = len(word) #so we know which tuple to pull from recentroot\n"," #print(recentroot)\n"," root = list(recentroot[l-1]) #you have to cast it as a list or else it is just referencing an element of recentroot and will change when that changes, which we don't want\n"," if transval(n,root)>bound: #if the new curvature was bigger than our bound, change the word according to the down function \n"," word = down(word)\n"," else: #if it was not bigger than our bound, replace/append the corresponding element of recentroot, add that curvature to the circlelist, \n"," #and change the word according to the over function, since we may go deeper into that part of the tree before we hit the bound\n"," root[n-1] = transval(n,root)\n"," if l>=len(recentroot):\n"," recentroot.append(root)\n"," else:\n"," recentroot[l]=root\n"," circlelist.append(root[n-1])\n"," word = over(word)\n","\n","#print(circlelist)\n","print(len(circlelist))\n","curvefit(50,bound,circlelist) #this will print our fractal dimension once the circlelist has been made"],"execution_count":32,"outputs":[{"output_type":"stream","text":["66991\n","fractal dimension: 1.3049498628054264\n"],"name":"stdout"}]}]} |