Subversion Repositories public iLand

/********************************************************************************************
**    iLand - an individual based forest landscape and disturbance model
**    http://iland.boku.ac.at
**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl
**
**    This program is free software: you can redistribute it and/or modify
**    it under the terms of the GNU General Public License as published by
**    the Free Software Foundation, either version 3 of the License, or
**    (at your option) any later version.
**
**    This program is distributed in the hope that it will be useful,
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**    GNU General Public License for more details.
**
**    You should have received a copy of the GNU General Public License
**    along with this program.  If not, see <http://www.gnu.org/licenses/>.
********************************************************************************************/
#include "abe_global.h"
#include "actthinning.h"

#include "fmstand.h"
#include "fmtreelist.h"
#include "fmstp.h"
#include "forestmanagementengine.h"
#include "fomescript.h"
#include "model.h"
#include "speciesset.h"
#include "species.h"

#include "tree.h"

#include <QJSValueIterator>
namespace ABE {

/** @class ActThinning
    @ingroup abe
    The ActThinning class implements a very general interface to thinning activties.

  */

// statics
QStringList ActThinning::mSyntaxCustom;


ActThinning::ActThinning(FMSTP *parent): Activity(parent)
{
    mBaseActivity.setIsScheduled(true); // use the scheduler
    mBaseActivity.setDoSimulate(true); // simulate per default
    mThinningType = Invalid;
    if (mSyntaxCustom.isEmpty())
        mSyntaxCustom = QStringList()  << Activity::mAllowedProperties
                                       << "percentile" << "removal" << "thinning"
                                       << "relative" << "remainingStems" << "minDbh"
                                       << "filter" << "targetVariable" << "targetRelative"
                                       << "targetValue" << "classes" << "onEvaluate";
}

QString ActThinning::type() const
{
    QString th;
    switch (mThinningType) {
    case Invalid: th="Invalid"; break;
    case FromBelow: th="from below"; break;
    case FromAbove: th="from above"; break;
    case Custom: th = "custom"; break;
    case Selection: th = "selection"; break;
    }

    return QString("thinning (%1)").arg(th);
}

void ActThinning::setup(QJSValue value)
{
    Activity::setup(value); // setup base events
    mThinningType = Invalid;
    QString th_type = FMSTP::valueFromJs(value, "thinning").toString();
    if (th_type=="fromBelow") mThinningType=FromBelow;
    else if (th_type=="fromAbove") mThinningType=FromAbove;
    else if (th_type=="custom") mThinningType=Custom;
    else if (th_type=="selection") mThinningType=Selection;
    else
        throw IException(QString("Setup of thinning: invalid thinning type: %1").arg(th_type));

    switch (mThinningType) {
    case Custom: setupCustom(value); break;
    case Selection: setupSelective(value); break;
    default: throw IException("No setup defined for thinning type");
    }

}

bool ActThinning::evaluate(FMStand *stand)
{
    bool return_value = true;
    switch (mThinningType) {
    case Custom:
        for (int i=0;i<mCustomThinnings.count();++i)
            return_value = return_value && evaluateCustom(stand, mCustomThinnings[i]);
        return return_value; // false if one fails

    case Selection:
        return evaluateSelective(stand);
    default: throw IException("ActThinning::evaluate: not available for thinning type");
    }
    return false;
}

bool ActThinning::execute(FMStand *stand)
{
    if (stand->trace()) qCDebug(abe) << stand->context() << "execute  activity" << name() << ":" << type();
    if (events().hasEvent(QStringLiteral("onExecute"))) {
        // switch off simulation mode
        stand->currentFlags().setDoSimulate(false);
        // execute this event
        bool result =  Activity::execute(stand);
        stand->currentFlags().setDoSimulate(true);
        return result;
    } else {
        // default behavior: process all marked trees (harvest / cut)
        if (stand->trace()) qCDebug(abe) << stand->context() << "activity" << name() << "remove all marked trees.";
        FMTreeList trees(stand);
        trees.removeMarkedTrees();
        return true;
    }

}

void ActThinning::setupCustom(QJSValue value)
{
    events().setup(value, QStringList() << "onEvaluate");
    mCustomThinnings.clear();
    if (value.hasProperty("thinnings") && value.property("thinnings").isArray()) {
        QJSValueIterator it(value.property("thinnings"));
        while (it.hasNext()) {
            it.next();
            if (it.name()==QStringLiteral("length"))
                continue;
            mCustomThinnings.push_back(SCustomThinning());
            setupSingleCustom(it.value(), mCustomThinnings.back());
        }
    } else {
        mCustomThinnings.push_back(SCustomThinning());
        setupSingleCustom(value, mCustomThinnings.back());
    }
}

void ActThinning::setupSelective(QJSValue value)
{
    mSelectiveThinning.N = FMSTP::valueFromJs(value, "N", "400").toInt();
}

// setup of the "custom" thinning operation
void ActThinning::setupSingleCustom(QJSValue value, SCustomThinning &custom)
{
    FMSTP::checkObjectProperties(value, mSyntaxCustom, "setup of 'custom' thinning:" + name());

    custom.usePercentiles = FMSTP::boolValueFromJs(value, "percentile", true);
    custom.removal = FMSTP::boolValueFromJs(value, "removal", true);
    custom.relative = FMSTP::boolValueFromJs(value, "relative", true);
    custom.remainingStems = FMSTP::valueFromJs(value, "remainingStems", "0").toInt();
    custom.minDbh = FMSTP::valueFromJs(value, "minDbh", "0").toNumber();
    QJSValue filter = FMSTP::valueFromJs(value, "filter", "");
    if (filter.isString())
        custom.filter = filter.toString();
    else
        custom.filter = QString();
    custom.targetVariable = FMSTP::valueFromJs(value, "targetVariable", "stems").toString();
    if (custom.targetVariable != "stems" &&
            custom.targetVariable != "basalArea" &&
            custom.targetVariable != "volume")
        throw IException(QString("setup of custom Activity: invalid targetVariable: %1").arg(custom.targetVariable));

    custom.targetRelative = FMSTP::boolValueFromJs(value, "targetRelative", true);
    custom.targetValue = FMSTP::valueFromJs(value, "targetValue", "30").toNumber();
    if (custom.targetRelative && (custom.targetValue>100. || custom.targetValue<0.))
        throw IException(QString("setup of custom Activity: invalid relative targetValue (0-100): %1").arg(custom.targetValue));

    QJSValue values = FMSTP::valueFromJs(value, "classes",  "", "setup custom acitvity");
    if (!values.isArray())
        throw IException("setup of custom activity: the 'classes' is not an array.");
    custom.classValues.clear();
    custom.classPercentiles.clear();
    QJSValueIterator it(values);
    while (it.hasNext()) {
        it.next();
        if (it.name()==QStringLiteral("length"))
            continue;
        custom.classValues.push_back(it.value().toNumber());
    }
    if (custom.classValues.size()==0)
        throw IException("setup of custom thinnings: 'classes' has no elements.");

    // check if sum is 100 for relative classes
    if (custom.relative) {
        double sum=0.;
        for (int i=0;i<custom.classValues.size();++i)
            sum+=custom.classValues[i];
        if (fabs(sum-100.)>0.000001)
            throw IException("setup of custom thinnings: 'classes' do not add up to 100 (relative=true).");
    }

    // span the range between 0..100: from e.g. 10,20,30,20,20 -> 0,10,30,60,80,100
    double f = 100. / custom.classValues.size();
    double p = 0.;
    for (int i=0;i<custom.classValues.size();++i, p+=f)
        custom.classPercentiles.push_back(qRound(p));
    custom.classPercentiles.push_back(100);
}

bool ActThinning::evaluateCustom(FMStand *stand, SCustomThinning &custom)
{
    // fire onEvaluate event and collect probabilities
    QJSValue eval_result = events().run(QStringLiteral("onEvaluate"), stand);
    if (eval_result.isBool() && eval_result.toBool()==false)
        return false; // do nothing
    bool species_selective = false;

    if( eval_result.isObject()) {
        // expecting a list of probabilities....
        // create list if not present
        if (mSpeciesSelectivity.isEmpty()) {
            foreach(const Species *s, GlobalSettings::instance()->model()->speciesSet()->activeSpecies())
                mSpeciesSelectivity[s] = 1.;
        }
        // fetch from javascript
        double rest_val = eval_result.property(QStringLiteral("rest")).isNumber() ?  eval_result.property(QStringLiteral("rest")).toNumber() : 1.;
        foreach(const Species *s, mSpeciesSelectivity.keys()) {
            mSpeciesSelectivity[s] = limit( eval_result.property(s->id()).isNumber() ? eval_result.property(s->id()).toNumber() : rest_val, 0., 1.);
        }
        species_selective = true;
    }


    FMTreeList trees(stand);
    QString filter = custom.filter;
    if (custom.minDbh>0.) {
        if (!filter.isEmpty())
            filter += " and ";
        filter += QString("dbh>%1").arg(custom.minDbh);
    }

    if (!filter.isEmpty())
        trees.load(filter);
    else
        trees.loadAll();

    if (custom.remainingStems>0 && custom.remainingStems>=trees.trees().size())
        return false;

    if (trees.trees().size()==0)
        return false;

    // remove harvest flags.
    clearTreeMarks(&trees);

    // sort always by target variable (if it is stems, then simply by dbh)
    bool target_dbh = custom.targetVariable=="stems";
    if (target_dbh)
        trees.sort("dbh");
    else
        trees.sort(custom.targetVariable);

    // count trees and values (e.g. volume) in the defined classes
    QVector<double> values = custom.classValues; // make a copy
    QVector<double> tree_counts = custom.classValues; // make a copy
    QVector<int> percentiles = custom.classPercentiles; // make a copy

    QVector<QPair<Tree*, double> >::const_iterator it;
    for (int i=0;i<values.count();++i) {
        tree_counts[i]=0.;
    }
    int class_index = 0;
    int n=0;

    percentiles.first()=0;
    double tree_count=trees.trees().count();
    double total_value = 0.;
    for (it = trees.trees().constBegin(); it!=trees.trees().constEnd(); ++it, ++n) {
        if (n/tree_count*100. > custom.classPercentiles[class_index+1]) {
            ++class_index;
            percentiles[class_index] = n; // then n'th tree
        }
        tree_counts[class_index]++;
        total_value += target_dbh?1.:it->second; // e.g., sum of volume in the class, or simply count
    }
    while (++class_index<percentiles.size())
        percentiles[class_index]=n+1;

    double target_value=0.;
    if (custom.targetRelative)
        target_value = custom.targetValue * total_value / 100.;
    else
        target_value = custom.targetValue * stand->area();

    if (!custom.relative) {
        // TODO: does not work now!!! redo!!
        // class values are given in absolute terms, e.g. 40m3/ha.
        // this needs to be translated to relative classes.
        // if the demand in a class cannot be met (e.g. planned removal of 40m3/ha, but there are only 20m3/ha in the class),
        // then the "miss" is distributed to the other classes (via the scaling below).
        for (int i=0;i<values.size();++i) {
            if (values[i]>0){
                if (values[i]<=custom.classValues[i]*stand->area()) {
                    values[i] = 1.; // take all from the class
                } else {
                    values[i] = custom.classValues[i]*stand->area() / values[i];
                }
            }
        }
        // scale to 100
        double sum=0.;
        for (int i=0;i<values.size();++i)
            sum+=values[i];
        if (sum>0.){
            for (int i=0;i<values.size();++i)
                values[i] *= 100./sum;
        }
    }

    // *****************************************************************
    // ***************    Main loop
    // *****************************************************************
    for (int i=0;i<values.size();++i)
        values[i] = 0;

    bool finished = false;
    int cls;
    double p;
    int removed_trees = 0;
    double removed_value = 0.;
    int no_tree_found = 0;
    bool target_value_reached=false;
    do {
        // look up a random number: it decides in which class to select a tree:
        p = nrandom(0,100);
        for (cls=0;cls<values.size();++cls) {
            if (p < custom.classPercentiles[cls+1])
                break;
        }
        // select a tree:
        int tree_idx = selectRandomTree(&trees, percentiles[cls], percentiles[cls+1]-1, species_selective);
        if (tree_idx>=0) {
            // stop harvesting, when the target size is reached: if the current tree would surpass the limit,
            // a random number decides whether the tree should be included or not.
            double tree_value = target_dbh?1.:trees.trees()[tree_idx].second;
            if (custom.targetValue>0.) {
                if (removed_value + tree_value > target_value) {
                    if (drandom()>0.5 || target_value_reached)
                        break;
                    else
                        target_value_reached = true;
                }

            }
            trees.remove_single_tree(tree_idx, true);
            removed_trees++;
            removed_value += tree_value;
            values[cls]++;

        } else {
            // tree_idx = -1: no tree found in list, -2: tree found but is not selected
            no_tree_found += tree_idx == -1 ? 100 : 1; // empty list counts much more
            if (no_tree_found > 1000)
                finished=true;
        }
        // stop harvesting, when the minimum remaining number of stems is reached
        if (trees.trees().size()-removed_trees <= custom.remainingStems*stand->area())
            finished = true;

        if (custom.targetValue>0. && removed_value > target_value)
            finished = true;

    } while (!finished);

    if (stand->trace()) {
        qCDebug(abe) << stand->context() << "custom-thinning: removed" << removed_trees << ". Reached cumulative 'value' of:" << removed_value << "(planned value:" << target_value << "). #of no trees found:" << no_tree_found << "; stand-area:" << stand->area();
        for (int i=0;i<values.count();++i)
            qCDebug(abe) << stand->context() << "class" << i << ": removed" << values[i] << "of" << percentiles[i+1]-percentiles[i];
    }

    return true;

}

int ActThinning::selectRandomTree(FMTreeList *list, const int pct_min, const int pct_max, const bool selective)
{
    // pct_min, pct_max: the indices of the first and last tree in the list to be looked for, including pct_max
    // seek a tree in the class 'cls' (which has not already been removed);
    int idx;
    if (pct_max<pct_min)
        return -1;
    // search randomly for a couple of times
    for (int i=0;i<5;i++) {
        idx = irandom(pct_min, pct_max);
        Tree *tree = list->trees().at(idx).first;
        if (!tree->isDead() && !tree->isMarkedForHarvest() && !tree->isMarkedForCut())
            return selectSelectiveSpecies(list, selective, idx);
    }
    // not found, now walk in a random direction...
    int direction = 1;
    if (drandom()>0.5) direction=-1;
    // start in one direction from the last selected random position
    int ridx=idx;
    while (ridx>=pct_min && ridx<pct_max) {
        Tree *tree = list->trees().at(ridx).first;
        if (!tree->isDead() && !tree->isMarkedForHarvest() && !tree->isMarkedForCut())
            return selectSelectiveSpecies(list, selective, ridx);

        ridx+=direction;
    }
    // now look in the other direction
    direction = -direction;
    ridx = idx;
    while (ridx>=pct_min && ridx<pct_max) {
        Tree *tree = list->trees().at(ridx).first;
        if (!tree->isDead() && !tree->isMarkedForHarvest() && !tree->isMarkedForCut())
            return selectSelectiveSpecies(list, selective, ridx);

        ridx+=direction;
    }

    // no tree found in the entire range
    return -1;


}

int ActThinning::selectSelectiveSpecies(FMTreeList *list, const bool is_selective, const int index)
{
    if (!is_selective)
        return index;
    // check probability for species [0..1, 0: 0% chance to take a tree of that species] against a random number
    if (mSpeciesSelectivity[list->trees()[index].first->species()] < drandom())
        return index; // take the tree

    // a tree was found but is not going to be removed
    return -2;
}

void ActThinning::clearTreeMarks(FMTreeList *list)
{
    QVector<QPair<Tree*, double> >::const_iterator it;
    for (it=list->trees().constBegin(); it!=list->trees().constEnd(); ++it) {
        Tree *tree = it->first;
        if (tree->isMarkedForHarvest())
            tree->markForHarvest(false);
        if (tree->isMarkedForCut())
            tree->markForCut(false);

    }
}

bool ActThinning::evaluateSelective(FMStand *stand)
{
    markCropTrees(stand);
    return true;
}

bool ActThinning::markCropTrees(FMStand *stand)
{
    // tree list from current exeution context
    FMTreeList *treelist = ForestManagementEngine::instance()->scriptBridge()->treesObj();
    treelist->setStand(stand);
    treelist->loadAll();
    clearTreeMarks(treelist);

    // get the 2x2m grid for the current stand
    Grid<float> &grid = treelist->localGrid();
    // clear (except the out of "stand" pixels)
    for (float *p=grid.begin(); p!=grid.end(); ++p)
        if (*p > -1.f)
            *p = 0.f;

    int target_n = mSelectiveThinning.N * stand->area();

    if (target_n>=treelist->trees().count())
        target_n = treelist->trees().count();

    int max_target_n = qMax(target_n * 1.5, treelist->trees().count()/2.);
    if (max_target_n>=treelist->trees().count())
        max_target_n = treelist->trees().count();
    // we have 2500 px per ha (2m resolution)
    // if each tree dominates its Moore-neighborhood, 2500/9 = 267 trees are possible (/ha)
    // if *more* trees should be marked, some trees need to be on neighbor pixels:
    // pixels = 2500 / N; if 9 px are the Moore neighborhood, the "overprint" is N*9 / 2500.
    // N*)/2500 -1 = probability of having more than zero overlapping pixels
    double overprint = (mSelectiveThinning.N * 9) / double(cPxPerHectare) - 1.;

    // order the list of trees according to tree height
    treelist->sort("-height");

    // start with a part of N and 0 overlap
    int n_found = 0;
    int tests=0;
    for (int i=0;i<target_n/3;i++) {
        float f=testPixel(treelist->trees().at(i).first->position(), grid); ++tests;
        if (f==0.f) {
            setPixel(treelist->trees().at(i).first->position(), grid);
            treelist->trees()[i].first->markCropTree(true);
            ++n_found;
        }
    }
    // continue with a higher probability --- incr
    for (int run=0;run<4;++run) {
        for (int i=0; i<max_target_n;++i) {
            if (treelist->trees().at(i).first->isMarkedAsCropTree())
                continue;

            float f=testPixel(treelist->trees().at(i).first->position(), grid); ++tests;

            if ( (f==0.f) ||
                 (f<=2.f && drandom()<overprint) ||
                 (run==1 && f<=4 && drandom()<overprint) ||
                 (run==2 && drandom()<overprint) ||
                 (run==3) ) {

                setPixel(treelist->trees().at(i).first->position(), grid);
                treelist->trees()[i].first->markCropTree(true);
                ++n_found;
                if (n_found == target_n)
                    break;
            }
        }
        if (n_found==target_n)
            break;
    }

    // now mark the competitors:
    // competitors are trees up to 75th percentile of the tree population that
    int n_competitor=0;
    for (int run=0;run<3;++run) {
        for (int i=0; i<max_target_n;++i) {
            Tree *tree = treelist->trees().at(i).first;
            if (tree->isMarkedAsCropTree() || tree->isMarkedAsCropCompetitor())
                continue;

            float f=testPixel(treelist->trees().at(i).first->position(), grid); ++tests;

            if ( (f>12.f) ||
                 (run==1 && f>8) ||
                 (run==2 && f>4) ) {
                tree->markCropCompetitor(true);
                n_competitor++;
            }
        }
    }


    if (FMSTP::verbose()) {
        qCDebug(abe) << stand->context() << "Thinning::markCropTrees: marked" << n_found << "(plan:" << target_n << ") from total" << treelist->trees().count()
                     << ". Tests performed:" << tests << "marked as competitors:" << n_competitor;
    }
    return n_found==target_n;

}

float ActThinning::testPixel(const QPointF &pos, Grid<float> &grid)
{
    // check Moore neighborhood
    int x=grid.indexAt(pos).x();
    int y=grid.indexAt(pos).y();

    float sum = 0.f;
    sum += grid.isIndexValid(x-1,y-1) ? grid.valueAtIndex(x-1, y-1) : 0;
    sum += grid.isIndexValid(x,y-1) ? grid.valueAtIndex(x, y-1) : 0;
    sum += grid.isIndexValid(x+1,y-1) ? grid.valueAtIndex(x+1, y-1) : 0;

    sum += grid.isIndexValid(x-1,y) ? grid.valueAtIndex(x-1, y) : 0;
    sum += grid.isIndexValid(x,y) ? grid.valueAtIndex(x, y) : 0;
    sum += grid.isIndexValid(x+1,y) ? grid.valueAtIndex(x+1, y) : 0;

    sum += grid.isIndexValid(x-1,y+1) ? grid.valueAtIndex(x-1, y+1) : 0;
    sum += grid.isIndexValid(x,y+1) ? grid.valueAtIndex(x, y+1) : 0;
    sum += grid.isIndexValid(x+1,y+1) ? grid.valueAtIndex(x+1, y+1) : 0;

    return sum;
}

void ActThinning::setPixel(const QPointF &pos, Grid<float> &grid)
{
    // check Moore neighborhood
    int x=grid.indexAt(pos).x();
    int y=grid.indexAt(pos).y();

    if (grid.isIndexValid(x-1,y-1)) grid.valueAtIndex(x-1, y-1)++;
    if (grid.isIndexValid(x,y-1)) grid.valueAtIndex(x, y-1)++;
    if (grid.isIndexValid(x+1,y-1)) grid.valueAtIndex(x+1, y-1)++;

    if (grid.isIndexValid(x-1,y)) grid.valueAtIndex(x-1, y)++;
    if (grid.isIndexValid(x,y)) grid.valueAtIndex(x, y) += 3; // more impact on center pixel
    if (grid.isIndexValid(x+1,y)) grid.valueAtIndex(x+1, y)++;

    if (grid.isIndexValid(x-1,y+1)) grid.valueAtIndex(x-1, y+1)++;
    if (grid.isIndexValid(x,y+1)) grid.valueAtIndex(x, y+1)++;
    if (grid.isIndexValid(x+1,y+1)) grid.valueAtIndex(x+1, y+1)++;
}


} // namespace