#!/usr/bin/env python
# encoding: utf-8
"""
Title: Error models for simulating forest inventory errors

*inventory_error_models.py*

Each model definition must be of the form:
    def model_name(arg, oind)
in which
    arg - contains all the data being passed between the simulator
          and the prediction model
    oind - is the index number for the object for which the prediction is
           being made

The arg data object has the following structure:
    arg.variables[variable_index, object_index]
        - the values of the explaining variables of the model. 
          The number and order of these must match with the variable 
          declarations in the variables section of the model definition in
          the model xml file
    arg.parameters[parameter_index]
        - the model parameter values. The order in which these are used must 
          match the order given in the parameters section of the model xml
          file
    arg.num_of_res_objs[object_index]
      - number of result objects. This is always 1 unless the prediction
        model really creates new objects; e.g. in stem distribution models
        each diameter or height class is a new object
    arg.mem[variable_index, result_index]
      - stores the model results. Again the results must be stored here
        in the same order as defined in the result/variables section of
        the model xml file
    arg.errors[object_index] 
        - for storing possible error messages. Before actually 
          evaluating the model, you should test for error situations and
          set error messages detailing why the model can't be evalueated and
          then exit
When getting or setting the attribute values of arg, oind is always used for
the object_index; e.g. the value of the first explaining variable of the 
model is accessed by arg.variables[0][oind]. As can be seen from the example
indices always start from 0.

The model must return an integer as its return value. 1 if there were no
errors and 0 if errors occurred.

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

import math
import numpy


def Normal_inventory_error_BA_Hdom(arg, oind):
    """Generate normally distributed forest inventory errors to attributes
    BA and H_dom
    
    e_tilde = N(e_mean, e_stdev)
    """
    ret = 1
    BA = arg.variables[0, oind]
    H_dom = arg.variables[1, oind]
    bias = arg.parameters[0]
    stdev = arg.parameters[1]
    if stdev > 0.0:
        e = numpy.random.normal(bias, stdev, 2)
        e_BA, e_H_dom = e[0], e[1]
    else:
        e_BA, e_H_dom = 0.0, 0.0
    arg.mem[0, oind] = BA + BA * (e_BA / 100.0)
    arg.mem[1, oind] = H_dom + H_dom * (e_H_dom / 100.0)
    return ret

def Normal_inventory_error_d(arg, oind):
    """Generate normally distributed forest inventory errors to diameter
    
    e_tilde = N(e_mean, e_stdev)
    """
    ret = 1
    d = arg.variables[0, oind]
    bias = arg.parameters[0]
    stdev = arg.parameters[1]
    if stdev > 0.0:
        e_d = numpy.random.normal(bias, stdev, 1)[0]
    else:
        e_d = 0.0
    arg.mem[0, oind] = d + d * (e_d / 100.0)
    return ret

def Normal_inventory_error_h(arg, oind):
    """Generate normally distributed forest inventory errors to height
    
    e_tilde = N(e_mean, e_stdev)
    """
    ret = 1
    h = arg.variables[0, oind]
    bias = arg.parameters[0]
    stdev = arg.parameters[1]
    if stdev > 0.0:
        e_h = numpy.random.normal(bias, stdev, 1)[0]
    else:
        e_h = 0.0
    arg.mem[0, oind] = h + h * (e_h / 100.0)
    return ret

def Tree_extraction_error(arg, oind):
    """Generate random tree extraction error, ie. remove a tree with
    a given probability in order to simulate forest inventory omission
    errors.
    
    Tree omission probabilities were estimated from Evo LIDAR dataset.
    """
    ret = 1
    h = arg.variables[0, oind]  # tree height
    H_gM = arg.variables[1, oind]  # stand mean height
    prob_model = arg.parameters[0]  # probability model no
    
    # calculate tree's relative height and the extraction probability
    h_rel = h / H_gM
    if prob_model == 1:
        z = -2.255 + 3.305 * h_rel
    elif prob_model == 2:
        z = -1.2732 + 2.4174 * h_rel
    elif prob_model == 3:
        z = -1.0 + 4.0 * h_rel
    elif prob_model == 4:
        z = 0.5 + 4.0 * h_rel
    else:
        z = 4.0 + 4.0 * h_rel
    p = 1.0 / (1.0 + math.exp(-z))
    
    # draw a random number x from uniform distribution and if x <= p,
    # tree is extracted (=found from the RS data)
    x = numpy.random.uniform(low=0, high=1)
    if x <= p:
        omitted = 0
    else:
        omitted = 1
    arg.mem[0, oind] = omitted
    
    return ret


def main():
    pass

if __name__ == '__main__':
    main()
