Refactor ba_simlulate_sw

For #66 #48

Also revise factor finder for sb

This has introduced a discrepancy in the simulation results originating
in factor finder

The Sw factor is now 3.75 for plot 1049300, before it was 1.1352

Probably because of:

-    initial_age_sw = kwargs['startTageSw']
-    years_to_bh_sw = kwargs['y2bh_Sw']
-    species_comp_sw = kwargs['SC_Sw']
     site_index = kwargs['SI_bh_Sw']
-    present_density = kwargs['N_bh_SwT']
+    densities = kwargs['densities']

now, the species_composition is variable isntead of fixed at the time of
data

provided the time series of species_composition from the densities array
passes through the data in the appropriate year, this hsould not be
changed

it may not pass through it; and this was peculiar for aspen because
aspen density for young plots is more variable than other species

one option is to use a flag in the simulation functions to use constant
values for species_comp and present density (the values from the data year)

pygypsy/basal_area_factor.py

@@ -94,11 +94,8 @@ def estimate_basal_area_factor_sb(**kwargs):
     '''
     LOGGER.debug('Getting basal area factor for black spruce')
     initial_age = kwargs['startTage']
-    initial_age_sb = kwargs['startTageSb']
-    years_to_bh_sb = kwargs['y2bh_Sb']
-    species_comp_sb = kwargs['SC_Sb']
     site_index = kwargs['SI_bh_Sb']
-    present_density = kwargs['N_bh_SbT']
+    densities = kwargs['densities']
     density_at_bh_age = kwargs['N0_Sb']
     basal_area_at_bh_age = kwargs['BA_Sb0']
     present_basal_area = kwargs['BA_SbT']
@@ -109,10 +106,8 @@ def estimate_basal_area_factor_sb(**kwargs):
     iter_count = 0
 
     while not within_tolerance:
-        ba_est = sim_basal_area_sb(initial_age, initial_age_sb, years_to_bh_sb,
-                                   species_comp_sb, site_index, present_density,
-                                   density_at_bh_age, basal_area_at_bh_age,
-                                   factor)[-1]
+        ba_est = sim_basal_area_sb(initial_age, site_index, density_at_bh_age,
+                                   basal_area_at_bh_age, factor, densities)[-1]
 
         if abs(present_basal_area - ba_est) < tolerance:
             within_tolerance = True
@@ -163,11 +158,8 @@ def estimate_basal_area_factor_sw(**kwargs):
     '''
     LOGGER.debug('Getting basal area factor for white spruce')
     initial_age = kwargs['startTage']
-    initial_age_sw = kwargs['startTageSw']
-    years_to_bh_sw = kwargs['y2bh_Sw']
-    species_comp_sw = kwargs['SC_Sw']
     site_index = kwargs['SI_bh_Sw']
-    present_density = kwargs['N_bh_SwT']
+    densities = kwargs['densities']
     density_at_bh_age = kwargs['N0_Sw']
     sdf_aw = kwargs['SDF_Aw0']
     sdf_pl = kwargs['SDF_Pl0']
@@ -181,11 +173,10 @@ def estimate_basal_area_factor_sw(**kwargs):
     factor1 = 1.5* factor
 
     while not within_tolerance:
-        ba_est = sim_basal_area_sw(initial_age, initial_age_sw, years_to_bh_sw,
-                                   species_comp_sw, site_index, present_density,
+        ba_est = sim_basal_area_sw(initial_age, site_index,
                                    density_at_bh_age, sdf_aw, sdf_pl, sdf_sb,
                                    basal_area_at_bh_age,
-                                   factor)[-1]
+                                   factor, densities)[-1]
         if abs(present_basal_area - ba_est) < tolerance:
             within_tolerance = True
         else:

pygypsy/basal_area_simulate.py

@@ -2,6 +2,8 @@
 #pylint: disable=no-member
 # TODO: names - e.g. densities is actually densities and other varsiables
 # TODO: shouldn't bhage be a constant? is it actuall years since bh_age?
+# TODO: is it really density at bh age or is it density at time of data?
+
 import logging
 import numpy as np
 
@@ -130,24 +132,18 @@ def sim_basal_area_sb(initial_age, site_index, density_at_bh_age,
     return basal_area_arr
 
 
-def sim_basal_area_sw(initial_age, intial_age_sw, years_to_bh, spec_comp,
-                      site_index, present_density, density_at_bh_age, sdf_aw,
+def sim_basal_area_sw(initial_age, site_index, density_at_bh_age, sdf_aw,
                       sdf_pl, sdf_sb, basal_area_at_bh_age, correction_factor,
-                      use_correction_factor_future=False,
+                      densities, use_correction_factor_future=False,
                       stop_at_initial_age=True):
     '''Simlulate basal area forward in time for White Spruce
     It created the trajectory of basal area from bhage up to the inventory year
     given a correction factor that is being optimized
 
     :param float initial_age: Clock that uses the oldest species as a reference to
                               become the stand age
-    :param float intial_age: species specific age counted independently
-    :param float years_to_bh: time elapseed in years from zero to breast height
-                              age of sp Sw
-    :param float spec_comp: proportion of basal area for Sw
     :param float site_index: site index of species Sw
     :param float basal_area_at_bh_age: basal area of Sw at breast height age
-    :param float present_density: density of species Sw at time T
     :param float density_at_bh_age: initial density of species Sw at breast height age
     :param float sdf_pl: Stand Density Factor of species Pl
     :param float sdf_aw: Stand Density Factor of species Aw
@@ -166,19 +162,21 @@ def sim_basal_area_sw(initial_age, intial_age_sw, years_to_bh, spec_comp,
     basal_area_arr = np.zeros(max_age)
     basal_area_temp = basal_area_at_bh_age
 
-    while year < max_age:
+    for i, spec_comp_dict in enumerate(densities[0: max_age]):
         sc_factor = correction_factor \
                     if year < initial_age or use_correction_factor_future\
                     else 1
-        tage = intial_age_sw - initial_age
-        bh_age_sw = tage - years_to_bh
+        bh_age_sw = spec_comp_dict['bhage_Sw']
+        spec_proportion = spec_comp_dict['SC_Sw']
+        present_density = spec_comp_dict['N_bh_SwT']
 
         if density_at_bh_age > 0:
             if bh_age_sw > 0:
-                spec_comp = spec_comp * sc_factor
+                spec_proportion = spec_proportion * sc_factor
                 basal_area_increment = incr.increment_basal_area_sw(
-                    spec_comp, site_index, present_density, density_at_bh_age,
-                    bh_age_sw, sdf_aw, sdf_pl, sdf_sb, basal_area_temp
+                    spec_proportion, site_index, present_density,
+                    density_at_bh_age, bh_age_sw, sdf_aw, sdf_pl,
+                    sdf_sb, basal_area_temp
                 )
                 basal_area_temp = basal_area_temp + basal_area_increment
                 new_basal_area = basal_area_temp
@@ -191,10 +189,7 @@ def sim_basal_area_sw(initial_age, intial_age_sw, years_to_bh, spec_comp,
             basal_area_temp = 0
             new_basal_area = 0
 
-        basal_area_arr[year] = new_basal_area
-
-        year += 1
-        intial_age_sw += 1
+        basal_area_arr[i] = new_basal_area
 
     return basal_area_arr
 

pygypsy/forward_simulation.py

@@ -261,16 +261,16 @@ def simulate_forwards_df(plot_df, utiliz_params=None):
         #       increments that's something for another day; because it is
         #       complicated by the whole factor business
         BA_0_to_data_Aw_arr = sim_basal_area_aw(
-            startTage, SI_bh_Aw, N0_Aw, BA_Aw0, SDF_Aw0, f_Aw, densities,
+            startTage, SI_bh_Aw,  N0_Aw, BA_Aw0, SDF_Aw0, f_Aw, densities,
             use_correction_factor_future=True, stop_at_initial_age=False
         )
         BA_0_to_data_Sb_arr = sim_basal_area_sb(
             startTage, SI_bh_Sb, N0_Sb,
             BA_Sb0, f_Sb, densities, stop_at_initial_age=False
         )
         BA_0_to_data_Sw_arr = sim_basal_area_sw(
-            startTage, startTageSw, y2bh_Sw, SC_Sw, SI_bh_Sw, N_bh_SwT, N0_Sw,
-            SDF_Aw0, SDF_Pl0, SDF_Sb0, BA_Sw0, f_Sw, stop_at_initial_age=False
+            startTage, SI_bh_Sw, N0_Sw, SDF_Aw0, SDF_Pl0, SDF_Sb0, BA_Sw0,
+            f_Sw, densities, stop_at_initial_age=False
         )
         BA_0_to_data_Pl_arr = sim_basal_area_pl(
             startTage, startTagePl, y2bh_Pl, SC_Pl, SI_bh_Pl, N_bh_PlT, N0_Pl,