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

*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 Basal_area_error_coeff(arg, oind):
    """
    Method: Basal_area_error_coeff
    
    Calculates random error coefficient to basal area estimate
    """
    ret = 1
    # error component u depicts the between stands variation in growth
    u = arg.variables[0, oind]
    age = arg.variables[1, oind]
    timestep = arg.variables[2, oind]
    
    if age > 10 and age <= 60:
        ya = 0.8
        yb = 0.1
        xa = 10
        xb = 60
        y = ya + ((age - xa) * (yb - ya)) / (xb - xa)
        u = u * y
    elif age > 60:
        u = u * 0.1
    
    # generate Gaussian random error component e depicting annual variation
    # in actual growth
    bias = 0.0
    std = (arg.parameters[0] ** 0.5) * (timestep / 5.0)
    e = numpy.random.normal(bias, std, 1)[0]
    
    # data back to SIMO
    arg.num_of_res_objs[oind] = 1
    arg.mem[0, oind] = arg.rect_factor * (u + e)
    return ret

def Dominant_height_error_coeff(arg, oind):
    """
    Method: Dominant_height_error_coeff
    
    Calculates random error coefficient to dominant height estimate
    """
    ret = 1
    # error component u depicts the between stands variation in growth
    u = arg.variables[0, oind]
    age = arg.variables[1, oind]
    timestep = arg.variables[2, oind]
    
    if age > 10 and age <= 60:
        ya = 0.8
        yb = 0.1
        xa = 10
        xb = 60
        y = ya + ((age - xa) * (yb - ya)) / (xb - xa)
        u = u * y
    elif age > 60:
        u = u * 0.1
    
    # generate Gaussian random error component e depicting annual variation
    # in actual growth
    bias = 0.0
    std = (arg.parameters[0] ** 0.5) * (timestep / 5.0)
    e = numpy.random.normal(bias, std, 1)[0]
    
    # data back to SIMO
    arg.num_of_res_objs[oind] = 1
    arg.mem[0, oind] = arg.rect_factor * (u + e)
    return ret

def Age_error_coeff_spruce(arg, oind):
    """
    Method: Age_error_coeff
    
    Calculates random error coefficient to age estimate
    """
    ret = 1
    # error component u depicts the between stands variation in growth
    u = arg.variables[0, oind]
    age = arg.variables[1, oind]
    timestep = arg.variables[2, oind]
    
    if age > 20 and age <= 80:
        ya = 0.9
        yb = 0.2
        xa = 20
        xb = 80
        y = ya + ((age - xa) * (yb - ya)) / (xb - xa)
        u = u * y
    elif age > 80:
        u = u * 0.2
    
    # generate Gaussian random error component e depicting annual variation
    # in actual growth
    bias = 0.0
    std = (arg.parameters[0] ** 0.5) * (timestep / 5.0)
    e = numpy.random.normal(bias, std, 1)[0]
    
    # data back to SIMO
    arg.num_of_res_objs[oind] = 1
    arg.mem[0, oind] = arg.rect_factor * (u + e)
    return ret

def Age_error_coeff_birch(arg, oind):
    """
    Method: Age_error_coeff
    
    Calculates random error coefficient to age estimate
    """
    ret = 1
    # error component u depicts the between stands variation in growth
    u = arg.variables[0, oind]
    age = arg.variables[1, oind]
    timestep = arg.variables[2, oind]
    
    if age > 20 and age <= 80:
        ya = 0.9
        yb = 0.2
        xa = 20
        xb = 80
        y = ya + ((age - xa) * (yb - ya)) / (xb - xa)
        u = u * y
    elif age > 80:
        u = u * 0.2
    
    # generate Gaussian random error component e depicting annual variation
    # in actual growth
    bias = 0.0
    std = (arg.parameters[0] ** 0.5) * (timestep / 5.0)
    e = numpy.random.normal(bias, std, 1)[0]
    
    # data back to SIMO
    arg.num_of_res_objs[oind] = 1
    arg.mem[0, oind] = arg.rect_factor * (u + e)
    return ret

def inventory_check(arg, oind):
    """
    Model for launching new inventory that nullifies the growth model errors.
    """
    ret = 1
    # model variables
    inventory_step = arg.variables[0, oind]
    time_step = arg.variables[1, oind]
    time_since_inventory = arg.variables[2, oind]
    # do inventory if time_since_inventory is equal to or greater than
    # inventory step
    do_inventory = 0
    if time_since_inventory >= inventory_step:
        time_since_inventory = 0
        do_inventory = 1
    # update time_since_inventory
    time_since_inventory += time_step
    # data back to SIMO
    arg.num_of_res_objs[oind] = 1
    arg.mem[0, oind] = time_since_inventory
    arg.mem[1, oind] = do_inventory
    return ret


def initial_growth_projection_errors(arg, oind):
    """Generate initial multinormal growth error coefficients for stand level
    growth models
    """
    ret = 1
    timestep = arg.variables[0, oind]
    f = timestep / 5.0
    E = _initial_random_error(1)
    arg.num_of_res_objs[oind] = 1
    arg.mem[0, oind] = E[0,0] * f
    arg.mem[1, oind] = E[0,1] * f
    arg.mem[2, oind] = E[0,2]
    arg.mem[3, oind] = E[0,3]
    return ret


def _initial_random_error(N):
   """Generate initial multinormal growth error coefficients for stand level
   growth models
   """
   sigma = numpy.array([
       [0.003600000, 0.002050122, -0.004285193, -0.004628952],
       [0.002050122, 0.004670000 ,-0.004880649, -0.005272176],
       [-0.004285193, -0.004880649, 0.007970000, 0.006887484],
       [-0.004628952, -0.005272176, 0.006887484, 0.009300000]],
       dtype=float)
   mu = numpy.zeros(4, dtype=float)
   return numpy.random.multivariate_normal(mu, sigma, N)


def main():
    pass

if __name__ == '__main__':
    main()
