# -*- coding: cp1252 -*- def diameter_at_height(hx, treeheight, d13, b): """This function returns the diameter at height hx""" # check the sanity of hx, the height of interest if (hx < 0) or (hx > treeheight): print "diameter_at_height(), problem with hx!" return 0 # Do conversions to float to assure "correct" decimal math hx = float(hx); d13 = float(d13); treeheight = float(treeheight) # Solve the diameter at 20% height, using the d13 as reference x = (1-1.3/treeheight) fx = b[0]*x + b[1]*x**2 + b[2]*x**3 + b[3]*x**5 + b[4]*x**8+ b[5]*x**13 + b[6]*x**21 + b[7]*x**34 d2h=d13/fx # Now, solve for diameter dx at height hx x = (1-hx/treeheight) fx = b[0]*x + b[1]*x**2 + b[2]*x**3 + b[3]*x**5 + b[4]*x**8+ b[5]*x**13 + b[6]*x**21 + b[7]*x**34 return d2h*fx def height_for_diameter(dx, treeheight, d13, b): """Finds thru iteration height at which diameter is dx""" # Data type changes to float and initialization of needed vars h = float(treeheight); g = 1.0; step = 1.0 d13 = float(d13) dx = float(dx) # Check the sanity of the parameters if ((0 > dx) or (dx > 85.0) or (1.4 > h) or (h > 40.0)): print "d13:" , dx, " h:" , h, " ? arguments in height for diameter" return 0 # Do a coarse 1-m step-search down from the top and stop when diameter # is greater than dx l = h / step # number of coarse steps hr = h # hr stores the current height if (l * step == h): hr = hr + 0.01 # start point for iteration slightly above top # Coarse search, down from the top, by 1 m step for i in range (1, int(l+1),1): hr = hr - step a = diameter_at_height(hr, h, d13, b) # if the diameter is greater than dx, stop coarse search if (a) > dx: break # Fine tuning, coarse solution (height) now stored in variable hr g = 1.0 # direction to go to has to be reversed, must return up the stem diff = 1 # assign value for accuracy variable while (abs(diff)>0.001) : # The wanted accuracy is diff, we continue until this acc. is reached step = step / 2.0 # reduce the step size at every round of loop hr = hr + g * step # compute new candidate position in the direction of g (-1 or 1) dr = diameter_at_height(hr, h, d13, b) # check diameter for hr diff = dx-dr # difference if (dr > dx): # check if we still are approching from below g = 1.0 if (dr <= dx): # if not, bo back down from above g = -1.0 return hr def log_volume(h1, h2, h, d13, b): """Computes the volume between heights h_low and h_high""" if (h1 > h2): (h1,h2) = (h2, h1) # swap the order if (h1 < 0 or h2 > h): print "Parameter problem in log_volume" return -1 B1 = b[0]; B2 = b[1]; B3 = b[2]; B4 = b[3] B5 = b[4]; B6 = b[5]; B7 = b[6]; B8 = b[7] x = (1-float(h1)/h) # Evaluate integral at low end V1 = 2./37*B2*B8*x**37+1./20*B8*B4*x**40+2./23*B1*B7*x**23+1./24*B6*B8*x**48+ \ 1./12*B7*B2*x**24+2./15*B1*B6*x**15+1./4*B4*B2*x**8+1./7*B5*B4*x**14+ \ 1./2*B1*B2*x**4+1./5*B5*B1*x**10+1./19*B3*B8*x**38+2.0/19*B6*B4*x**19+2./9*B3*B4*x**9 \ +1./3*B2*B3*x**6+1./28*B7*B8*x**56+1./8*B2*B6*x**16+2./25*B7*B3*x**25+1./18*B8*B1*x**36 \ +1./15*B7*B5*x**30+1./11*B6*B5*x**22+1./6*B3*B5*x**12+2./35*B6*B7*x**35+1./27*(B6**2+2*B4*B7)*x**27 \ +1./43*(2*B8*B5+B7**2)*x**43+1./7*(2*B4*B1+B3**2)*x**7+1./5*(2*B1*B3+B2**2)*x**5+1./69*B8**2*x**69 \ +1./17*(B5**2+2*B6*B3)*x**17+1./11*(B4**2+2*B5*B2)*x**11+1./3*B1**2*x**3 x = (1-float(h2)/h) # Evaluate integral at high end V2 = 2./37*B2*B8*x**37+1./20*B8*B4*x**40+2./23*B1*B7*x**23+1./24*B6*B8*x**48+ \ 1./12*B7*B2*x**24+2./15*B1*B6*x**15+1./4*B4*B2*x**8+1./7*B5*B4*x**14+ \ 1./2*B1*B2*x**4+1./5*B5*B1*x**10+1./19*B3*B8*x**38+2./19*B6*B4*x**19+2./9*B3*B4*x**9 \ +1./3*B2*B3*x**6+1./28*B7*B8*x**56+1./8*B2*B6*x**16+2./25*B7*B3*x**25+1./18*B8*B1*x**36 \ +1./15*B7*B5*x**30+1./11*B6*B5*x**22+1./6*B3*B5*x**12+2./35*B6*B7*x**35+1./27*(B6**2+2*B4*B7)*x**27 \ +1./43*(2*B8*B5+B7**2)*x**43+1./7*(2*B4*B1+B3**2)*x**7+1./5*(2*B1*B3+B2**2)*x**5+1./69*B8**2*x**69 \ +1./17*(B5**2+2*B6*B3)*x**17+1./11*(B4**2+2*B5*B2)*x**11+1./3*B1**2*x**3 x = (1-1.3/h) fx = b[0]*x + b[1]*x**2 + b[2]*x**3 + b[3]*x**5 + b[4]*x**8+ b[5]*x**13 + b[6]*x**21 + b[7]*x**34 d2h=d13/fx return (V1-V2) * (d2h/100.0)**2 * h * math.pi/4. class Tree: def __init__(self): self.sp = None self.d13 = None self.h = None self.hc_dry = None self.hc_living = None self.b=[] def get_user_input(tree): tree.sp = 1 # input("Give me species, 1 = pine: , 2= spruce, 3= rest") tree.d13 = 33 # input("dbh, d13 in cm: ") tree.h = 27 # input("height, h in m: ") tree.hc_dry = 9 # input("Height of the lowest dry branch, m: ") tree.hc_living = 17 # input("Height of the lowest living branch, m:") return tree # 5. Assign taper curve's coefficients to b[], pine, spruce, birch def gimme_b(sp, b): if sp == 1: b.append(2.1288),b.append(-0.63157),b.append(-1.6082),b.append(2.4886) b.append(-2.4147),b.append(2.3619), b.append(-1.7539),b.append(1.0817) if sp == 2: b.append(2.3366), b.append(-3.2684) ,b.append(3.6513), b.append(-2.2608) b.append(0), b.append(2.1501),b.apend(-2.7412),b.append(1.8876) if sp == 3: b.append(0.93838), b.append(4.106), b.append(-7.8517), b.append(7.8993) b.append(-7.5018), b.append(6.3863), b.append(-4.3918), b.append(2.1604) return b def check_this_log(letter, h_low, h_high, t, d_lengths, d_prices, d_mindim): # This instructor checks if the log is valid, and assigns value == 0 # if not, otherwise assigns and returns value in €. d_top = diameter_at_height(h_high, t.h, t.d13, t.b) if letter == "A" or letter == "B" or letter == "C" or letter == "D" : # its a saw log if d_top > d_mindim["Saw_1"]: # Check if the diameter criteria is met if t.hc_dry > h_high: # No branches in this log value= d_prices["Saw_1"] * log_volume(h_low, h_high, t.h, t.d13, t.b) return value, "Kossu" elif t.hc_dry < h_high and t.hc_living > h_high: # Just dry branches? value= d_prices["Saw_2"] * log_volume(h_low, h_high, t.h, t.d13, t.b) return value, "Jallu" elif t.hc_living < h_high: # lowest price category value= d_prices["Saw_3"] * log_volume(h_low, h_high, t.h, t.d13, t.b) return value, "Viski" else: return 0, "0" if letter == "E" : # its a pulp wood log if d_top > d_mindim["Pulp"]: value= d_prices["Pulp"] * log_volume(h_low, h_high, t.h, t.d13, t.b) return value, "Keppana" else: return 0, "0" if letter == "F" : # its a Fuel wood log if d_top > d_mindim["Fuel"]: value= d_prices["Fuel"] * log_volume(h_low, h_high, t.h, t.d13, t.b) return value, "Vissy" else: return 0, "0" return 0, "0" def get_the_word(i, M): # i to binary string Bs = bin(i) # remove the '0b' prefix Bs = Bs[2:len(Bs)] # It needs to have 2*M characters, padd from the left with zeroes L_padd = 2*M - len(Bs) # This many zeroes are needed if L_padd > 0: for c1 in range(L_padd): Bs = "0" + Bs # Padd from the left L_padd times # Just check that the length now is ok if len(Bs) != 2*M : print "problem in making the word!" # Now start the conversion into base-4, pairwise manner in the binary string ABCD = "ABCD" word = "" # This will be outputted finally having the M-character 4-letter word for i in range(len(Bs)-1, 0, -2): # Get the pairs starting from the end B4num = int(Bs[i-1]+Bs[i], 2) # Convert a binary base-2 string into integer, 0,1,2, or 3 word = word + ABCD[B4num] # Concatenate the corresponding letter 0=A, 1=B, 2=C, 3=D return word import random def get_the_word_thru_lottery(i, M): # This one does Monte Carlo style lottery for the M-character word word = "" for i in range(M): rand = random.random() if rand <= 1./6. : word += "A" if rand > 1./6. and rand <= 2./6. : word += "B" if rand > 2./6. and rand <= 3./6. : word += "C" if rand > 3./6. and rand <= 4./6. : word += "D" if rand > 4./6. and rand <= 5./6.: word += "E" if rand > 5./6. : word += "F" return word # EXERCISE 1; MAIN PROGRAM STARS HERE *************** import math # required for pi, in math.pi b=[] # list that holds the coefficients of the taper curve t=Tree() # 1. The data structures that hold the rules: # ************** NEW ************************ d_lengths = {'A': 4.0, 'B': 4.6, 'C': 5.2, 'D': 6.1, "E": 3.0, "F" : 2.5} # ************** four saw log lenghts now *** d_prices = {'Saw_1':50.0,'Saw_2':47.0,'Saw_3':38.0,'Pulp':28.0,'Fuel':40.0} d_mindim = {'Saw_1':18.0,'Saw_2':18.0,'Saw_3':18.0,'Pulp': 9.0,'Fuel': 4.0} # 2. The minimum diameter of any assortment minimum_diameter = min(d_mindim.values()) # 3. The minimum length of any log minimum_length = min(d_lengths.values()) # 4. Get the tree attributes, force sp to be pine, ==1 t = get_user_input(t) # 5. retrive the coeffiecients t.b = gimme_b(t.sp, b) # 6. Find the upper end of merchantable timber, where minimum diameter occurs h_max = height_for_diameter(minimum_diameter, t.h, t.d13, t.b) # 7. How many logs at maximum (worst case) h_stump = 0.1 length_merc = h_max - h_stump # M is the number of logs (length of bucking order-words) M = int(length_merc/minimum_length) # Main loop that is executed 4**M times (altertives, full enumeration) i = 0; current_maximum = -1 percentage = 4 # ************ NEW ********** while(i < 6**M * percentage/100.0): # ********** NEW ************ # Felling of the tree h_low = h_stump # Next, we need the bucking instructions, M-character, 4-letter word # word = get_the_word(i, M) # hiphurraa! # *********** NEW ******************** word = get_the_word_thru_lottery(i, M) # MONTE-CARLO STYLE #************************************* tot_val = 0; tot_vol = 0 check_word= ""; quality_word = "" for letter in word: # Loop for cutting the logs, the height of the upper end h_high = h_low + d_lengths[letter] if h_high < h_max: # Do a log! value, Quality = check_this_log(letter, h_low, h_high, t, d_lengths, d_prices, d_mindim) if value == 0: break check_word = check_word + letter quality_word = quality_word + Quality + "," vol = log_volume(h_low, h_high, t.h, t.d13, t.b) tot_vol += vol tot_val += value h_low = h_high if tot_val > current_maximum : print "Val %5.2f € Vol %5.3f m3 %s %s" % (tot_val, tot_vol, check_word, quality_word) current_maximum = tot_val i += 1 exit () #for htest in range(0,int(h),1): # print htest, diameter_at_height(htest, h, d13, b)