#!/usr/bin/env python
# encoding: utf-8
"""
Title: harvestop

Harves operation implementations

*harvestop.py*

Copyright (c) 2007 Simosol Oy.
Distributed under the GNU General Public License 2.0.
"""

import math
import numpy
from collections import defaultdict

from constants import STRATUMLABEL, EPSILON
from operations_ext import compute_treestem, bucking_algorithm
from carbonlifecycle import ProductLifeCycle, energy_conversion

class HarvestOperation(object):
    def __init__(self, harvestdata, cache_obj):
        super(HarvestOperation, self).__init__()
        self.hdata = harvestdata
        self.stem_cache = cache_obj
        if not self.stem_cache.initialized:
            divisor = self.hdata.cdata['divisor']
            log_reduction = int(self.hdata.cdata['log_reduction_type'])
            self.stem_cache.initialize(divisor, log_reduction)
        self.BA_current = self.hdata.get_ba_from_treelist()
        self.N_current = self.hdata.get_n_from_treelist()
        self.N_reserve_trees = 0.0
        self.target_N_reserve = self.hdata.cdata.get('n_reserve')
        self.VOL_reserve_trees = 0.0
        self.target_VOL_reserve = self.hdata.cdata.get('vol_reserve')
        self.stopnext = False
        self.N_seed_trees = 0.0
        self.target_N_seed = self.hdata.cdata.get('N_seedtrees')
        self.seedtree_sp = self.hdata.cdata.get('regen_sp')
        self.cut_timber = 0.0
        self.decimals = int(self.hdata.cdata['cache_resolution'])

    def cut_dclass_trees(self, dclassind, props, scale=None,
                         volres=False, biomass=True, trackchanges=True):
        """
        Remove proportion (prop) of stems in each diameter class.
        Estimate biomass, if not tracking changes, the diameter class
        modification and removal flags are not set in the data
        """
        # valueprop is the proportion of removed stems calculated in volume
        # and cash results, dataprop is to proportion removed from the data
        # the difference comes from the reserve trees
        valueprop, dataprop = props
        # make sure that the proportions are not negative (in that case, the
        # number of trees would actually increase!)
        valueprop = 0.0 if valueprop < 0.0 else valueprop
        dataprop = 0.0 if dataprop < 0.0 else dataprop
        strata = self.hdata.stratumdata
        ind = self.hdata.ind
        cdata = self.hdata.cdata
        targetvar = cdata['target_var']
        res = self.hdata.result
        nind = ind['tree_n']
        tree, istratum, species = self.get_dclass_info(dclassind)
        dbh = tree[ind['tree_d']]
        nstems = tree[nind]
        if nstems == 0.:
            return species, 0., 0., 0.
        nremoved = nstems * dataprop
        nstanding = nstems * (1.0 - dataprop)
        tree[nind] -= nremoved
        if tree[nind] < EPSILON or dataprop >= (1.0 - EPSILON):
            tree[nind] = 0.
            cdata['n_children'][istratum] -= 1
        stemvolume = res['dclass_vol'][dclassind]
        timbervolume = 0
        # Calculate the volumes of different timber assorments in the current
        # diameter class; reserve trees are not included in the volumes
        # eventhough they are removed from the data; hence the use of valueprop
        valuen = nstems * valueprop
        stimber = res['stratum_timber'][istratum]
        for assort, values in res['dclass_timber'][dclassind].items():
            stimber[assort][0] += values[0] * valuen
            if volres:
                stimber[assort][1] += values[0]
            else:
                stimber[assort][1] += values[1] * valuen
            timbervolume += values[0]
        scrapvol = stemvolume - timbervolume
        res['stratum_scrap'][istratum] += scrapvol * valuen
        if biomass:
            self._calc_biomass(istratum, tree, valuen, scrapvol)
        # set tracking variables
        cutba = math.pi * (dbh * 0.005) ** 2 * valuen
        self.BA_current -= cutba
        self.N_current -= valuen
        cutv_com = 0
        if targetvar == 'v_com':
            stratum_scale = scale[istratum]
            cutv_com = timbervolume * valuen * stratum_scale
            self.cut_timber += timbervolume * valuen * stratum_scale
        # set modified attribute values in original data
        # tag objects modified/deleted
        if trackchanges:
            if tree[nind] > EPSILON and tree[ind['tree_tag']] != 2:
                tree[ind['tree_tag']] = 1.0
            else:
                tree[ind['tree_tag']] = 2.0
            if 'stratum_N' in ind:
                strata[istratum, ind['stratum_N']] -= nremoved
                if strata[istratum, ind['stratum_N']] < EPSILON:
                    strata[istratum, ind['stratum_N']] = 0.0
            if 'stratum_ba' in ind:
                strata[istratum, ind['stratum_ba']] -= cutba
                if strata[istratum, ind['stratum_ba']] < EPSILON:
                    strata[istratum, ind['stratum_ba']] = 0.0
        return species, cutba, valuen, cutv_com

    def buck(self, dclassind):
        """
        Buck (buck) the trees and returns tree volume
        """
        ind = self.hdata.ind
        cdata = self.hdata.cdata
        res = self.hdata.result
        nind = ind['tree_n']
        tree, istratum, species = self.get_dclass_info(dclassind)
        dbh = tree[ind['tree_d']]
        height = tree[ind['tree_h']]
        nstems = tree[nind]
        ptkey = (cdata['price_region'], species)
        timber = res['dclass_timber'][dclassind]
        log_reduction = tree[ind['tree_log_reduction']]
        stemvol = self._buck_stem(species, dbh, height, nstems, ptkey,
                                  timber, istratum, log_reduction)
        res['dclass_vol'][dclassind] = stemvol

    def cut_evenly(self, scale=None):
        """
        Iterate through the sorted (ascending) treelist and remove a
        proportion relative to tree's expansion factor (N) on each iteration
        until target, minimum BA or minimum N is reached

        As a result, no trees (size classes) should get removed totally
        """
        removalorder = self.hdata.cdata['removal_order']
        target = self.hdata.cdata['target_var']
        ind = self.hdata.ind
        self.bottom_targetvar = 0.0
        self.top_targetvar = 0.0
        cdata = self.hdata.cdata
        d_limit = cdata['d_cut_all_above']
        min_class_size = cdata['min_class_size']
        removed_dc = set()
        self.sp_targetvar_taken = defaultdict(float)
        stop = False

        # Calculate the total removal proportion (either basal area or
        # expansion factor N) for an estimate on how large proportion should
        # be removed on each iteration
        if target == 'ba':
            removal_prop = cdata['to_remove'] / self.BA_current
        elif target == 'n':
            removal_prop = cdata['to_remove'] / self.N_current
        base_prop = removal_prop * 0.2
        while not stop:
            for i in range(removalorder.size):
                dcind = removalorder[i]
                tree, istratum, sp = self.get_dclass_info(dcind)

                # if size class diameter is greater that diameter limit or
                # expansion factor is less than minimum, remove the whole class
                d = tree[ind['tree_d']]
                n = tree[ind['tree_n']]
                if d-d_limit >= -EPSILON or n - min_class_size <= EPSILON:
                    prop = 1.0
                else:
                    prop = base_prop

                # check if removing the current proportion goes past the target
                if target == 'ba':
                    ba = self._get_dclass_ba(tree, ind)
                    stop = (self.BA_current - ba * prop) - cdata['target'] \
                                                                     <= EPSILON
                else:
                    n = tree[ind['tree_n']]
                    stop = (self.N_current - n * prop) - cdata['target'] \
                                                                     <= EPSILON
                if stop:
                    break
                if dcind in removed_dc:
                    continue

                # mark the whole class as removed if proportion is one
                if abs(prop - 1.0) <= EPSILON:
                    removed_dc.add(dcind)

                # cut the given proportion from the diameter class
                ret = self.cut_dclass_trees(dcind, (prop, prop), scale=scale)
                sp, cutba, cutn, cutv_com = ret
                if target == 'ba':
                    cuttargetvar = cutba
                else:
                    cuttargetvar = cutn

                self.bottom_targetvar += cuttargetvar
                self.sp_targetvar_taken[sp] += cuttargetvar

            if stop:
                # remove even proportion of trees from each size class to
                # reach the target level
                target_value = self.hdata.cdata['target']
                if target == 'ba':
                    min_N = self.hdata.cdata['min_N_after_thinning']
                    BA_c = self.BA_current
                    if BA_c - target_value > EPSILON and \
                        self.N_current - min_N > EPSILON:
                        prop = (BA_c - target_value) / BA_c
                    else:
                        prop = 0.0
                elif target == 'n':
                    N_c = self.N_current
                    if N_c - target_value > EPSILON:
                        prop = (N_c - target_value) / N_c
                    else:
                        prop = 0.0
                if prop > EPSILON:
                    for dcind in removalorder:
                        self._get_even_proportion_to_remove(dcind,prop,target)

    def cut_from_top_and_bottom(self, scale=None):
        """
        Iterate through the sorted (ascending) treelist and remove tree stems
        until target or minimum BA or minimum N is reached
        """
        removalorder = self.hdata.cdata['removal_order']
        targetvar = self.hdata.cdata['target_var']
        if targetvar!='v_com':
            target = self.hdata.cdata['target']
        ind = self.hdata.ind
        self.bottom_targetvar = 0.0
        self.top_targetvar = 0.0
        removed_dc = set()
        self.sp_targetvar_taken = defaultdict(float)
        stop = False
        onlybottom = False
        for ind in xrange(1, len(removalorder) + 1):
            if stop:
                break
            for cutind, whichend in ((ind - 1, 'bottom'), (-1 * ind, 'top')):
                if stop:
                    break
                dcind = removalorder[cutind]
                if dcind in removed_dc:
                    continue
                propdata = self._get_proportion_to_remove(dcind, whichend,
                                                          scale=scale)
                rprop, stop, skipclass = propdata
                if not skipclass:
                    sp, cutba, cutn, cutv_com = self.cut_dclass_trees(dcind, \
                                                   (rprop, rprop), scale=scale)
                    removed_dc.add(dcind)
                    if targetvar == 'ba':
                        cuttargetvar = cutba
                    elif targetvar == 'n':
                        cuttargetvar = cutn
                    elif targetvar == 'v_com':
                        cuttargetvar = cutv_com
                    if whichend == 'bottom':
                        self.bottom_targetvar += cuttargetvar
                    else:
                        self.top_targetvar += cuttargetvar
                    self.sp_targetvar_taken[sp] += cuttargetvar
            if cutind == -len(removalorder):
                stop = True
            if stop:
                # Target not reached. This is caused by the restrictions
                # (species etc). To reach the target trees should be removed
                # from every diameter class.
                if targetvar=='v_com':
                    # Assumed there is enough removable timber in every
                    # diameter class
                    v_harv = self.hdata.cdata['v_harvested']
                    prop = (v_harv - self.cut_timber) / v_harv
                elif targetvar=='ba':
                    min_N = self.hdata.cdata['min_N_after_thinning']
                    if (self.BA_current - target) > EPSILON and \
                       (self.N_current - min_N) > EPSILON:
                        prop = (self.BA_current - target)/self.BA_current
                    else:
                        prop = 0.0
                elif targetvar == 'n':
                    if (self.N_current - target) > EPSILON:
                        prop = (self.N_current - target)/self.N_current
                    else:
                        prop = 0.0
                if prop > EPSILON:
                    for dcind in removalorder:
                        self._get_even_proportion_to_remove(dcind, prop,
                                                            targetvar)

    def get_dclass_info(self, dclassind):
        """
        Get the tree object, stratum index and species
        """
        trees = self.hdata.treedata
        strata = self.hdata.stratumdata
        ind = self.hdata.ind
        tree = trees[dclassind]
        parentid = tree[ind['tree_parentid']]
        istratum = int(numpy.where(strata[:,ind['stratum_id']] == parentid)[0])
        species = strata[istratum, ind['stratum_sp']]
        return tree, istratum, species

    def reduce_prop_for_reserve_trees(self, dcind, prop):
        """
        if the reserve tree stock is not full, put aside trees if
        they are big enough, those trees are not counted into the harvest
        results, but they are removed from the data
        """
        N_target = self.target_N_reserve
        VOL_target = self.target_VOL_reserve
        if N_target == 0 or VOL_target == 0:
            return prop, 1.0
        N_reserved = self.N_reserve_trees
        VOL_reserved = self.VOL_reserve_trees
        stopnext = self.stopnext
        ind = self.hdata.ind
        tree, istratum, sp = self.get_dclass_info(dcind)
        dbh = tree[ind['tree_d']]
        n = tree[ind['tree_n']]
        stemvolume = self.hdata.result['dclass_vol'][dcind]
        dlimit = self.hdata.cdata['reserve_tree_d_limit']
        potVOL = prop*stemvolume*n
        potN = prop*n
        VOLcond = (VOL_reserved + potVOL) - VOL_target < -EPSILON
        Ncond = (N_reserved + potN) - N_target < -EPSILON
        if dbh - dlimit >= -EPSILON and not stopnext:
            if not VOLcond and not Ncond:
                # Reserve part of the current diameter class
                needed_N = N_target - N_reserved
                reserveprop_N = prop - (needed_N / n)
                needed_VOL = VOL_target - VOL_reserved
                reserveprop_VOL = prop - (needed_VOL / (stemvolume * n))
                if reserveprop_N - reserveprop_VOL >= -EPSILON:
                    reserveprop = reserveprop_N
                else:
                    reserveprop = reserveprop_VOL
                self.N_reserve_trees += (needed_N / (stemvolume * n)) * n
                self.VOL_reserve_trees += (needed_VOL / n) * stemvolume * n
                self.stopnext = True
            elif not Ncond:
                # Reserve part of the current diameter class
                needed = N_target - N_reserved
                reserveprop = prop - (needed / n)
                self.N_reserve_trees = N_target
                self.VOL_reserve_trees += (needed / n) * stemvolume * n
                self.stopnext = True
            elif not VOLcond:
                # Reserve part of the current diameter class
                needed = VOL_target - VOL_reserved
                reserveprop = prop - (needed / (stemvolume * n))
                self.N_reserve_trees += (needed / (stemvolume * n)) * n
                self.VOL_reserve_trees = VOL_target
                self.stopnext = True

            else:
                # Reserve the whole diameter class
                self.N_reserve_trees += potN
                self.VOL_reserve_trees += potVOL
                reserveprop = 0.0
        else:
            reserveprop = prop

        return reserveprop, 1.0

    def reduce_prop_for_seed_and_reserve_trees(self, dcind, prop):
        """
        leave a given number of dominant trees of the given species
        as well as possible reserve trees of any species
        """
        ind = self.hdata.ind
        N_target = self.target_N_seed
        N_reserved = self.N_seed_trees
        tree, istratum, sp = self.get_dclass_info(dcind)
        dbh = tree[ind['tree_d']]
        n = tree[ind['tree_n']]
        # keep the seed trees relatively close to dominant diameter
        dlim = 0.9 * self.hdata.cdata['ddom']
        seedtreespecies=self.hdata.cdata['seedtreespecies']
        if self.seedtree_sp not in seedtreespecies:
            self.seedtree_sp = self.hdata.cdata['seedtreesp2']
        # note that trees should be iterated from the biggest to the smallest
        if N_reserved - N_target < -EPSILON and dbh - dlim >= -EPSILON \
           and sp == self.seedtree_sp:
            if (prop * n) - (N_target - N_reserved) <= EPSILON:
                # Reserve the whole diameter class
                self.N_seed_trees += prop * n
                prop = 0.0
            else:
                # Reserve part of the current diameter class
                prop = 1.0-((N_target - N_reserved) / n * prop)
                self.N_seed_trees = N_target
        reserveprop, dummy = self.reduce_prop_for_reserve_trees(dcind, prop)
        return reserveprop, prop

    def store_harvest_results(self, results, cash_flow, biomass=True,
                              tagmods=True, restype='operation_result'):
        """
        store harvest results into operation result structure, update stratum
        volume and stratum treatment tags and calculate cashflow.
        """
        c_flux = {}
        # NB! This assumes that AREA is used as operation result weight
        results['Area'] = 1.0
        if int(self.hdata.cdata['c_life_cycle']) == 1:
            c_flux['Carbon_flux'] = defaultdict(float)
            c_flux['Carbon_flux_energy'] = defaultdict(float)
        hdind = self.hdata.ind
        svol = self.hdata.result['stratum_vol']
        if len(svol)==0:
            return
        stimber = self.hdata.result['stratum_timber']
        sscrap = self.hdata.result['stratum_scrap']
        if 'data' in restype:
            self._store_harvest_litter(results)
        nchildren = self.hdata.cdata['n_children']
        ao = self.hdata.cdata['ao']
        removedvolume = {}  # volume removed from each stratum
        remvolmerch = 0  # merchantable wood volume removed from each stratum
        scale, removedvolume = self._set_result_scaling(svol)
        # get the species list as operation results should be stored by
        # species
        splist = self.hdata.data[STRATUMLABEL][:,hdind['stratum_sp']]
        spdict = dict(enumerate(splist))
        # add volumes by stratum and assortment to resultlist
        for dataind, factor in enumerate(scale):
            for assort in ao:
                key = (spdict[dataind],assort)
                values = stimber[dataind][assort]
                results['Volume'][key] += scale[dataind] * values[0]
                removedvolume[dataind] += scale[dataind] * values[0]
                # add money by stratum and assortment to resultlist
                results['Income'][key] += factor * values[1]
                cash_flow[0] += factor * values[1]
            #update remvolmerch
            remvolmerch += removedvolume[dataind]
            # add scrapwood by stratum to resultlist
            results['Scrapwood'][key[:1]] += factor * sscrap[dataind]
            removedvolume[dataind] += factor * sscrap[dataind]
            if biomass:
                self._store_harvested_biomass(results, spdict, dataind, c_flux)
        if int(self.hdata.cdata['c_life_cycle']) == 1:
            self._predict_c_lifecycle(c_flux, results)
        # update stratum variables and update total volume of removed wood
        get_ba = self.hdata.get_ba_from_treelist
        get_n = self.hdata.get_n_from_treelist
        strata = self.hdata.stratumdata
        for sind, stratum in enumerate(strata):
            if not sind in nchildren:
                continue
            # tag for removal or modification
            if nchildren[sind] == 0:
                strata[sind, hdind['stratum_tag']] = 2
            elif tagmods:
                strata[sind, hdind['stratum_tag']] = 1
                strata[sind, hdind['stratum_vol']] -= removedvolume[sind]
        return remvolmerch

    def _predict_c_lifecycle(self, c_flux, results):
        """
        Predicts the time and amount of carbon released from the extracted
        wood and stores it as a date-array in operation results
        """
        area = self.hdata.cdata['area']
        mwh2c_ref = self.hdata.cdata['mwh_to_kg_c']
        c_in_kg = self.hdata.cdata['c_in_kg']
        pl = ProductLifeCycle(area, results, c_flux, mwh2c_ref, c_in_kg)
        life_cycle_period = self.hdata.cdata['life_cycle_period']
        pl.c_release(life_cycle_period)
        res = []
        res_energy = []
        for yr, flux in c_flux['Carbon_flux'].iteritems():
            res.append((flux, yr))
        for yr, flux in c_flux['Carbon_flux_energy'].iteritems():
            res_energy.append((flux, yr))
        results['Carbon_flux'] = res
        results['Carbon_flux_energy'] = res_energy


    def _store_harvest_litter(self, results):
        """
        Biomass of the trees being harvested that is left on the site
        """
        bmbranchall = 0.0
        for val in self.hdata.result['bm_branch_all'].values():
            bmbranchall += val
        results['operation_data_results'][0] = bmbranchall
        bmstumpall = 0.0
        for val in self.hdata.result['bm_stump_all'].values():
            bmstumpall += val
        results['operation_data_results'][1] = bmstumpall
        bmrootall = 0.0
        for val in self.hdata.result['bm_root_all'].values():
            bmrootall += val
        results['operation_data_results'][2] = bmrootall
        bmscrapall = 0.0
        for val in self.hdata.result['bm_scrap_all'].values():
            bmscrapall += val
        results['operation_data_results'][3] = bmscrapall

    def _store_harvested_biomass(self, results, spdict, dataind, c_flux):
        """
        Biomass of the tree parts that is harvested for energy production
        """
        bmtimber = self.hdata.result['bm_timber']
        bm_ao = self.hdata.cdata['bm_ao']
        stump_biomass_to_energy = self.hdata.cdata['stump_biomass_to_energy']
        branch_biomass_to_energy = self.hdata.cdata['branch_biomass_to_energy']
        for bm_label in bm_ao:
            assort = bm_label.value
            label = bm_label.label
            key = (spdict[dataind],assort)
            values = bmtimber[dataind][assort]
            results['Biomass'][key] += values[0]
            # add energy of biomass by stratum and assortment to resultlist
            if label == 'stump_roots':
                biomass_to_energy = stump_biomass_to_energy
            else:
                biomass_to_energy = branch_biomass_to_energy
            results['Energy'][key] += values[0] * biomass_to_energy
            # add income of biomass by stratum and assortment to resultlist
            results['Income_biomass'][key] += values[1]
            # Carbon release from energy wood harvested
            if int(self.hdata.cdata['c_life_cycle']) == 1:
                mwh2c_ref = self.hdata.cdata['mwh_to_kg_c']
                c_in_kg = self.hdata.cdata['c_in_kg']
                c_rel, c_rel_energy = energy_conversion(values[0], mwh2c_ref,
                                                        c_in_kg)
                c_flux['Carbon_flux'][0] += c_rel
                c_flux['Carbon_flux_energy'][0] += c_rel_energy

    def store_stand_value_results(self, results):
        """
        set the operation result data values by looping over the strata
        in the order they are in the incoming data
        """
        svol = self.hdata.result['stratum_vol']
        if len(svol)==0:
            return
        stimber = self.hdata.result['stratum_timber']
        bmbranch = self.hdata.result['bm_branch']
        bmstump_root = self.hdata.result['bm_stump_root']
        bmfoliage = self.hdata.result['bm_foliage']
        stump_biomass_to_energy = self.hdata.cdata['stump_biomass_to_energy']
        branch_biomass_to_energy = self.hdata.cdata['branch_biomass_to_energy']
        sscrap = self.hdata.result['stratum_scrap']
        ao = self.hdata.cdata['assortment_order']
        pulplabel = ao[-1]
        scale, rvol = self._set_result_scaling(svol)
        ind = 0
        logvol = 0.0
        pulpvol = 0.0
        for dataind, factor in enumerate(scale):
            # Go through each grade in the current stratum
            for assort in ao:
                values = stimber[dataind][assort]
                results[ind] += values[1] * factor
                if assort!=pulplabel:
                    logvol += values[0]
                else:
                    pulpvol += values[0]
            results[ind+1] = logvol / svol[dataind] * 100
            results[ind+2] = logvol * factor
            results[ind+3] = pulpvol * factor
            results[ind+4] = (logvol + pulpvol) * factor # commercial volume
            results[ind+5] = bmbranch[dataind]
            results[ind+6] = bmstump_root[dataind]
            results[ind+7] = bmfoliage[dataind]
            results[ind+8] = bmbranch[dataind] * branch_biomass_to_energy
            results[ind+9] = bmstump_root[dataind] * stump_biomass_to_energy
            results[ind+10] = bmfoliage[dataind] * branch_biomass_to_energy
            results[ind+11] = sscrap[dataind] * factor
            ind += 11

    def store_tree_value_results(self, results):
        """
        set the operation result data values by looping over the strata
        in the order they are in the incoming data
        """
        svol = self.hdata.result['stratum_vol']
        treevol = self.hdata.result['dclass_vol']
        stimber = self.hdata.result['dclass_timber']
        ao = self.hdata.cdata['assortment_order']
        pulplabel = ao[-1]
        scale, rvol = self._set_result_scaling(svol)
        ind = 0
        for dataind, factor in enumerate(scale):
            # Go through each tree in the current stratum
            for treeid in stimber:
                # Go through each grade in the current stratum
                logvol = 0
                pulpvol = 0
                for assort in ao:
                    values = stimber[treeid][assort]
                    results[ind] += values[1] * factor
                    if assort!=pulplabel:
                        logvol += values[0]
                    else:
                        pulpvol += values[0]
                results[ind+1] = logvol / treevol[treeid] * 100
                results[ind+2] = logvol * factor
                results[ind+3] = pulpvol * factor
                # commercial volume
                results[ind+4] = (logvol + pulpvol) * factor
                # scrapvol
                results[ind+5] = (treevol[treeid]-(logvol + pulpvol)) * factor
                ind += 6

    def store_wood_assortment_results(self, results):
        """
        set the operation result data values by looping over the strata
        in the order they are in the incoming data
        """
        svol = self.hdata.result['stratum_vol']
        stimber = self.hdata.result['stratum_timber']
        ao = self.hdata.cdata['assortment_order']
        pulplabel = ao[-1]
        scale, rvol = self._set_result_scaling(svol)
        ind = 0
        logvol = 0.0
        pulpvol = 0.0
        logvol_plot = 0.0
        pulpvol_plot = 0.0
        for dataind, factor in enumerate(scale):
            # Go through each grade in the current stratum
            for assort in ao:
                values = stimber[dataind][assort]
                if assort!=pulplabel:
                    logvol += values[0]
                    logvol_plot += values[1]
                else:
                    pulpvol += values[0]
                    pulpvol_plot += values[1]
            results[ind] = logvol * factor
            results[ind+1] = pulpvol * factor
            results[ind+2] = logvol_plot * factor
            results[ind+3] = pulpvol_plot * factor
            ind += 4

    def _bucking_algorithm(self, sp, divisor, treestem, pricematrix,
                            timber, pulplabel, log_reduction):
        """
        Bucking algorithm for optimizing the value of tree stem
        """
        divisor = int(divisor * 100)
        nsegments = len(treestem)  # Number of segments
        nclasses = numpy.size(pricematrix, 0)  # of length/diameter-classes
        treestem[:, 1] *= 1000
        volumes = numpy.zeros(nsegments)  # cumulative volume at height h
        values = numpy.zeros(nsegments)  # cumulative value at height h
        grades = numpy.zeros(nsegments)  # grade at height h
        lastlogpos = numpy.zeros(nsegments, int) # position of log bottom at h
        ret = bucking_algorithm(nsegments, nclasses, divisor, treestem,
                                pricematrix, volumes, values, grades,
                                lastlogpos)
        # Find the optimal starting (first cut from the top down) point
        imaxvalue = numpy.argmax(values)
        # Check that the maximum value log is graded
        if abs(grades[imaxvalue]) < EPSILON:
            # tree is too small and is considered as scrapwood
            return 1
        # Cut the tree to optimal length logs and divide the values and
        # volumes to corresponding assortments
        while imaxvalue > 0:
            assortment = grades[imaxvalue]
            timber[assortment][0] += volumes[imaxvalue] - \
                             volumes[lastlogpos[imaxvalue]]
            timber[assortment][1] += values[imaxvalue] - \
                             values[lastlogpos[imaxvalue]]
            imaxvalue = lastlogpos[imaxvalue]
            if abs(imaxvalue - lastlogpos[imaxvalue]) < EPSILON:
                break  # imaxvalue = 0
        # Do the tree level log reduction for each of the assortments
        if self.hdata.cdata['log_reduction_type'] in (1,3):
            timber= self._log_reduction(timber, sp, pulplabel, log_reduction)
        # Do the regional log reduction for each of the assortments
        if self.hdata.cdata['log_reduction_type'] in (2,3):
            timber= self._regional_log_reduction(timber, sp, pulplabel)
        return 1

    def _buck_stem(self, sp, dbh, height, nstems, ptkey, timber, parentid,
                   log_reduction):
        """
        Construct tree stem with taper curve and buck the stem
        """
        pm = self.hdata.cdata['price_matrix']
        divisor = self.hdata.cdata['divisor']
        log_reduction_type = self.hdata.cdata['log_reduction_type']
        stratumvol = self.hdata.result['stratum_vol']
        use_cache = bool(int(self.hdata.cdata['use_cache'])) and \
                self.stem_cache.use_cache
        pricematrix, mindiam, minlen, pulplabel = pm[ptkey]
        sp = int(sp)

        buck = True
        # check if the current tree (or a tree similar enough) is already in
        # the bucking cache - if tree is cached, use the cached values,
        # otherwise buck the stem and store stem volume and timber assortment
        # volumes to the cache
        if use_cache:

            # if using cache, tree dbh and height should be rounded before
            # bucking
            dbh = round(dbh, self.decimals)
            height = round(height, self.decimals)
            stemvol, buck = self.stem_cache.check(timber, sp, dbh, height,
                                    pricematrix, log_reduction, self.decimals)
            treestem = None
            if stemvol == 0.0 and not buck:
                # tree is too small, return just zero
                return 0.0

        if buck:
            # if using cache, tree dbh and height should be rounded before
            # bucking
            if use_cache:
                dbh = round(dbh, self.decimals)
                height = round(height, self.decimals)

            # Calculate the minimum number of segments
            minsegm = int(minlen / (100. * divisor) + 0.5)
            dbh_in_meters = dbh / 100.  # Conversion from cm to m

            # Compute tree stem form and volume with a taper curve function
            vals = self._taper_curve(dbh_in_meters, height, sp, divisor)
            stemvol, treestem = vals
            if treestem is None:
                # Something wrong when calculating tree stem taper and volume
                return 0.0

        # update stratum volume dictionary
        if not parentid in stratumvol:
            stratumvol[parentid] = 0
        stratumvol[parentid] += stemvol * nstems

        if buck:
            # Check that minimun length log top diameter is more than min diam
            if treestem is not None:
                if len(treestem) < minsegm or \
                    (treestem[minsegm-1][1] - (mindiam / 1000)) < -EPSILON:

                    self.stem_cache.update(stemvol, timber, sp, dbh,
                                 height, pricematrix, log_reduction,
                                 self.decimals, True)
                    return stemvol  # buck only big enough trees

            # run the actual bucking algorithm
            if use_cache:
                # if using cache, use a temporary timber volume & value
                # container for the current tree
                tree_timber = dict((a, [0.0, 0.0]) for a in timber.keys())
                ret = self._bucking_algorithm(sp, divisor, treestem,
                                              pricematrix, tree_timber,
                                              pulplabel, log_reduction)

                # update bucking cache with the bucked tree information
                self.stem_cache.update(stemvol, timber, sp, dbh,
                             height, pricematrix, log_reduction, self.decimals)

                # store the values from the temporary timber dictionary to the
                # actual timber dictionary
                for key, val in tree_timber.items():
                    timber[key][0] += val[0]
                    timber[key][1] += val[1]

            else:
                ret = self._bucking_algorithm(sp, divisor, treestem,
                                              pricematrix, timber, pulplabel,
                                              log_reduction)
        return stemvol

    def _calc_biomass(self, istrat, tree, nstems, scrapvol):
        """
        Calculate the biomass in tons of stumps and branches from the removed
        stems of the current diameter class.
        """
        res = self.hdata.result
        cd = self.hdata.cdata
        ind = self.hdata.ind
        stumpdiam = tree[ind['tree_d_s']]
        res['bm_stump_all'][istrat] += tree[ind['tree_bm_stump']] \
                                        * nstems / 1000
        res['bm_root_all'][istrat] += tree[ind['tree_bm_roots']] \
                                        * nstems / 1000
        res['bm_scrap_all'][istrat] += scrapvol * tree[ind['tree_density']] \
                                        * nstems / 1000
        res['bm_branch_all'][istrat] += tree[ind['tree_bm_branches']] \
                                        * nstems / 1000
        if 'bm_timber' in res:
            bm_pm = self.hdata.cdata['bm_price_matrix']
            species = self.hdata.stratumdata[istrat, ind['stratum_sp']]
            bm_ptkey = species
            bm_pricem = bm_pm[bm_ptkey]
        if cd['stump_incl_prop'] > EPSILON and \
           stumpdiam - cd['min_stump_diam'] >= -EPSILON:
            value_stump = tree[ind['tree_bm_stump']] * nstems \
                            * cd['stump_incl_prop'] / 1000
            value_root = tree[ind['tree_bm_roots']] * nstems \
                            * cd['stump_incl_prop'] / 1000
            value_stump_root = value_stump + value_root
            res['bm_root'][istrat] += value_root
            res['bm_stump'][istrat] += value_stump
            res['bm_stump_root'][istrat] += value_stump_root
            if 'bm_timber' in res:
                res['bm_timber'][istrat][1][0] += value_stump_root
                res['bm_timber'][istrat][1][1] += value_stump_root \
                                                    * bm_pricem[0][0][1]
        if cd['biomass_incl_prop'] > EPSILON:
            value_branch = tree[ind['tree_bm_branches']] * nstems \
                            * cd['biomass_incl_prop'] / 1000
            res['bm_branch'][istrat] += value_branch
            if 'bm_timber' in res:
                res['bm_timber'][istrat][2][0] += value_branch
                res['bm_timber'][istrat][2][1] += value_branch \
                                                    * bm_pricem[0][1][1]
            # By default, all of the stemwood which is not assigned as timber
            # or pulp is considered to be energywood.
            value_scrap = scrapvol * tree[ind['tree_density']] * nstems \
                            * cd['biomass_incl_prop'] / 1000
            res['bm_branch'][istrat] += value_scrap
            res['bm_scrap'][istrat] += value_scrap
            if 'bm_timber' in res:
                res['bm_timber'][istrat][2][0] += value_scrap
                res['bm_timber'][istrat][2][1] += value_scrap \
                                                    * bm_pricem[0][1][1]
            # Biomass of foliage
            value_foliage = tree[ind['tree_bm_foliage']] * nstems \
                            * cd['biomass_incl_prop'] / 1000
            res['bm_foliage'][istrat] += value_foliage
            if 'bm_timber' in res:
                res['bm_timber'][istrat][3][0] += value_foliage
                res['bm_timber'][istrat][3][1] += value_foliage \
                                                    * bm_pricem[0][2][1]

    def _check_distr_restr(self, whichend, prop, pot_removal, stopnext):
        cdata = self.hdata.cdata
        skipclass = False
        if whichend == 'bottom':
            limitstr = 'remove_from_below_percentage'
            alreadyremoved = self.bottom_targetvar
        else:
            limitstr = 'remove_from_above_percentage'
            alreadyremoved = self.top_targetvar
        stopp = alreadyremoved / cdata['to_remove'] + EPSILON
        if stopp - cdata[limitstr] >= -EPSILON or abs(cdata[limitstr])<EPSILON:
            skipclass = True
        else:
            potrem = alreadyremoved + pot_removal * prop
            potremprop = potrem / cdata['to_remove']
            if potremprop>cdata[limitstr]:
                remaining = cdata['to_remove'] \
                            * cdata[limitstr] - alreadyremoved
                prop = remaining / pot_removal
                if whichend == 'bottom':
                    togo = self.top_targetvar + alreadyremoved + \
                           pot_removal * prop + EPSILON
                    if togo-cdata['to_remove'] < -EPSILON:
                        stopnext = False
                else:
                    togo = self.bottom_targetvar + alreadyremoved + \
                           pot_removal * prop + EPSILON
                    if togo - cdata['to_remove'] < -EPSILON:
                        stopnext = False
        return prop, skipclass, stopnext

    def _check_sp_restr(self, prop, sp, istratum, pot_removal, targetvar):
        """
        species specific BA or N removal constraints
        NB: for commercial volume target variable sp restriction can't be
        enforced at the moment
        """
        cdata = self.hdata.cdata
        spprop = cdata['sp_proportion'][sp]
        skipclass = False
        ind = self.hdata.ind
        strata = self.hdata.stratumdata
        if targetvar=='ba':
            spba = strata[istratum, ind['stratum_ba']]
            sptargetba = spprop * cdata['target']
            if spba>sptargetba:
                maxbaremoval = prop * pot_removal
                if (spba - maxbaremoval) - sptargetba < -EPSILON:
                    prop = (spba - sptargetba) / maxbaremoval
            else:
                skipclass = True
        elif targetvar == 'n':
            spn = strata[istratum, ind['stratum_N']]
            sptargetn = spprop * cdata['target']
            if spn > sptargetn:
                maxnremoval = prop * pot_removal
                if (spn - maxnremoval) < sptargetn:
                    prop = (spn - sptargetn) / maxnremoval
            else:
                skipclass = True
        return prop, skipclass

    def _get_dclass_ba(self, tree, ind):
        dbh, nstems = tree[ind['tree_d']], tree[ind['tree_n']]
        return math.pi * (dbh / 100 / 2) ** 2 * nstems

    def _get_proportion_to_remove(self, dcind, whichend, scale=None):
        """
        Should the particular diameter class be cut, and if, how much of it?
        """
        ind = self.hdata.ind
        cdata = self.hdata.cdata
        tree, istratum, sp = self.get_dclass_info(dcind)
        n = tree[ind['tree_n']]
        targetvar = cdata['target_var']
        # first the stand level constraints
        if targetvar == 'v_com':
            stratum_scale = scale[istratum]
            res = self.hdata.result
            timbervolume = 0
            for assort, values in res['dclass_timber'][dcind].iteritems():
                timbervolume += values[0] * n * stratum_scale
            # Potentially removable timber volume
            pot_removal = timbervolume
            Vcond = timbervolume - (cdata['v_harvested'] - self.cut_timber) \
                    < -EPSILON
            if Vcond:
                prop = 1.0
                stopnext = False
            else:
                prop = (cdata['v_harvested'] - self.cut_timber) / timbervolume
                stopnext = True
        else:
            if targetvar == 'ba':
                pot_removal = self._get_dclass_ba(tree, ind)
                BAcond = (self.BA_current - pot_removal) - cdata['target'] \
                            > EPSILON
                Ncond = (self.N_current - n) - cdata['min_N_after_thinning'] \
                            > EPSILON
            elif targetvar == 'n':
                pot_removal = n
                BAcond = True
                Ncond = (self.N_current - n) - cdata['target'] > EPSILON
            if BAcond and Ncond:
                prop = 1.0
                stopnext = False
            elif not BAcond and not Ncond:
                propBA = (self.BA_current - cdata['target']) / pot_removal
                propN = (self.N_current - cdata['min_N_after_thinning']) / n
                if propN <= propBA:
                    prop = propN
                else:
                    prop = propBA
                stopnext = False
            elif not BAcond:
                prop = (self.BA_current - cdata['target']) / pot_removal
                stopnext = True
            elif not Ncond:
                if targetvar == 'ba':
                    prop = (self.N_current - cdata['min_N_after_thinning']) / n
                elif targetvar == 'n':
                    prop = (self.N_current - cdata['target']) / n
                stopnext = True
        # distribution and species specific constraints
        skipclass = False
        if tree[ind['tree_d']] - cdata['d_cut_all_below'] >= -EPSILON:
            prop, skipclass, stopnext = self._check_distr_restr(whichend, \
                                                   prop, pot_removal, stopnext)
            if not skipclass and sp in cdata['sp_proportion']:
                prop, skipclass = self._check_sp_restr(prop, sp, istratum,
                                                        pot_removal, targetvar)
        if skipclass:
            prop = 0.0
            stopnext = False
        return prop, stopnext, skipclass

    def _height_of_stump(self, species, d, h):
        """
        calculate the heigth of stump (the lowest cutting point) of tree stem.
        """
        d = d*100 # conversion from m to cm
        if species == 1:
            hs = (0.4456 * d + 0.0952 * h) / 100  # pine
        elif species == 2:
            hs = (0.5089 * d + 0.56 * h) / 100  # spruce
        elif species >= 3:
            hs = (0.4862 * d + 0.4979 * h) / 100  # birch & alder
        if hs < 0.1:
            hs = 0.1
        return hs

    def _log_reduction(self, timber, species, pulplabel, log_reduction):
        """
        reduce the bucked volume and value of logs with a given reduction
        """
        volumered = 0.0  # total volume removed from log assortments
        redp = log_reduction
        totprice = 0.0
        totvol = 0.0
        for assortment, values in timber.iteritems():
            if assortment == pulplabel:
                continue
            volumered += values[0] * (1 -redp)
            values[0] *= redp  # reduce the volume by percentage
            values[1] *= redp  # reduce the value by percentage
            totprice += values[1]
            totvol += values[0]
        if volumered > EPSILON and timber[pulplabel][0] > 0:
            increase = 1 + volumered / timber[pulplabel][0]
            timber[pulplabel][0] *= increase
            timber[pulplabel][1] *= increase
        elif volumered > EPSILON:
            timber[pulplabel][0] = volumered
            timber[pulplabel][1] = totprice/totvol
        return timber

    def _regional_log_reduction(self, timber, species, pulplabel):
        """
        reduce the bucked volume and value of logs with a given reduction
        """
        volumered = 0.0 # total volume removed from log or pulp assortments
        volumeadd = 0.0 # total volume to be add to log assortments
        log_assortment_count = 0 # counts number of log assortments
        redp = self.hdata.cdata['regional_log_reduction']
        totprice = 0.0
        totvol = 0.0
        if redp <= 1:
            for assortment, values in timber.iteritems():
                if assortment == pulplabel:
                    continue
                volumered += values[0] * (1 - redp)
                values[0] *= redp  # reduce the volume by percentage
                values[1] *= redp  # reduce the value by percentage
                totprice += values[1]
                totvol += values[0]
            #increase pulp volume
            if volumered > EPSILON and timber[pulplabel][0] > 0:
                increase = 1 + volumered / timber[pulplabel][0]
                timber[pulplabel][0] *= increase
                timber[pulplabel][1] *= increase
            elif volumered > EPSILON:
                timber[pulplabel][0] = volumered
                timber[pulplabel][1] = totprice/totvol
        else: # redp > 1, ie. add log volume and reduce pulp volume
            for assortment, values in timber.iteritems():
                if assortment == pulplabel:
                    continue
                log_assortment_count += 1
                if timber[pulplabel][0] > EPSILON and \
                   timber[assortment][0] > EPSILON:
                    add_volume = (timber[assortment][0] * redp) \
                                 - timber[assortment][0]
                    if timber[pulplabel][0] - add_volume < -EPSILON:
                        add_volume = timber[pulplabel][0]
                        redp = (add_volume + timber[assortment][0]) \
                                / timber[assortment][0]
                    volumeadd += add_volume
            for assortment, values in timber.iteritems():
                if assortment == pulplabel:
                    if timber[pulplabel][0] > EPSILON:
                        redp_pulp = 1 - volumeadd / timber[pulplabel][0]
                        volumered += values[0] * (1 - redp_pulp)
                        values[0] *= redp_pulp  # reduce the volume by %
                        values[1] *= redp_pulp  # reduce the value by %
                        totprice += values[1]
                        totvol += values[0]
            #increase log volume
            if volumered > EPSILON:
                #Divide removed pulp volume for log assotments
                volumered = volumered / log_assortment_count
                for assortment, values in timber.iteritems():
                    if assortment == pulplabel or abs(timber[assortment][0]) \
                                                                    <= EPSILON:
                        continue
                    increase = 1 + volumered / timber[assortment][0]
                    timber[assortment][0] *= increase
                    timber[assortment][1] *= increase

        return timber

    def _set_result_scaling(self, stratumvol):
        """
        scaling factor for each stratum results based on the stratum volume
        """
        strata = self.hdata.stratumdata
        scale = [None for s in strata]
        removedvolume = {}
        volind = self.hdata.ind['stratum_vol']
        for dataind, stratum in enumerate(strata):
            removedvolume[dataind] = 0
            scale[dataind] = 1  # calculate the volume scaling factor
            if dataind in stratumvol:
                if stratumvol[dataind] >= -EPSILON:
                    scale[dataind] = stratum[volind] / stratumvol[dataind]
        return scale, removedvolume

    def assign_scale(self, stratumvol):
        """
        scaling factor for each stratum results based on the stratum volume
        """
        strata = self.hdata.stratumdata
        scale = [None for s in strata]
        volind = self.hdata.ind['stratum_vol']
        for dataind, stratum in enumerate(strata):
            scale[dataind] = 1  # calculate the volume scaling factor
            if dataind in stratumvol:
                if stratumvol[dataind] >= -EPSILON:
                    scale[dataind] = stratum[volind] / stratumvol[dataind]
        return scale

    def _taper_curve(self, dbh, height, species, divisor):
        """
        Construct tree stem with taper curve method. The actual taper curve
        algorithm is implemented in C (compute_treestem) for performance
        """
        hstump = self._height_of_stump(species, dbh, height)
        piececount = int(math.floor((height - hstump) / divisor) + 1)
        if piececount < 1:
            return 0, None
        treestem = numpy.zeros((piececount, 3))
        stemvol = compute_treestem(treestem, species, piececount, dbh,
                                   height, divisor)
        #For trees shorter than 2 m taper curve returns illogical values
        if height < 2:
            stemvol = 0.001
        return stemvol, treestem

    def _get_even_proportion_to_remove(self, dcind, prop, targetvar):
        """
        How much of every diameter class should be cut, so that target
        basal area or stem number will be reached?
        """
        ind = self.hdata.ind
        cdata = self.hdata.cdata
        tree, istratum, sp = self.get_dclass_info(dcind)
        n = tree[ind['tree_n']]
        # first the stand level constraints
        if targetvar != 'v_com':
            if targetvar == 'ba':
                pot_removal = self._get_dclass_ba(tree, ind)
                BAcond = (self.BA_current - pot_removal) - cdata['target'] \
                            > EPSILON
                Ncond = (self.N_current - n) - cdata['min_N_after_thinning'] \
                            > EPSILON
            elif targetvar == 'n':
                pot_removal = n
                BAcond = True
                Ncond = (self.N_current - n) - cdata['target'] > EPSILON
            if pot_removal>EPSILON:
                if not BAcond and not Ncond:
                    propBA = (self.BA_current - cdata['target']) / pot_removal
                    propN = (self.N_current-cdata['min_N_after_thinning']) / n
                    if propN <= propBA:
                        prop = propN
                    else:
                        prop = propBA
                    stopnext = False
                elif not BAcond:
                    prop = (self.BA_current - cdata['target']) / pot_removal
                    stopnext = True
                elif not Ncond:
                    if targetvar == 'ba':
                        prop = (self.N_current-cdata['min_N_after_thinning'])/n
                    elif targetvar == 'n':
                        prop = (self.N_current - cdata['target']) / n
                    stopnext = True
        # cut the same proportion from every diameter class
        res = self.cut_dclass_trees(dcind, (prop, prop))

