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/>.

********************************************************************************************/

/** @class Model

  Main object of the iLand model composited of various sub models / sub components.

  @ingroup core

  The class Model is the top level container of iLand. The Model holds a collection of ResourceUnits, links to SpeciesSet and Climate.

  ResourceUnit are grid cells with (currently) a size of 1 ha (100x100m). Many stand level processes (NPP produciton, WaterCycle) operate on this

  level.

  The Model also contain the landscape-wide 2m LIF-grid (http://iland.boku.ac.at/competition+for+light).

  */

#include "global.h"

#include "model.h"

#include "sqlhelper.h"

#include "xmlhelper.h"

#include "debugtimer.h"

#include "environment.h"

#include "timeevents.h"

#include "helper.h"

#include "resourceunit.h"

#include "climate.h"

#include "speciesset.h"

#include "standloader.h"

#include "tree.h"

#include "management.h"

#include "modelsettings.h"

#include "standstatistics.h"

#include "mapgrid.h"

#include "modelcontroller.h"

#include "modules.h"

#include "dem.h"

#include "outputmanager.h"

#include "forestmanagementengine.h"

#include <QtCore>

#include <QtXml>

/** iterate over all trees of the model. return NULL if all trees processed.

  Usage:

  @code

  AllTreeIterator trees(model);

  while (Tree *tree = trees.next()) { // returns NULL when finished.

     tree->something(); // do something

  }

  @endcode  */

Tree *AllTreeIterator::next()

{

    if (!mTreeEnd) {

        // initialize to first ressource unit

        mRUIterator = mModel->ruList().constBegin();

        // fast forward to the first RU with trees

        while (mRUIterator!=mModel->ruList().constEnd()) {

            if ((*mRUIterator)->trees().count()>0)

                break;

            ++mRUIterator;

        }

            // finished if all RU processed

        if (mRUIterator == mModel->ruList().constEnd())

            return NULL;

        mTreeEnd = &((*mRUIterator)->trees().back()) + 1; // let end point to "1 after end" (STL-style)

        mCurrent = &((*mRUIterator)->trees().front());

    }

    if (mCurrent==mTreeEnd) {

        ++mRUIterator; // switch to next RU (loop until RU with trees is found)

        while (mRUIterator!=mModel->ruList().constEnd()) {

            if ((*mRUIterator)->trees().count()>0) {

                break;

            }

            ++mRUIterator;

        }

        if (mRUIterator == mModel->ruList().constEnd()) {

            mCurrent = NULL;

            return NULL; // finished!!

        }else {

            mTreeEnd = &((*mRUIterator)->trees().back()) + 1;

            mCurrent = &((*mRUIterator)->trees().front());

        }

    }

    return mCurrent++;

}

Tree *AllTreeIterator::nextLiving()

{

    while (Tree *t = next())

        if (!t->isDead()) return t;

    return NULL;

}

Tree *AllTreeIterator::current() const

{

    return mCurrent?mCurrent-1:NULL;

}

ModelSettings Model::mSettings;

Model::Model()

{

    initialize();

    GlobalSettings::instance()->setModel(this);

    GlobalSettings::instance()->resetScriptEngine(); // clear the script

    QString dbg="running in release mode.";

    DBGMODE( dbg="running in debug mode."; );

    qDebug() << dbg;

}

Model::~Model()

{

    clear();

    GlobalSettings::instance()->setModel(NULL);

}

/** Initial setup of the Model.

  */

void Model::initialize()

{

   mSetup = false;

   GlobalSettings::instance()->setCurrentYear(0);

   mGrid = 0;

   mHeightGrid = 0;

   mManagement = 0;

   mABEManagement = 0;

   mEnvironment = 0;

   mTimeEvents = 0;

   mStandGrid = 0;

   mModules = 0;

   mDEM = 0;

}

/** sets up the simulation space.

*/

void Model::setupSpace()

{

    XmlHelper xml(GlobalSettings::instance()->settings().node("model.world"));

    double cellSize = xml.value("cellSize", "2").toDouble();

    double width = xml.value("width", "100").toDouble();

    double height = xml.value("height", "100").toDouble();

    double buffer = xml.value("buffer", "60").toDouble();

    mModelRect = QRectF(0., 0., width, height);

    qDebug() << QString("setup of the world: %1x%2m with cell-size=%3m and %4m buffer").arg(width).arg(height).arg(cellSize).arg(buffer);

    QRectF total_grid(QPointF(-buffer, -buffer), QPointF(width+buffer, height+buffer));

    qDebug() << "setup grid rectangle:" << total_grid;

    if (mGrid)

        delete mGrid;

    mGrid = new FloatGrid(total_grid, cellSize);

    mGrid->initialize(1.f);

    if (mHeightGrid)

        delete mHeightGrid;

    mHeightGrid = new HeightGrid(total_grid, cellSize*5);

    mHeightGrid->wipe(); // set all to zero

    Tree::setGrid(mGrid, mHeightGrid);

    // setup the spatial location of the project area

    if (xml.hasNode("location")) {

        // setup of spatial location

        double loc_x = xml.valueDouble("location.x");

        double loc_y = xml.valueDouble("location.y");

        double loc_z = xml.valueDouble("location.z");

        double loc_rot = xml.valueDouble("location.rotation");

        setupGISTransformation(loc_x, loc_y, loc_z, loc_rot);

        qDebug() << "setup of spatial location: x/y/z" << loc_x << loc_y << loc_z << "rotation:" << loc_rot;

    } else {

        setupGISTransformation(0., 0., 0., 0.);

    }

    // load environment (multiple climates, speciesSets, ...

    if (mEnvironment)

        delete mEnvironment;

    mEnvironment = new Environment();

    if (xml.valueBool("environmentEnabled", false)) {

        QString env_file = GlobalSettings::instance()->path(xml.value("environmentFile"));

        bool grid_mode = (xml.value("environmentMode")=="grid");

        QString grid_file = GlobalSettings::instance()->path(xml.value("environmentGrid"));

        if (grid_mode && QFile::exists(grid_file))

            mEnvironment->setGridMode(grid_file);

        if (!mEnvironment->loadFromFile(env_file))

            return;

        // retrieve species sets and climates:

        mSpeciesSets << mEnvironment->speciesSetList();

        mClimates << mEnvironment->climateList();

    } else {

        // load and prepare default values

        // (2) SpeciesSets: currently only one a global species set.

        SpeciesSet *speciesSet = new SpeciesSet();

        mSpeciesSets.push_back(speciesSet);

        speciesSet->setup();

        // Climate...

        Climate *c = new Climate();

        mClimates.push_back(c);

        mEnvironment->setDefaultValues(c, speciesSet);

    } // environment?

    // time series data

    if (xml.valueBool(".timeEventsEnabled", false)) {

        mTimeEvents = new TimeEvents();

        mTimeEvents->loadFromFile(GlobalSettings::instance()->path(xml.value("timeEventsFile"), "script"));

    }

    // simple case: create ressource units in a regular grid.

    if (xml.valueBool("resourceUnitsAsGrid")) {

        mRUmap.setup(QRectF(0., 0., width, height),100.); // Grid, that holds positions of resource units

        mRUmap.wipe();

        bool mask_is_setup = false;

        if (xml.valueBool("standGrid.enabled")) {

            QString fileName = GlobalSettings::instance()->path(xml.value("standGrid.fileName"));

            mStandGrid = new MapGrid(fileName,false); // create stand grid index later

            if (mStandGrid->isValid()) {

                for (int i=0;i<mStandGrid->grid().count();i++) {

                    const int &grid_value = mStandGrid->grid().constValueAtIndex(i);

                    mHeightGrid->valueAtIndex(i).setValid( grid_value > -1 );

                    if (grid_value>-1)

                        mRUmap.valueAt(mStandGrid->grid().cellCenterPoint(i)) = (ResourceUnit*)1;

                    if (grid_value < -1)

                        mHeightGrid->valueAtIndex(i).setForestOutside(true);

                }

            }

            mask_is_setup = true;

        } else {

            if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {

                // in the case we have no stand grid but only a large rectangle (without the torus option)

                // we assume a forest outside

                for (int i=0;i<mHeightGrid->count();++i) {

                    const QPointF &p = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));

                    if (p.x() < 0. || p.x()>width || p.y()<0. || p.y()>height) {

                        mHeightGrid->valueAtIndex(i).setForestOutside(true);

                        mHeightGrid->valueAtIndex(i).setValid(false);

                    }

                }

            }

        }

        ResourceUnit **p; // ptr to ptr!

        ResourceUnit *new_ru;

        int ru_index = 0;

        for (p=mRUmap.begin(); p!=mRUmap.end(); ++p) {

            QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));

            mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.

            if (!mStandGrid || !mStandGrid->isValid() || *p>0) {

                // create resource units for valid positions only

                new_ru = new ResourceUnit(ru_index++); // create resource unit

                new_ru->setClimate( mEnvironment->climate() );

                new_ru->setSpeciesSet( mEnvironment->speciesSet() );

                new_ru->setup();

                new_ru->setID( mEnvironment->currentID() ); // set id of resource unit in grid mode

                new_ru->setBoundingBox(r);

                mRU.append(new_ru);

                *p = new_ru; // save in the RUmap grid

            }

        }

        if (mStandGrid && mStandGrid->isValid())

            mStandGrid->createIndex();

        // now store the pointers in the grid.

        // Important: This has to be done after the mRU-QList is complete - otherwise pointers would

        // point to invalid memory when QList's memory is reorganized (expanding)

//        ru_index = 0;

//        for (p=mRUmap.begin();p!=mRUmap.end(); ++p) {

//            *p = mRU.value(ru_index++);

//        }

        qDebug() << "created a grid of ResourceUnits: count=" << mRU.count() << "number of RU-map-cells:" << mRUmap.count();

        calculateStockableArea();

        // setup of the project area mask

        if (!mask_is_setup && xml.valueBool("areaMask.enabled", false) && xml.hasNode("areaMask.imageFile")) {

            // to be extended!!! e.g. to load ESRI-style text files....

            // setup a grid with the same size as the height grid...

            FloatGrid tempgrid((int)mHeightGrid->cellsize(), mHeightGrid->sizeX(), mHeightGrid->sizeY());

            QString fileName = GlobalSettings::instance()->path(xml.value("areaMask.imageFile"));

            loadGridFromImage(fileName, tempgrid); // fetch from image

            for (int i=0;i<tempgrid.count(); i++)

                mHeightGrid->valueAtIndex(i).setValid( tempgrid.valueAtIndex(i)>0.99 );

            qDebug() << "loaded project area mask from" << fileName;

        }

        // list of "valid" resource units

        QList<ResourceUnit*> valid_rus;

        foreach(ResourceUnit* ru, mRU)

            if (ru->id()!=-1)

                valid_rus.append(ru);

        // setup of the digital elevation map (if present)

        QString dem_file = xml.value("DEM");

        if (!dem_file.isEmpty()) {

            mDEM = new DEM(GlobalSettings::instance()->path(dem_file));

            // add them to the visuals...

            GlobalSettings::instance()->controller()->addGrid(mDEM, "DEM height", GridViewRainbow, 0, 1000);

            GlobalSettings::instance()->controller()->addGrid(mDEM->slopeGrid(), "DEM slope", GridViewRainbow, 0, 3);

            GlobalSettings::instance()->controller()->addGrid(mDEM->aspectGrid(), "DEM aspect", GridViewRainbow, 0, 360);

            GlobalSettings::instance()->controller()->addGrid(mDEM->viewGrid(), "DEM view", GridViewGray, 0, 1);

        }

        // setup of external modules

        mModules->setup();

        if (mModules->hasSetupResourceUnits()) {

            for (ResourceUnit **p=mRUmap.begin(); p!=mRUmap.end(); ++p) {

                QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));

                mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.

                mModules->setupResourceUnit( *p );

            }

        }

        // setup of scripting environment

        ScriptGlobal::setupGlobalScripting();

        // setup the helper that does the multithreading

        threadRunner.setup(valid_rus);

        threadRunner.setMultithreading(GlobalSettings::instance()->settings().valueBool("system.settings.multithreading"));

        threadRunner.print();

    } else  {

        throw IException("resourceUnitsAsGrid MUST be set to true - at least currently :)");

    }

    mSetup = true;

}

/** clear() frees all ressources allocated with the run of a simulation.

  */

void Model::clear()

{

    mSetup = false;

    qDebug() << "Model clear: attempting to clear" << mRU.count() << "RU, " << mSpeciesSets.count() << "SpeciesSets.";

    // clear ressource units

    qDeleteAll(mRU); // delete ressource units (and trees)

    mRU.clear();

    qDeleteAll(mSpeciesSets); // delete species sets

    mSpeciesSets.clear();

    // delete climate data

    qDeleteAll(mClimates);

    // delete the grids

    if (mGrid)

        delete mGrid;

    if (mHeightGrid)

        delete mHeightGrid;

    if (mManagement)

        delete mManagement;

    if (mEnvironment)

        delete mEnvironment;

    if (mTimeEvents)

        delete mTimeEvents;

    if (mStandGrid)

        delete mStandGrid;

    if (mModules)

        delete mModules;

    if (mDEM)

        delete mDEM;

    if (mABEManagement)

        delete mABEManagement;

    mGrid = 0;

    mHeightGrid = 0;

    mManagement = 0;

    mEnvironment = 0;

    mTimeEvents = 0;

    mStandGrid  = 0;

    mModules = 0;

    mDEM = 0;

    mABEManagement = 0;

    GlobalSettings::instance()->outputManager()->close();

    qDebug() << "Model ressources freed.";

}

/** Setup of the Simulation.

  This really creates the simulation environment and does the setup of various aspects.

  */

void Model::loadProject()

{

    DebugTimer dt("load project");

    GlobalSettings *g = GlobalSettings::instance();

    g->printDirecories();

    const XmlHelper &xml = g->settings();

    g->clearDatabaseConnections();

    // database connections: reset

    GlobalSettings::instance()->clearDatabaseConnections();

    // input and climate connection

    // see initOutputDatabase() for output database

    QString dbPath = g->path( xml.value("system.database.in"), "database");

    GlobalSettings::instance()->setupDatabaseConnection("in", dbPath, true);

    dbPath = g->path( xml.value("system.database.climate"), "database");

    GlobalSettings::instance()->setupDatabaseConnection("climate", dbPath, true);

    mSettings.loadModelSettings();

    mSettings.print();

    // random seed: if stored value is <> 0, use this as the random seed (and produce hence always an equal sequence of random numbers)

    uint seed = xml.value("system.settings.randomSeed","0").toUInt();

    RandomGenerator::setup(RandomGenerator::ergMersenneTwister, seed); // use the MersenneTwister as default

    // linearization of expressions: if true *and* linearize() is explicitely called, then

    // function results will be cached over a defined range of values.

    bool do_linearization = xml.valueBool("system.settings.expressionLinearizationEnabled", false);

    Expression::setLinearizationEnabled(do_linearization);

    if (do_linearization)

        qDebug() << "The linearization of certains expressions is enabled (performance optimization).";

    // log level

    QString log_level = xml.value("system.settings.logLevel", "debug").toLower();

    if (log_level=="debug") setLogLevel(0);

    if (log_level=="info") setLogLevel(1);

    if (log_level=="warning") setLogLevel(2);

    if (log_level=="error") setLogLevel(3);

    // snag dynamics / soil model enabled? (info used during setup of world)

    changeSettings().carbonCycleEnabled = xml.valueBool("model.settings.carbonCycleEnabled", false);

    // class size of snag classes

    Snag::setupThresholds(xml.valueDouble("model.settings.soil.swdDBHClass12"),

                          xml.valueDouble("model.settings.soil.swdDBHClass23"));

    // setup of modules

    if (mModules)

        delete mModules;

    mModules = new Modules();

    setupSpace();

    if (mRU.isEmpty())

        throw IException("Setup of Model: no resource units present!");

    // (3) additional issues

    // (3.1) load javascript code into the engine

    ScriptGlobal::loadScript(g->path(xml.value("system.javascript.fileName"),"script"));

    // (3.2) setup of regeneration

    changeSettings().regenerationEnabled = xml.valueBool("model.settings.regenerationEnabled", false);

    if (settings().regenerationEnabled) {

        foreach(SpeciesSet *ss, mSpeciesSets)

            ss->setupRegeneration();

    }

    Sapling::setRecruitmentVariation(xml.valueDouble("model.settings.seedDispersal.recruitmentDimensionVariation",0.1));

    // (3.3) management

    bool use_abe = xml.valueBool("model.management.abeEnabled");

    if (use_abe) {

        // use the agent based forest management engine

        mABEManagement = new AMIE::ForestManagementEngine();

        // setup of ABE after loading of trees.

    } else {

        // use the standard management

        QString mgmtFile = xml.value("model.management.file");

        if (!mgmtFile.isEmpty() && xml.valueBool("model.management.enabled")) {

            mManagement = new Management();

            QString path = GlobalSettings::instance()->path(mgmtFile, "script");

            mManagement->loadScript(path);

            qDebug() << "setup management using script" << path;

        }

    }

}

ResourceUnit *Model::ru(QPointF coord)

{

    if (!mRUmap.isEmpty() && mRUmap.coordValid(coord))

        return mRUmap.valueAt(coord);

    return ru(); // default RU if there is only one

}

void Model::initOutputDatabase()

{

    GlobalSettings *g = GlobalSettings::instance();

    QString dbPath = g->path(g->settings().value("system.database.out"), "output");

    // create run-metadata

    int maxid = SqlHelper::queryValue("select max(id) from runs", g->dbin()).toInt();

    maxid++;

    QString timestamp = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss");

    SqlHelper::executeSql(QString("insert into runs (id, timestamp) values (%1, '%2')").arg(maxid).arg(timestamp), g->dbin());

    // replace path information

    dbPath.replace("$id$", QString::number(maxid));

    dbPath.replace("$date$", timestamp);

    // setup final path

   g->setupDatabaseConnection("out", dbPath, false);

}

/// multithreaded running function for the resource unit level establishment

ResourceUnit *nc_sapling_growth_establishment(ResourceUnit *unit)

{

    try {

        // define a height map for the current resource unit on the stack

        float sapling_map[cPxPerRU*cPxPerRU];

        // set the map and initialize it:

        unit->setSaplingHeightMap(sapling_map);

        // (1) calculate the growth of (already established) saplings (populate sapling map)

        QList<ResourceUnitSpecies*>::const_iterator rus;

        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)

            (*rus)->calclulateSaplingGrowth();

        // (2) calculate the establishment probabilities of new saplings

        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)

            (*rus)->calculateEstablishment();

    } catch (const IException& e) {

        GlobalSettings::instance()->controller()->throwError(e.message());

    }

    unit->setSaplingHeightMap(0); // invalidate again

    return unit;

}

/// multithreaded execution of the carbon cycle routine

ResourceUnit *nc_carbonCycle(ResourceUnit *unit)

{

    try {

        // (1) do calculations on snag dynamics for the resource unit

        unit->calculateCarbonCycle();

        // (2) do the soil carbon and nitrogen dynamics calculations (ICBM/2N)

    } catch (const IException& e) {

        GlobalSettings::instance()->controller()->throwError(e.message());

    }

    return unit;

}

/// beforeRun performs several steps before the models starts running.

/// inter alia: * setup of the stands

///             * setup of the climates

void Model::beforeRun()

{

    // setup outputs

    // setup output database

    if (GlobalSettings::instance()->dbout().isOpen())

        GlobalSettings::instance()->dbout().close();

    initOutputDatabase();

    GlobalSettings::instance()->outputManager()->setup();

    GlobalSettings::instance()->clearDebugLists();

    // initialize stands

    {

    DebugTimer loadtrees("load trees");

    StandLoader loader(this);

    loader.processInit();

    }

    // load climate

    {

    if (logLevelDebug()) qDebug() << "attempting to load climate..." ;

    DebugTimer loadclim("load climate");

    foreach(Climate *c, mClimates)

        if (!c->isSetup())

            c->setup();

    }

    { DebugTimer loadinit("load standstatistics");

    if (logLevelDebug()) qDebug() << "attempting to calculate initial stand statistics (incl. apply and read pattern)..." ;

    Tree::setGrid(mGrid, mHeightGrid);

    // debugCheckAllTrees(); // introduced for debugging session (2012-04-06)

    applyPattern();

    readPattern();

    // force the compilation of initial stand statistics

    createStandStatistics();

    }

    // initalization of ABE

    if (mABEManagement) {

        DebugTimer t("ABESetup");

        mABEManagement->setup();

    }

    // outputs to create with inital state (without any growth) are called here:

    GlobalSettings::instance()->outputManager()->execute("stand"); // year=0

    GlobalSettings::instance()->outputManager()->execute("landscape"); // year=0

    GlobalSettings::instance()->outputManager()->execute("sapling"); // year=0

    GlobalSettings::instance()->outputManager()->execute("tree"); // year=0

    GlobalSettings::instance()->outputManager()->execute("dynamicstand"); // year=0

    GlobalSettings::instance()->setCurrentYear(1); // set to first year

}

/** Main model runner.

  The sequence of actions is as follows:

  (1) Load the climate of the new year

  (2) Reset statistics for resource unit as well as for dead/managed trees

  (3) Invoke Management.

  (4) *after* that, calculate Light patterns

  (5) 3PG on stand level, tree growth. Clear stand-statistcs before they are filled by single-tree-growth. calculate water cycle (with LAIs before management)

  (6) execute Regeneration

  (7) invoke disturbance modules

  (8) calculate statistics for the year

  (9) write database outputs

  */

void Model::runYear()

{

    DebugTimer t("Model::runYear()");

    GlobalSettings::instance()->systemStatistics()->reset();

    RandomGenerator::checkGenerator(); // see if we need to generate new numbers...

    // initalization at start of year for external modules

    mModules->yearBegin();

    // execute scheduled events for the current year

    if (mTimeEvents)

        mTimeEvents->run();

    // load the next year of the climate database

    foreach(Climate *c, mClimates)

        c->nextYear();

    // reset statistics

    foreach(ResourceUnit *ru, mRU)

        ru->newYear();

    foreach(SpeciesSet *set, mSpeciesSets)

        set->newYear();

    // management classic

    if (mManagement) {

        DebugTimer t("management");

        mManagement->run();

        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();

    }

    // ... or ABE (the agent based variant)

    if (mABEManagement) {

        DebugTimer t("ABERun");

        mABEManagement->run();

        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();

    }

    // process a cycle of individual growth

    applyPattern(); // create Light Influence Patterns

    readPattern(); // readout light state of individual trees

    grow(); // let the trees grow (growth on stand-level, tree-level, mortality)

    // regeneration

    if (settings().regenerationEnabled) {

        // seed dispersal

        DebugTimer tseed("Regeneration and Establishment");

        foreach(SpeciesSet *set, mSpeciesSets)

            set->regeneration(); // parallel execution for each species set

        GlobalSettings::instance()->systemStatistics()->tSeedDistribution+=tseed.elapsed();

        // establishment

        { DebugTimer t("saplingGrowthEstablishment");

        executePerResourceUnit( nc_sapling_growth_establishment, false /* true: force single thraeded operation */);

        GlobalSettings::instance()->systemStatistics()->tSaplingAndEstablishment+=t.elapsed();

        }

    }

    // calculate soil / snag dynamics

    if (settings().carbonCycleEnabled) {