text:
#=
Functions for biological rates :
     - metabolic rate
     - growth rate
     - handling time
     - attack rate

When temperature is not included, the functions :
    - no_effect_x for metabolic rate
    - no_effect_r for growth rate
    - no_effect_handlingt for handling time
    - no_effect_attackr for attack rate

return the default values for temperature independent rates.

When temperature is included, they are 4 different functions of temperature dependence :

1) Extended Eppley function:
    - extended_eppley_r for growth rate
    - extended_eppley_x for metabolic rate
2) Exponential Boltzmann-Arrhenius function:
    - exponential_BA_r for growth rate
    - exponential_BA_x for metabolic rate
    - exponential_BA_attackr for attack rate
	- exponential_BA_handlingt for handling time
3) Extended Boltzmann-Arrhenius function (Johnson-Lewin):
    - extended_BA_r for growth rate
    - extended_BA_x for metabolic rate
    - extended_BA_attackr for attack rate
4) Gaussian (inverted Gaussian) function
    - gaussian_r for growth rate
    - gaussian_x for metabolic rate
    - gaussian_attackr for attack rate
	- gaussian_handlingt for handling time

These functions are called in biological_rates.jl, where there are wrapped in one function for each temperature dependence functions:
 	- NoEffectTemperature for the functions with no effect of temperature
	- ExtendedEppley for the extended Eppley functions
	- ExponentialBA for the exponential Boltzmann-Arrhenius functions
	- ExtendedBA for the extended Boltzmann-Arrhenius functions
	- Gaussian for gaussian (inverted Gaussian) functions

Which are then used in model_parameters, as : growthrate = ExtendedEppley(:r), where the rate is specified as a key argument.
This allows to directly specify the function to use for each rate and the set of parameters associated with this function within the parameters.
See the documentation of biological_rates.jl for more details.

The functions of temperature dependence can be called with the default parameters provided :
	- ex: extended_eppley_r()
or with specified parameters :
	- ex : extended_eppley_r(passed_temp_parameters = T_parameters)
where T_parameters are the parameters associated with the extended_eppley_r function.

Internally the function extended_eppley_r (or whichever is chosen) takes 3 arguments:

- species body mass (standardized by the smaller producer species),
- the temperature T in Kelvin,
- the set of parameters p.

It is called in model_parameters as r = growthrate(1.0, 295.0, p).

(All rates scale with bodymass to an exponent β, and parameters of the functional response also scale with resource bodymass to an exponent β_resource.)

model_parameters still returns a Dict, containing:

- x, the metabolic rate (unchanged)
- r, the intrinsic producer's growth rate, which is now specified for each species
- y, the maximum consumption rate, which is 1/ht, with ht being the handling time
- Γ, the half saturation density, which is 1/(ar * ht), with ar being the attack rate
=#

"""
**No effect of temperature on metabolic rate**

This function is by default called as an argument in `NoEffectTemperature` to define a temperature independent metabolic rate.
It returns the default values of metabolic rate when there is no effect of temperature (see Delmas et al 2017).
The function accounts for different metabolic rate according to the type of species (vertebrate, invertebrate, producer)

| Parameter       | Meaning             | Default values | Reference                       |
|:----------------|:--------------------|:---------------|:--------------------------------|
| a_vertebrate    | allometric constant | 0.88           | Brose, Williams & Martinez 2006 |
| a_invertebrate  | allometric constant | 0.3141         | Brose, Williams & Martinez 2006 |
| a_producer      | allometric constant | 0.138          | Brose, Williams & Martinez 2006 |

Internally the function takes 3 arguments:

- species body mass (standardized by the smaller producer species),
- the temperature T,
- the set of parameters p.

"""
function no_effect_x(default_temp_parameters = (a_vertebrate = 0.88, a_invertebrate = 0.3141, a_producer = 0.138); passed_temp_parameters...)
    if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end
    return (bodymass, T, p) ->  (temperature_param.a_vertebrate .* (p[:vertebrates] .& .!p[:is_producer]) + temperature_param.a_invertebrate * (.!p[:vertebrates] .& .!p[:is_producer]) + temperature_param.a_producer .* p[:is_producer]) .* (bodymass.^-0.25)
end

"""
**No effect of temperature on growth rate**

This function is by default called as an argument in `NoEffectTemperature` to define a temperature independent growth rate.
It returns a default value for growth rate when there is no effect of temperature (see Delmas et al 2017).

| Parameter       | Meaning             | Default values | Reference                       |
|:----------------|:--------------------|:---------------|:--------------------------------|
| r               | growth rate         | 1              | Brose, Williams & Martinez 2006 |

Internally the function takes 3 arguments (unused in this case):

- species body mass (standardized by the smaller producer species),
- the temperature T,
- the set of parameters p.

"""
function no_effect_r(default_temp_parameters = (r = 1.0,); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end
    return (bodymass, T, p) -> repeat([temperature_param.r], size(p[:A], 1))
end

"""
**No effect of temperature on handling time**

This function is by default called as an argument in `NoEffectTemperature` to define a temperature independent handling time.
It returns a default value for handling time which is defined by 1/y, y being the maximum consumption rate.
The function accounts for different maximum consumption rate values according to the type of species (vertebrate or invertebrate).

| Parameter       | Meaning                                    | Default values | Reference                       |
|:----------------|:-------------------------------------------|:---------------|:--------------------------------|
| y_vertebrate    | maximum consumption rate for vertebrates   | 4.0            | Brose, Williams & Martinez 2006 |
| y_invertebrate  | maximum consumption rate for invertebrates | 8.0            | Brose, Williams & Martinez 2006 |

Internally the function takes 3 arguments (unused in this case):

- species body mass (standardized by the smaller producer species),
- the temperature T,
- the set of parameters p.

"""
function no_effect_handlingt(default_temp_parameters = (y_vertebrate = 4.0, y_invertebrate = 8.0); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end
    return (bodymass, T, p) ->  1 ./ (temperature_param.y_vertebrate .* (p[:vertebrates] .& .!p[:is_producer]) + temperature_param.y_invertebrate * (.!p[:vertebrates] .& .!p[:is_producer]))
end

"""
**No effect of temperature on attack rate**

This function is by default called as an argument in `NoEffectTemperature` to define a temperature independent attack rate.
It returns a default value for handling time which is defined by 1/(Γ*h), Γ being the half saturation density.

| Parameter  | Meaning                   | Default values | Reference                       |
|:-----------|:--------------------------|:---------------|:--------------------------------|
| Γ          | half-saturation density   | 0.5            | Brose, Williams & Martinez 2006 |

Internally the function takes 3 arguments (unused in this case):

- species body mass (standardized by the smaller producer species),
- the temperature T,
- the set of parameters p.

"""
function no_effect_attackr(default_temp_parameters = (Γ = 0.5,); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end
    return (bodymass, T, p) -> 1 ./ (temperature_param.Γ .* p[:ht])
end

"""
TODO
**Option 1 : Extended Eppley function for growth rate**

This function can be called with the keywords :growth, :growthrate or :r as an argument in `ExtendedEppley`, itself called in `model_parameters`, to define an extended Eppley function (Eppley 1972, Thomas et al. 2012) for producers growth rate.
Example : model_parameters(A, growthrate = ExtendedEppley(:r))

| Parameter       | Meaning                                          | Default values| Reference            |
|:----------------|:-------------------------------------------------|:--------------|:---------------------|
| maxrate_0       | Maximum rate at 273.15 degrees Kelvin            | 0.81          | Eppley 1972          |
| eppley_exponent | Exponential rate of increase                     | 0.0631        | Eppley 1972          |
| z		          | Location of the inflexion point of the function  | 298.15        | NA                   |
| range           | Thermal breadth                                  | 35            | NA                   |
| β               | Allometric exponent                              | -0.25         | Gillooly et al. 2002 |

Default values are given as an example.

The function can be called with the default parameters provided as an example:
	- extended_eppley_r()
Parameters can also be specified:
	- extended_eppley_r(passed_temp_parameters = (maxrate_0 = 0.5, eppley_exponent = 0.0631, z = 315.0, β = -0.25, range = 35))
"""
function extended_eppley_r(default_temp_parameters = (maxrate_0 = 0.81, eppley_exponent = 0.0631, z = 298.15, β = -0.25, range = 35); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

    z = temperature_param.z - 273.15

    return (bodymass, T, p) -> bodymass.^temperature_param.β .* temperature_param.maxrate_0 .* exp(temperature_param.eppley_exponent .* (T.-273.15)) * (1 .- (((T.-273.15) .- z) ./ (temperature_param.range./2)).^2)
end

"""
TODO
**Option 1 : Extended Eppley function for metabolic rate**

This function can be called with the keywords :metabolism, :x or :metabolicrate as an argument in `ExtendedEppley`, itself called in `model_parameters`, to define an extended Eppley function (Eppley 1972, Thomas et al. 2012) for metabolic rate.
Example : model_parameters(A, metabolicrate = ExtendedEppley(:x))

| Parameter                     | Meaning                                                                           | Default values| Reference            |
|:------------------------------|:----------------------------------------------------------------------------------|:--------------|:---------------------|
| maxrate_0_producer            | Maximum rate at 273.15 degrees Kelvin for producers                   | 0.81          | Eppley 1972          |
| maxrate_0_invertebrate        | Maximum rate at 273.15 degrees Kelvin for invertebrates               | 0.81          | Eppley 1972          |
| maxrate_0_vertebrate          | Maximum rate at 273.15 degrees Kelvin for vertebrates                 | 0.81          | Eppley 1972          |
| eppley_exponent_producer      | Exponential rate of increase for producers                            | 0.0631        | Eppley 1972          |
| eppley_exponent_invertebrate  | Exponential rate of increase for invertebrates                        | 0.0631        | Eppley 1972          |
| eppley_exponent_vertebrate    | Exponential rate of increase for vertebrates                          | 0.0631        | Eppley 1972          |
| z_producer  		            | Location of the inflexion point of the function for producers    		| 298.15        | NA                   |
| z_invertebrate                | Location of the inflexion point of the function for invertebrates		| 298.15        | NA                   |
| z_producer_vertebrate         | Location of the inflexion point of the function for vertebrates  		| 298.15        | NA                   |
| range_producer                | Thermal breadth for producers                                         | 35            | NA                   |
| range_invertebrate            | Thermal breadth for invertebrates                                     | 35            | NA                   |
| range_vertebrate              | Thermal breadth for vertebrates                                       | 35            | NA                   |
| β_producer                    | Allometric exponent for producers                                     | -0.25         | Gillooly et al. 2002 |
| β_invertebrate                | Allometric exponent for invertebrates                                 | -0.25         | Gillooly et al. 2002 |
| β_vertebrate                  | Allometric exponent for vertebrates                                   | -0.25         | Gillooly et al. 2002 |

Default values are given as an example (note that they are initially provided for phytoplankton growth rate in Eppley 1972)

The function can be called with the default parameters provided as an example:
	- extended_eppley_x()
Parameters can also be specified (for producers, invertebrates and vertebrates):
	- extended_eppley_x(passed_temp_parameters = (maxrate_0_producer = 0.81, maxrate_0_invertebrate = 0.81, maxrate_0_vertebrate = 0.81,
	                                     eppley_exponent_producer = 0.0631, eppley_exponent_invertebrate = 0.0631, eppley_exponent_vertebrate = 0.0631,
	                                     z_producer = 310.15, z_invertebrate = 310.15, z_vertebrate = 298.15,
	                                     range_producer = 35, range_invertebrate = 35, range_vertebrate = 35,
	                                     β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25))

"""
function extended_eppley_x(default_temp_parameters = (maxrate_0_producer = 0.81, maxrate_0_invertebrate = 0.81, maxrate_0_vertebrate = 0.81,
                                     eppley_exponent_producer = 0.0631, eppley_exponent_invertebrate = 0.0631, eppley_exponent_vertebrate = 0.0631,
                                     z_producer = 298.15, z_invertebrate = 298.15, z_vertebrate = 298.15,
                                     range_producer = 35, range_invertebrate = 35, range_vertebrate = 35,
                                     β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

    return (bodymass, T, p) -> for i in 1:1
                                    maxrate_0_all = temperature_param.maxrate_0_producer .* p[:is_producer] .+ temperature_param.maxrate_0_vertebrate .* p[:vertebrates] .+ temperature_param.maxrate_0_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                    eppley_exponent_all = temperature_param.eppley_exponent_producer .* p[:is_producer] .+ temperature_param.eppley_exponent_vertebrate .* p[:vertebrates] .+ temperature_param.eppley_exponent_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                    z_all = temperature_param.z_producer .* p[:is_producer] .+ temperature_param.z_vertebrate .* p[:vertebrates] .+ temperature_param.z_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                    z_all = z_all .- 273.15
                                    range_all = temperature_param.range_producer .* p[:is_producer] .+ temperature_param.range_vertebrate .* p[:vertebrates] .+ temperature_param.range_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                    β_all = temperature_param.β_producer .* p[:is_producer] .+ temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])

                                    return bodymass.^β_all .* maxrate_0_all .* exp.(eppley_exponent_all .* (T.-273.15)) .* (1 .- (((T.-273.15) .- z_all) ./ (range_all./2)).^2)
                                end
end

"""
**Option 2 : Exponential Boltzmann-Arrhenius function for growth rate**

This function can be called with the keywords :growth, :growthrate or :r as an argument in `ExponentialBA`, itself called in `model_parameters`, to define an exponential Boltzmann-Arrhénius function (Gillooly et al. 2001, Brown et al. 2004) for producers growth rate.
Example : model_parameters(A, growthrate = ExponentialBA(:r))

| Parameter         | Meaning                               | Default values  | Reference                             |
|:------------------|:--------------------------------------|:----------------|:--------------------------------------|
| norm_constant     | scaling coefficient                   | exp(-15.68)*4e6 | Ehnes et al. 2011, Binzer et al. 2012 |
| activation_energy | activation energy                     | -0.84           | Ehnes et al. 2011, Binzer et al. 2012 |
| T0                | normalization temperature (K)         | 293.15          | Binzer et al. 2012, Binzer et al. 2012|
| β                 | allometric exponent                   | -0.25           | Ehnes et al. 2011                     |
| k                 | Boltzmann norm_constant               | 8.617e-5        |                                       |

Default values are given as an example.

The function can be called with the default parameters:
	- exponential_BA_r()
Parameters can also be specified:
	- exponential_BA_r(passed_temp_parameters = (norm_constant = exp(-15.68)*4e6, activation_energy = -0.72, T0 = 293.15, β = -0.25))
"""
function exponential_BA_r(default_temp_parameters = (norm_constant = exp(-15.68)*4e6, activation_energy = -0.84, T0 = 293.15, β = -0.25); passed_temp_parameters...)
    k = 8.617e-5
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end
    return (bodymass, T, p) -> temperature_param.norm_constant .* (bodymass .^temperature_param.β) .* exp.(temperature_param.activation_energy .* (temperature_param.T0 .- T) ./ (k .* T .* temperature_param.T0))
end

"""
**Option 2 : Exponential Boltzmann-Arrhenius function for metabolic rate**

This function can be called with the keywords :metabolism, :x or :metabolicrate as an argument in `ExponentialBA`, itself called in `model_parameters`, to define an exponential Boltzmann-Arrhénius function (Gillooly et al. 2001, Brown et al. 2004) for metabolic rate.
Example : model_parameters(A, metabolicrate = ExponentialBA(:x))

| Parameter                      | Meaning                                         | Default values   | Reference                             |
|:-------------------------------|:------------------------------------------------|:-----------------|:--------------------------------------|
| norm_constant_producer         | scaling coefficient for producers               | exp(-16.54)*4e6  | Ehnes et al. 2011, Binzer et al. 2012 |
| norm_constant_invertebrate     | scaling coefficient for invertebrates           | exp(-16.54)*4e6  | Ehnes et al. 2011, Binzer et al. 2012 |
| norm_constant_vertebrate       | scaling coefficient for vertebrates             | exp(-16.54)*4e6  | Ehnes et al. 2011, Binzer et al. 2012 |
| activation_energy_producer     | activation energy for producers                 | -0.69            | Ehnes et al. 2011, Binzer et al. 2012 |
| activation_energy_invertebrate | activation energy for invertebrates             | -0.69            | Ehnes et al. 2011, Binzer et al. 2012 |
| activation_energy_vertebrate   | activation energy for vertebrates               | -0.69            | Ehnes et al. 2011, Binzer et al. 2012 |
| T0_producer                    | normalization temperature (K) for producers     | 293.15           | Binzer et al. 2012					  |
| T0_invertebrate                | normalization temperature (K) for invertebrates | 293.15           | Binzer et al. 2012					  |
| T0_vertebrate                  | normalization temperature (K) for vertebrates   | 293.15           | Binzer et al. 2012					  |
| β_producer                     | allometric exponent for producers               | -0.31            | Ehnes et al. 2011                     |
| β_invertebrate                 | allometric exponent for invertebrates           | -0.31            | Ehnes et al. 2011                     |
| β_vertebrate                   | allometric exponent for vertebrates             | -0.31            | Ehnes et al. 2011                     |

Default values are given as an example.

The function can be called with the default parameters:
	- exponential_BA_x()
Parameters can also be specified:
	- exponential_BA_x(passed_temp_parameters = (norm_constant_producer = exp(-16.54)*4e6, norm_constant_invertebrate = exp(-16.54)*4e6, norm_constant_vertebrate = exp(-16.54)*4e6,
	                                   activation_energy_producer = -0.69, activation_energy_invertebrate = -0.69, activation_energy_vertebrate = -0.69,
	                                   T0_producer = 300.15, T0_invertebrate = 300.15, T0_vertebrate = 293.15,
	                                   β_producer = -0.31, β_invertebrate = -0.31, β_vertebrate = -0.31))
"""
function exponential_BA_x(default_temp_parameters = (norm_constant_producer = exp(-16.54)*4e6, norm_constant_invertebrate = exp(-16.54)*4e6, norm_constant_vertebrate = exp(-16.54)*4e6,
                                   activation_energy_producer = -0.69, activation_energy_invertebrate = -0.69, activation_energy_vertebrate = -0.69,
                                   T0_producer = 293.15, T0_invertebrate = 293.15, T0_vertebrate = 293.15,
                                   β_producer = -0.31, β_invertebrate = -0.31, β_vertebrate = -0.31); passed_temp_parameters...)
    k=8.617e-5
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

    return (bodymass, T, p) -> for i in 1:1
                                norm_constant_all = temperature_param.norm_constant_producer .* p[:is_producer] .+ temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                activation_energy_all = temperature_param.activation_energy_producer .* p[:is_producer] .+ temperature_param.activation_energy_vertebrate .* p[:vertebrates] .+ temperature_param.activation_energy_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                T0_all = temperature_param.T0_producer .* p[:is_producer] .+ temperature_param.T0_vertebrate .* p[:vertebrates] .+ temperature_param.T0_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                β_all = temperature_param.β_producer .* p[:is_producer] .+ temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])

                                return norm_constant_all .* (bodymass .^β_all) .* exp.(activation_energy_all .* (T0_all .- T ) ./ (k .* T .* T0_all))
                            end
end

"""
**Option 2 : Exponential Boltzmann-Arrhenius function for attack rate**

This function can be called with the keywords :attackrate as an argument in `ExponentialBA`, itself called in `model_parameters`, to define an exponential Boltzmann-Arrhénius function (Gillooly et al. 2001, Brown et al. 2004) for attack rate.
Example : model_parameters(A, attackrate = ExponentialBA(:attackrate))

| Parameter                      | Meaning                                         | Default values | Reference                             |
|:-------------------------------|:------------------------------------------------|:---------------|:--------------------------------------|
| norm_constant_invertebrate     | scaling coefficient for invertebrate            | exp(-13.1)*4e6 | Rall et al. 2012, Binzer et al. 2015  |
| norm_constant_vertebrate       | scaling coefficient for vertebrate              | exp(-13.1)*4e6 | Rall et al. 2012, Binzer et al. 2015  |
| activation_energy_invertebrate | activation energy for invertebrates             | -0.38          | Rall et al. 2012, Binzer et al. 2015  |
| activation_energy_vertebrate   | activation energy for vertebrates               | -0.38          | Rall et al. 2012, Binzer et al. 2015  |
| T0_invertebrate                | normalization temperature (K) for invertebrates | 293.15         | Rall et al. 2012, Binzer et al. 2015  |
| T0_vertebrate                  | normalization temperature (K) for vertebrates   | 293.15         | Rall et al. 2012, Binzer et al. 2015  |
| β_producer                     | allometric exponent for invertebrate            | 0.25           | Rall et al. 2012, Binzer et al. 2015  |
| β_invertebrate                 | allometric exponent for invertebrate            | -0.8           | Rall et al. 2012, Binzer et al. 2015  |
| β_vertebrate                   | allometric exponent for vertebrates             | -0.8           | Rall et al. 2012, Binzer et al. 2015  |

Default values are given as an example.

The function can be called with the default parameters:
	- exponential_BA_attackr()
Parameters can also be specified:
	- exponential_BA_attackr(passed_temp_parameters = (norm_constant_vertebrate = exp(-13.1)*4e6, norm_constant_invertebrate = exp(-13.1)*4e6,
	                             activation_energy_vertebrate = -0.38, activation_energy_invertebrate = -0.38,
	                             T0_vertebrate = 293.15, T0_invertebrate = 293.15,
	                             β_producer = 0.25, β_vertebrate = -0.8, β_invertebrate = 0.8))
"""
function exponential_BA_attackr(default_temp_parameters = (norm_constant_vertebrate = exp(-13.1)*4e6, norm_constant_invertebrate = exp(-13.1)*4e6,
											  activation_energy_vertebrate = -0.38, activation_energy_invertebrate = -0.38,
											  T0_vertebrate = 293.15, T0_invertebrate = 293.15,
											  β_producer = 0.25, β_vertebrate = -0.8, β_invertebrate = -0.8); passed_temp_parameters...)
    k=8.617e-5
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

    return (bodymass, T, p) -> for i in 1:1
                                norm_constant_all = temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                activation_energy_all = temperature_param.activation_energy_vertebrate .* p[:vertebrates] .+ temperature_param.activation_energy_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                T0_all = temperature_param.T0_vertebrate .* p[:vertebrates] .+ temperature_param.T0_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                β_consumer = temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])

                                β_resource = zeros(Float64,(p[:S], p[:S]))
                                for consumer in 1:p[:S]
                                for resource in 1:p[:S]
                                  if p[:A][consumer, resource] == 1
                                    if p[:is_producer][resource]
                                        β_resource[consumer, resource] = temperature_param.β_producer
                                    elseif p[:vertebrates][resource]
                                        β_resource[consumer, resource] = temperature_param.β_vertebrate
                                    else
                                        β_resource[consumer, resource] = temperature_param.β_invertebrate
                                    end
                                 end
                                end
                                end
                                rate = norm_constant_all .* (bodymass .^β_consumer) .* (bodymass' .^β_resource) .* exp.(activation_energy_all .* (T0_all .- T) ./ (k .* T .* T0_all))
                                rate[isnan.(rate)] .= 0
                            return rate
                        end
end

"""
**Option 2 : Exponential Boltzmann-Arrhenius function for handling time**

This function can be called with the keywords :handlingtime as an argument in `ExponentialBA`, itself called in `model_parameters`, to define an exponential Boltzmann-Arrhénius function (Gillooly et al. 2001, Brown et al. 2004) for handling time.
Example : model_parameters(A, handlingtime = ExponentialBA(:handlingtime))

| Parameter                      | Meaning                                         | Default values | Reference                             |
|:-------------------------------|:------------------------------------------------|:---------------|:--------------------------------------|
| norm_constant_invertebrate     | scaling coefficient for invertebrate            | exp(9.66)*4e6  | Rall et al. 2012, Binzer et al. 2015  |
| norm_constant_vertebrate       | scaling coefficient for vertebrate              | exp(9.66)*4e6  | Rall et al. 2012, Binzer et al. 2015  |
| activation_energy_invertebrate | activation energy for invertebrates             | 0.26  	        | Rall et al. 2012, Binzer et al. 2015  |
| activation_energy_vertebrate   | activation energy for vertebrates               | 0.26           | Rall et al. 2012, Binzer et al. 2015  |
| T0_invertebrate                | normalization temperature (K) for invertebrates | 293.15         | Rall et al. 2012, Binzer et al. 2015  |
| T0_vertebrate                  | normalization temperature (K) for vertebrates   | 293.15         | Rall et al. 2012, Binzer et al. 2015  |
| β_producer                     | allometric exponent for producer                | -0.45          | Rall et al. 2012, Binzer et al. 2015  |
| β_invertebrate                 | allometric exponent for invertebrate            | 0.47           | Rall et al. 2012, Binzer et al. 2015  |
| β_vertebrate                   | allometric exponent for vertebrates             | 0.47           | Rall et al. 2012, Binzer et al. 2015  |

Default values are given as an example.

The function can be called with the default parameters:
	- exponential_BA_handlingt()
Parameters can also be specified:
	- exponential_BA_handlingt(passed_temp_parameters = (norm_constant_vertebrate = exp(9.66)*4e6, norm_constant_invertebrate = exp(9.66)*4e6,
                                              activation_energy_vertebrate = 0.26, activation_energy_invertebrate = 0.26,
                                              T0_vertebrate = 293.15, T0_invertebrate = 293.15,
                                              β_producer = -0.45, β_vertebrate = 0.47, β_invertebrate = 0.47))
"""
function exponential_BA_handlingt(default_temp_parameters = (norm_constant_vertebrate = exp(9.66)*4e6, norm_constant_invertebrate = exp(9.66)*4e6,
											  activation_energy_vertebrate = 0.26, activation_energy_invertebrate = 0.26,
											  T0_vertebrate = 293.15, T0_invertebrate = 293.15,
											  β_producer = -0.45, β_vertebrate = 0.47, β_invertebrate = 0.47); passed_temp_parameters...)
    k=8.617e-5
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

    return (bodymass, T, p) -> for i in 1:1
                                norm_constant_all = temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                activation_energy_all = temperature_param.activation_energy_vertebrate .* p[:vertebrates] .+ temperature_param.activation_energy_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                T0_all = temperature_param.T0_vertebrate .* p[:vertebrates] .+ temperature_param.T0_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                β_consumer = temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])

                                β_resource = zeros(Float64,(p[:S], p[:S]))
                                for consumer in 1:p[:S]
                                for resource in 1:p[:S]
                                  if p[:A][consumer, resource] == 1
                                    if p[:is_producer][resource]
                                        β_resource[consumer, resource] = temperature_param.β_producer
                                    elseif p[:vertebrates][resource]
                                        β_resource[consumer, resource] = temperature_param.β_vertebrate
                                    else
                                        β_resource[consumer, resource] = temperature_param.β_invertebrate
                                    end
                                 end
                                end
                                end
                                rate = norm_constant_all .* (bodymass .^β_consumer) .* (bodymass' .^β_resource) .* exp.(activation_energy_all .* (T0_all .- T) ./ (k .* T .* T0_all))
                                rate[isnan.(rate)] .= 0
                            return rate
                        end
end

"""
**Option 3 : Extended Boltzmann-Arrhenius function for growth rate**

This function can be called with the keywords :growth, :growthrate or :r as an argument in `ExponentialBA`, itself called in `model_parameters`, to define an extended Boltzmann-Arrhénius function (Pawar et al 2015, Gillooly et al. 2001, Brown et al. 2004) for producers growth rate.
Example : model_parameters(A, growthrate = ExtendedBA(:r))

| Parameter          | Meaning                                               | Default values | Reference            |
|:-------------------|:------------------------------------------------------|:---------------|----------------------|
| norm_constant      | scaling coefficient                                   | 1.8e9          | NA					 |
| activation_energy  | activation energy                                     | 0.53           | Dell et al 2011      |
| deactivation_energy| deactivation energy                                   | 1.15           | Dell et al 2011      |
| T_opt              | temperature at which trait value is maximal           | 298.15         | NA                   |
| β                  | allometric exponent                                   | -0.25          | Gillooly et al. 2002 |

Default values are given as an example.

The function can be called with the default parameters:
	- extended_BA_r()
Parameters can also be specified:
	- extended_BA_r(passed_temp_parameters = (norm_constant = 1.8e9, activation_energy = 0.53, deactivation_energy = 1.15, T_opt = 298.15, β = -0.25))
"""
function extended_BA_r(default_temp_parameters = (norm_constant = 1.8e9, activation_energy = 0.53, deactivation_energy = 1.15, T_opt = 298.15, β = -0.25); passed_temp_parameters...)
     if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end
	k = 8.617e-5 # Boltzmann constant
	Δenergy = temperature_param.deactivation_energy .- temperature_param.activation_energy

    return(bodymass, T, p) -> temperature_param.norm_constant .* bodymass .^(temperature_param.β) .* exp.(.-temperature_param.activation_energy ./ (k * T)) .* (1 ./ (1 + exp.(-1 / (k * T) .* (temperature_param.deactivation_energy .- (temperature_param.deactivation_energy ./ temperature_param.T_opt .+ k .* log(temperature_param.activation_energy ./ Δenergy)).*T))))
end

"""
**Option 3 : Extended Boltzmann-Arrhenius function for metabolic rate**

This function can be called with the keywords :metabolism, :x or :metabolicrate as an argument in `ExponentialBA`, itself called in `model_parameters`, to define an extended Boltzmann-Arrhénius function (Pawar et al 2015, Gillooly et al. 2001, Brown et al. 2004) for metabolic rate.
Example : model_parameters(A, metabolicrate = ExtendedBA(:x))

| Parameter                       | Meaning                                                       | Default values | Reference            |
|:--------------------------------|:--------------------------------------------------------------|:---------------|----------------------|
| norm_constant_producer          | scaling coefficient for producers		                      | 1.5e9          | NA                   |
| norm_constant_invertebrate      | scaling coefficient for invertebrates                         | 1.5e9          | NA                   |
| norm_constant_invertebrate      | scaling coefficient for vertebrates                           | 1.5e9          | NA                   |
| activation_energy_producer      | activation energy for producers                               | 0.53           | Dell et al 2011      |
| activation_energy_invertebrate  | activation energy for invertebrates                           | 0.53           | Dell et al 2011      |
| activation_energy_vertebrate    | activation energy for vertebrates                             | 0.53           | Dell et al 2011      |
| deactivation_energy_producer    | deactivation energy for producers                             | 1.15           | Dell et al 2011      |
| deactivation_energy_invertebrate| deactivation energy for invertebrates                         | 1.15           | Dell et al 2011      |
| deactivation_energy_vertebrate  | deactivation energy for vertebrates                           | 1.15           | Dell et al 2011      |
| T_opt_producer                  | temperature at which trait value is maximal for producers	  | 298.15         | NA                   |
| T_opt_invertebrate              | temperature at which trait value is maximal for invertebrates | 298.15         | NA                   |
| T_opt_vertebrate                | temperature at which trait value is maximal for vertebrates   | 298.15         | NA                   |
| β_producer                      | allometric exponent for producers                             | -0.25          | Gillooly et al. 2002 |
| β_invertebrate                  | allometric exponent for invertebrates                         | -0.25          | Gillooly et al. 2002 |
| β_vertebrate                    | allometric exponent for vertebrates                           | -0.25          | Gillooly et al. 2002 |

Default values are given as an example.

The function can be called with the default parameters:
	- extended_BA_x()
Parameters can also be specified:
	- extended_BA_x(passed_temp_parameters = (norm_constant_producer = 1.5e9, norm_constant_invertebrate = 1.5e9, norm_constant_vertebrate = 1.5e9,
	                                activation_energy_producer = 0.53, activation_energy_invertebrate = 0.53, activation_energy_vertebrate = 0.53,
	                                deactivation_energy_producer = 1.15, deactivation_energy_invertebrate = 1.15, deactivation_energy_vertebrate = 1.15,
	                                T_opt_producer = 298.15, T_opt_invertebrate = 298.15, T_opt_vertebrate = 298.15,
	                                β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25))
"""
function extended_BA_x(default_temp_parameters = (norm_constant_producer = 3e8, norm_constant_invertebrate = 3e8, norm_constant_vertebrate = 3e8,
                                activation_energy_producer = 0.53, activation_energy_invertebrate = 0.53, activation_energy_vertebrate = 0.53,
                                deactivation_energy_producer = 1.15, deactivation_energy_invertebrate = 1.15, deactivation_energy_vertebrate = 1.15,
                                T_opt_producer = 298.15, T_opt_invertebrate = 298.15, T_opt_vertebrate = 298.15,
                                β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25); passed_temp_parameters...)
     k = 8.617e-5 # Boltzmann constant
	 if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
 	else
	  temperature_param = default_temp_parameters
	end

     return(bodymass, T, p) -> for i in 1:1
                                 norm_constant_all = temperature_param.norm_constant_producer .* p[:is_producer] .+ temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 activation_energy_all = temperature_param.activation_energy_producer .* p[:is_producer] .+ temperature_param.activation_energy_vertebrate .* p[:vertebrates] .+ temperature_param.activation_energy_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 deactivation_energy_all = temperature_param.deactivation_energy_producer .* p[:is_producer] .+ temperature_param.deactivation_energy_vertebrate .* p[:vertebrates] .+ temperature_param.deactivation_energy_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 T_opt_all = temperature_param.T_opt_producer .* p[:is_producer] .+ temperature_param.T_opt_vertebrate .* p[:vertebrates] .+ temperature_param.T_opt_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 β_all = temperature_param.β_producer .* p[:is_producer] .+ temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 Δenergy = deactivation_energy_all .- activation_energy_all

                                 return  norm_constant_all .* bodymass .^ (β_all) .* exp.(.-activation_energy_all ./ (k * T)) .* (1 ./ (1 .+ exp.(-1 / (k * T) .* (deactivation_energy_all .- (deactivation_energy_all ./ T_opt_all .+ k .* log.(activation_energy_all ./ Δenergy)) .* T))))

                             end
end

"""
**Option 3 : Extended Boltzmann-Arrhenius function for attack rate**

This function can be called with the keywords :attackrate as an argument in `ExtendedBA`, itself called in `model_parameters`, to define an extended Boltzmann-Arrhénius function (Gillooly et al. 2001, Brown et al. 2004) for attack rate.
Example : model_parameters(A, attackrate = ExtendedBA(:attackrate))

| Parameter                       | Meaning                                                       | Default values | Reference            |
|:--------------------------------|:--------------------------------------------------------------|:---------------|----------------------|
| norm_constant_invertebrate      | scaling coefficient for invertebrates                         | 5e13           | Bideault et al 2019  |
| norm_constant_vertebrate        | scaling coefficient for vertebrates                           | 5e13           | Bideault et al 2019  |
| activation_energy_invertebrate  | activation energy for invertebrates                           | 0.8            | Dell et al 2011      |
| activation_energy_vertebrate    | activation energy for vertebrates                             | 0.8            | Dell et al 2011      |
| deactivation_energy_invertebrate| deactivation energy for invertebrates                         | 1.15           | Dell et al 2011      |
| deactivation_energy_vertebrate  | deactivation energy for vertebrates                           | 1.15           | Dell et al 2011      |
| T_opt_invertebrate              | temperature at which trait value is maximal for invertebrates | 298.15         | NA                   |
| T_opt_vertebrate                | temperature at which trait value is maximal for vertebrates   | 298.15         | NA                   |
| β_producer                      | allometric exponent for producers                             | 0.25           | Rall et al. 2012     |
| β_invertebrate                  | allometric exponent for invertebrates                         | 0.25           | Rall et al. 2012     |
| β_vertebrate                    | allometric exponent for vertebrates                           | 0.25           | Rall et al. 2012     |

Default values are given as an example.

The function can be called with the default parameters:
	- extended_BA_attackr()
Parameters can also be specified:
	- extended_BA_attackr(passed_temp_parameters = (norm_constant_invertebrate = 5e13, norm_constant_vertebrate = 5e13,
	                                   activation_energy_invertebrate = 0.8, activation_energy_vertebrate = 0.8,
	                                   deactivation_energy_invertebrate = 1.15, deactivation_energy_vertebrate = 1.15,
	                                   T_opt_invertebrate = 298.15, T_opt_vertebrate = 298.15,
	                                   β_producer = 0.25, β_invertebrate = 0.25, β_vertebrate = 0.25))
"""
function extended_BA_attackr(default_temp_parameters = (norm_constant_invertebrate = 5e13, norm_constant_vertebrate = 5e13,
                                   activation_energy_invertebrate = 0.8, activation_energy_vertebrate = 0.8,
                                   deactivation_energy_invertebrate = 1.15, deactivation_energy_vertebrate = 1.15,
                                   T_opt_invertebrate = 298.15, T_opt_vertebrate = 298.15,
                                   β_producer = 0.25, β_invertebrate = 0.25, β_vertebrate = 0.25); passed_temp_parameters...)
     k = 8.617e-5 # Boltzmann constant
	 if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

     return(bodymass, T, p) -> for i in 1:1
                                # parameters vary if the consumer is a vertebrate/invertebrate
                                norm_constant_all = temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                activation_energy_all = temperature_param.activation_energy_vertebrate .* p[:vertebrates] .+ temperature_param.activation_energy_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                deactivation_energy_all = temperature_param.deactivation_energy_vertebrate .* p[:vertebrates] .+ temperature_param.deactivation_energy_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                T_opt_all = temperature_param.T_opt_vertebrate .* p[:vertebrates] .+ temperature_param.T_opt_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                β_consumer = temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                Δenergy = deactivation_energy_all .- activation_energy_all
                                # β for resources
                                β_resource = zeros(Float64,(p[:S], p[:S]))
                                for consumer in 1:p[:S]
                                for resource in 1:p[:S]
                                  if p[:A][consumer, resource] == 1
                                    if p[:is_producer][resource]
                                        β_resource[consumer, resource] = temperature_param.β_producer
                                    elseif p[:vertebrates][resource]
                                        β_resource[consumer, resource] = temperature_param.β_vertebrate
                                    else
                                        β_resource[consumer, resource] = temperature_param.β_invertebrate
                                    end
                                 end
                                end
                                end
                                rate = norm_constant_all .* bodymass .^(β_consumer) .* bodymass' .^(β_resource) .* exp.(.-activation_energy_all ./ (k * T)) .* (1 ./ (1 .+ exp.(-1 / (k * T) .* (deactivation_energy_all .- (deactivation_energy_all ./ T_opt_all .+ k .* log.(activation_energy_all ./ Δenergy)).* T))))
                                rate[isnan.(rate)] .= 0
                                return  rate

                             end
end

"""
**Option 4 : Gaussian function for growth rate**

This function can be called with the keywords :growth, :growthrate or :r as an argument in `Gaussian`, itself called in `model_parameters`, to define a Gaussian function (Amarasekare et al 2015) for producers growth rate.
Example : model_parameters(A, growthrate = Gaussian(:r))

| Parameter    | Meaning                                        | Default values | Reference            |
|:-------------|:-----------------------------------------------|:---------------|:---------------------|
| norm_constant| minimal/maximal trait value                    | 1              | NA                   |
| range        | performance breath (width of function)         | 20             | Amarasekare 2015     |
| T_opt        | temperature at which trait value is maximal    | 298.15         | Amarasekare 2015     |
| β            | allometric exponent                            | -0.25          | Gillooly et al 2002  |

Default values are given as an example.

The function can be called with the default parameters:
	- gaussian_r()
Parameters can also be specified:
	- gaussian_r(passed_temp_parameters = (norm_constant = 1, range = 20, T_opt = 298.15, β = -0.25))
"""
function gaussian_r(default_temp_parameters = (norm_constant = 1, range = 20, T_opt = 298.15, β = -0.25); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end
    return(bodymass, T, p) -> bodymass .^ temperature_param.β .* temperature_param.norm_constant .* exp(.-(T .- temperature_param.T_opt) .^ 2 ./ (2 .*temperature_param.range .^ 2))
end

"""
**Option 4 : Gaussian function for metabolic rate**

This function can be called with the keywords :metabolism, :x or :metabolicrate as an argument in `Gaussian`, itself called in `model_parameters`, to define a Gaussian function (Amarasekare et al 2015) for metabolic rate.
Example : model_parameters(A, growthrate = Gaussian(:r))

| Parameter                  | Meaning                                                       | Default values | Reference            |
|:---------------------------|:--------------------------------------------------------------|:---------------|:---------------------|
| norm_constant_producer     | maximal trait value for producers    	                     | 0.2            | NA                   |
| norm_constant_invertebrate | maximal trait value for invertebrates  			             | 0.35           | NA                   |
| norm_constant_vertebrate   | maximal trait value for vertebrates     			             | 0.9            | NA                   |
| T_opt_producer             | temperature at which trait value is maximal for producers     | 298.15         | Amarasekare 2015     |
| T_opt_invertebrate         | temperature at which trait value is maximal for invertebrates | 298.15         | Amarasekare 2015     |
| T_opt_vertebrate           | temperature at which trait value is maximal for vertebrates   | 298.15         | Amarasekare 2015     |
| β_producer                 | allometric exponent for producers                             | -0.25          | Gillooly et al 2002  |
| β_invertebrate             | allometric exponent for invertebrates                         | -0.25          | Gillooly et al 2002  |
| β_vertebrate               | allometric exponent for vertebrates                           | -0.25          | Gillooly et al 2002  |
| range_producer             | performance breath (width of function) for producers          | 20             | Amarasekare 2015     |
| range_invertebrate         | performance breath (width of function) for invertebrates      | 20             | Amarasekare 2015     |
| range_vertebrate           | performance breath (width of function) for vertebrates        | 20             | Amarasekare 2015     |

Default values are given as an example.

The function can be called with the default parameters:
	- gaussian_x()
Parameters can also be specified:
	- gaussian_x(passed_temp_parameters = (norm_constant_producer = 0.2, norm_constant_invertebrate = 0.35, norm_constant_vertebrate = 0.9,
	                             range_producer = 20, range_invertebrate = 20, range_vertebrate = 20,
	                             T_opt_producer = 298.15, T_opt_invertebrate = 298.15, T_opt_vertebrate = 298.15,
	                             β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25))
"""
function gaussian_x(default_temp_parameters = (norm_constant_producer = 0.2, norm_constant_invertebrate = 0.35, norm_constant_vertebrate = 0.9,
                             range_producer = 20, range_invertebrate = 20, range_vertebrate = 20,
                             T_opt_producer = 298.15, T_opt_invertebrate = 298.15, T_opt_vertebrate = 298.15,
                             β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

     return(bodymass, T, p) -> for i in 1:1
                                 norm_constant_all = temperature_param.norm_constant_producer .* p[:is_producer] .+ temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 T_opt_all = temperature_param.T_opt_producer .* p[:is_producer] .+ temperature_param.T_opt_vertebrate .* p[:vertebrates] .+ temperature_param.T_opt_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 β_all = temperature_param.β_producer .* p[:is_producer] .+ temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 range_all = temperature_param.range_producer .* p[:is_producer] .+ temperature_param.range_vertebrate .* p[:vertebrates] .+ temperature_param.range_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                 return bodymass .^ β_all .* norm_constant_all .* exp.( .- (T .- T_opt_all) .^ 2 ./ (2 .* range_all .^ 2))
                               end
end

"""
**Option 4 : Gaussian function for attack rate**

This function can be called with the keywords :attackrate as an argument in `Gaussian`, itself called in `model_parameters`, to define a Gaussian function (Amarasekare et al 2015) for attack rate.
Example : model_parameters(A, attackrate = Gaussian(:attackrate))

| Parameter                  | Meaning                                                       | Default values | Reference            |
|:---------------------------|:--------------------------------------------------------------|:---------------|:---------------------|
| norm_constant_invertebrate | minimal/maximal trait value for invertebrates                 | 16             | NA                   |
| norm_constant_vertebrate   | minimal/maximal trait value for vertebrates                   | 16             | NA                   |
| range_invertebrate         | performance breath (width of function) for invertebrates      | 20             | Amarasekare 2015     |
| range_vertebrate           | performance breath (width of function) for vertebrates        | 20             | Amarasekare 2015     |
| T_opt_invertebrate         | temperature at which trait value is maximal                   | 298.15         | Amarasekare 2015     |
| T_opt_vertebrate           | temperature at which trait value is maximal for invertebrates | 298.15         | Amarasekare 2015     |
| β_producer                 | allometric exponent for producers                             | -0.25          | Gillooly et al 2002  |
| β_invertebrate             | allometric exponent for invertebrates                         | -0.25          | Gillooly et al 2002  |
| β_vertebrate               | allometric exponent for vertebrates                           | -0.25          | Gillooly et al 2002  |

The function can be called with the default parameters:
	- gaussian_attackr()
Parameters can also be specified:
	- gaussian_attackr(passed_temp_parameters = (norm_constant_invertebrate = 16, norm_constant_vertebrate = 16,
	                   range_invertebrate = 20, range_vertebrate = 20,
	                   T_opt_invertebrate = 298.15, T_opt_vertebrate = 298.15,
	                   β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25))
"""
function gaussian_attackr(default_temp_parameters = (norm_constant_invertebrate = 16, norm_constant_vertebrate = 16,
                                    range_invertebrate = 20, range_vertebrate = 20,
                                    T_opt_invertebrate = 295, T_opt_vertebrate = 295,
                                    β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

	return(bodymass, T, p) -> for i in 1:1
                                norm_constant_all = temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                T_opt_all = temperature_param.T_opt_vertebrate .* p[:vertebrates] .+ temperature_param.T_opt_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                β_consumer = temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                range_all = temperature_param.range_vertebrate .* p[:vertebrates] .+ temperature_param.range_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                # β for resources
                                β_resource = zeros(Float64,(p[:S], p[:S]))
                                for consumer in 1:p[:S]
                                for resource in 1:p[:S]
                                  if p[:A][consumer, resource] == 1
                                    if p[:is_producer][resource]
                                        β_resource[consumer, resource] = temperature_param.β_producer
                                    elseif p[:vertebrates][resource]
                                        β_resource[consumer, resource] = temperature_param.β_vertebrate
                                    else
                                        β_resource[consumer, resource] = temperature_param.β_invertebrate
                                    end
                                 end
                                end
                                end
                                rate = bodymass .^ β_consumer .* bodymass' .^ β_resource .* norm_constant_all .* exp.( .- (T .- T_opt_all) .^ 2 ./ (2 .* range_all .^ 2))
                                rate[isnan.(rate)] .= 0
                                return rate
                            end
end

"""
**Option 4 : Gaussian function for handling time**

This function can be called with the keyword :handlingtime as an argument in `Gaussian`, itself called in `model_parameters`, to define a Gaussian function (Amarasekare et al 2015) for handling time.
Example : model_parameters(A, handlingtime = Gaussian(:handlingtime))

| Parameter                  | Meaning                                                       | Default values | Reference            |
|:---------------------------|:--------------------------------------------------------------|:---------------|:---------------------|
| norm_constant_invertebrate | minimal/maximal trait value for invertebrates                 | 0.125          | NA                   |
| norm_constant_vertebrate   | minimal/maximal trait value for vertebrates                   | 0.125          | NA                   |
| range_invertebrate         | performance breath (width of function) for invertebrates      | 20             | Amarasekare 2015     |
| range_vertebrate           | performance breath (width of function) for vertebrates        | 20             | Amarasekare 2015     |
| T_opt_invertebrate         | temperature at which trait value is maximal                   | 298.15         | NA   				 |
| T_opt_vertebrate           | temperature at which trait value is maximal for invertebrates | 298.15         | NA					 |
| β_producer                 | allometric exponent for producers                             | -0.25          | Gillooly et al 2002  |
| β_invertebrate             | allometric exponent for invertebrates                         | -0.25          | Gillooly et al 2002  |
| β_vertebrate               | allometric exponent for vertebrates                           | -0.25          | Gillooly et al 2002  |

The function can be called with the default parameters:
	- gaussian_handlingt()
Parameters can also be specified:
	- gaussian_handlingt(passed_temp_parameters = (norm_constant_invertebrate = 0.5, norm_constant_vertebrate = 0.5,
							range_invertebrate = 20, range_vertebrate = 20,
							T_opt_invertebrate = 295, T_opt_vertebrate = 295,
							β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25))
"""
function gaussian_handlingt(default_temp_parameters = (norm_constant_invertebrate = 0.5, norm_constant_vertebrate = 0.5,
                                    range_invertebrate = 20, range_vertebrate = 20,
                                    T_opt_invertebrate = 295, T_opt_vertebrate = 295,
                                    β_producer = -0.25, β_invertebrate = -0.25, β_vertebrate = -0.25); passed_temp_parameters...)
	if length(passed_temp_parameters) != 0
	  tmpargs = passed_temp_parameters[:passed_temp_parameters]
	  temperature_param = merge(default_temp_parameters, tmpargs)
	else
	  temperature_param = default_temp_parameters
	end

	return(bodymass, T, p) -> for i in 1:1
                                norm_constant_all = temperature_param.norm_constant_vertebrate .* p[:vertebrates] .+ temperature_param.norm_constant_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                T_opt_all = temperature_param.T_opt_vertebrate .* p[:vertebrates] .+ temperature_param.T_opt_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                β_consumer = temperature_param.β_vertebrate .* p[:vertebrates] .+ temperature_param.β_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                range_all = temperature_param.range_vertebrate .* p[:vertebrates] .+ temperature_param.range_invertebrate .* (.!p[:vertebrates] .& .!p[:is_producer])
                                # β for resources
                                β_resource = zeros(Float64,(p[:S], p[:S]))
                                for consumer in 1:p[:S]
                                for resource in 1:p[:S]
                                  if p[:A][consumer, resource] == 1
                                    if p[:is_producer][resource]
                                        β_resource[consumer, resource] = temperature_param.β_producer
                                    elseif p[:vertebrates][resource]
                                        β_resource[consumer, resource] = temperature_param.β_vertebrate
                                    else
                                        β_resource[consumer, resource] = temperature_param.β_invertebrate
                                    end
                                 end
                                end
                                end
                                rate = bodymass .^ β_consumer .* bodymass' .^ β_resource .* norm_constant_all .* exp.((T .- T_opt_all) .^ 2 ./ (2 .* range_all .^ 2))
                                rate[isnan.(rate)] .= 0
                                return rate
                            end
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
export Stagewise

function Stagewise(prob::Problem{L, B, CS},
                   ϵ::Real,
                   subsolver::AbstractNextBasisSolver=defaultsolver(prob);
                   ad::Bool=false,
                   lasso::Bool=false,
                   tol::Real=1e-8,
                   maxiter::Int=1000) where {L, F, B<:AbBasis{F}, CS}
    ens = SWEnsemble(prob, ϵ, ϵ)
    enses = BBEnsemble{F,B}[]
    iter = 0
    λ = Inf
    currloss = trainerr(ens)
    # llog(txt="") = printfmtln("{:3d}) {:.3f} λ={:.3f} |β|={:.1f}ϵ({}) {}", iter, currloss, λ, sum(abs, ens.β)/ϵ, countnz(ens.β), txt)
    llog(t="") = nothing
    llog()

    # Get first ω
    objtv, grad = lossngrad(prob, ens.fit)
    ω, ρ, score = getnextω(prob, subsolver, grad)
    pushnew!(ens, ω)
    ad && push!(enses, BBEnsemble(ens))
    newloss = trainerr(ens)
    λ = (currloss-newloss)/ϵ
    currloss = newloss
    iter += 1
    llog()

    ϵlosses(d) = [(βi≈0. ? Inf : loss(prob, ens.fit+d*ϵ*z)) for (βi,z) in zip(ens.β, ens.zs)]

    for iter in 2:maxiter
        prevloss = currloss

        # Try decreasing |β|
        dec_losses = ϵlosses(-1)
        declossstar, distar = findmin(dec_losses)
        if declossstar - (lasso ? (λ*ϵ) : tol) ≤ currloss - tol       # W/o tol, guy just in will be out
            decrementβi!(ens, distar)
            currloss = declossstar
            llog("⇩ $distar ! $(ens.β[distar])")
        else

        # Get a new predictor
          objtv, grad = lossngrad(prob, ens.fit)
          ω, ρ, score = getnextω(prob, subsolver, grad)
          z = predictor(prob, ω)
          newzloss = loss(prob, ens.fit + ϵ*z)

        # Try incrementing old ones
          inc_losses = ϵlosses(1)
          inclossstar, instar = findmin(inc_losses)

          if currloss < min(inclossstar, newzloss)
              printfmtln("Can not continue at {:d} c={:.4f} < (d={:.4f}({}), i={:.4f}({}), n={:.4f}).\n Try decreasing ϵ!",
                iter, currloss, declossstar, distar, inclossstar, instar, newzloss)
              break
          end

          if inclossstar ≤ newzloss
            augmentβi!(ens, instar)
            currloss = inclossstar
            # llog("⇑ $instar !")
          else
            pushnew!(ens, ω, z)
            currloss = newzloss
            llog("new (ω, β)!")
          end
          λ = min(λ, (prevloss-currloss-tol)/ϵ)
        end

        if ad
            ens, currloss = backpropωs(prob, ens)
            push!(enses, BBEnsemble(ens))
            llog("ad {ω}")
        end
    end

    llog("\nEND")
    ad ? enses : ens
end

ext: jl


text:
# Time to wait for the Docker containers to launch and the workers to connect to the
# manager
const DOCKER_TIMEOUT = Minute(1)

"""
    DockerManager(num_workers; kwargs...)

A cluster manager which spawns workers via a locally running [Docker](https://docs.docker.com/)
daemon service. Typically used on a single machine to debug multi-machine Julia code.

In order to make a local Docker cluster you'll need to have an available Docker image that
has Julia, a version of AWSClusterManagers which includes DockerManager, and the docker
cli all baked into the image.

You can then create a Docker container which is capable of spawning additional Docker
containers via:

   docker run --network=host -v /var/run/docker.sock:/var/run/docker.sock --rm -it <image> julia

**Note**: The host networking is required for the containers to be able to intercommunicate.
The Docker host's UNIX socket needs to be forwarded so we can allow the container to
communicate with host Docker process.

## Arguments
- `num_workers::Int`: The number of workers to spawn

## Keywords
- `image::AbstractString`: The docker image to run.
- `timeout::Second`: The maximum number of seconds to wait for workers to become available
  before aborting.

## Examples
```julia
julia> addprocs(DockerManager(4, "myproject:latest"))
```
"""
DockerManager

struct DockerManager <: ContainerManager
    num_workers::Int
    image::AbstractString
    timeout::Second
    min_ip::IPv4
    max_ip::IPv4

    function DockerManager(
        num_workers::Integer,
        image::AbstractString,
        timeout::Period=DOCKER_TIMEOUT,
        min_ip::IPv4=ip"0.0.0.0",
        max_ip::IPv4=ip"255.255.255.255",
    )
        num_workers >= 0 || throw(ArgumentError("num workers must be non-negative"))

        # Workers by default inherit the defaults from the manager.
        if isempty(image)
            image = @mock image_id()
        end

        return new(num_workers, image, Second(timeout), min_ip, max_ip)
    end
end

function DockerManager(
    num_workers::Integer;
    image::AbstractString="",
    timeout::Union{Real,Period}=DOCKER_TIMEOUT,
    min_ip::IPv4=ip"0.0.0.0",
    max_ip::IPv4=ip"255.255.255.255",
)
    return DockerManager(num_workers, image, timeout, min_ip, max_ip)
end

launch_timeout(mgr::DockerManager) = mgr.timeout
desired_workers(mgr::DockerManager) = mgr.num_workers, mgr.num_workers

function Base.:(==)(a::DockerManager, b::DockerManager)
    return (a.num_workers == b.num_workers && a.image == b.image && a.timeout == b.timeout)
end

function spawn_containers(mgr::DockerManager, override_cmd::Cmd)
    # Requires that the `docker` is installed
    cmd = `docker run --detach --network=host --rm $(mgr.image)`
    cmd = `$cmd $override_cmd`

    # Docker only allow us to spawn a job at a time
    for id in 1:(mgr.num_workers)
        container_id = @mock read(cmd, String)
        notice(LOGGER, "Spawning container: $container_id")
    end
end

# Determine the image ID of the currently running container
function image_id(container_id::AbstractString=container_id())
    json = JSON.parse(read(`docker container inspect $container_id`, String))
    return last(split(json[1]["Image"], ':'))
end


text:
"""
Generic plugins are plugins that add any number of patterns to the generated package's
`.gitignore`, and have at most one associated file to generate.

# Attributes
* `gitignore::Vector{AbstractString}`: Array of patterns to be added to the `.gitignore` of
  generated packages that use this plugin.
* `src::Union{AbstractString, Nothing}`: Path to the file that will be copied into the generated
  package repository. If set to `nothing`, no file will be generated. When this defaults
  to an empty string, there should be a default file in `defaults` that will be copied.
  That file's name is usually the same as the plugin's name, except in all lowercase and
  with the `.yml` extension. If this is not the case, an `interactive` method needs to be
  implemented to call `interactive(; file="file.ext")`.
* `dest::AbstractString`: Path to the generated file, relative to the root of the generated
  package repository.
* `badges::Vector{Badge}`: Array of [`Badge`](@ref)s containing information used to
  create Markdown-formatted badges from the plugin. Entries will be run through
  [`substitute`](@ref), so they may contain placeholder values.
* `view::Dict{String, Any}`: Additional substitutions to make in both the plugin's badges
  and its associated file. See [`substitute`](@ref) for details.

# Example
```julia
struct MyPlugin <: GenericPlugin
    gitignore::Vector{AbstractString}
    src::Union{AbstractString, Nothing}
    dest::AbstractString
    badges::Vector{Badge}
    view::Dict{String, Any}

    function MyPlugin(; config_file::Union{AbstractString, Nothing}="")
        if config_file != nothing
            config_file = if isempty(config_file)
                joinpath(DEFAULTS_DIR, "my-plugin.toml")
            elseif isfile(config_file)
                abspath(config_file)
            else
                throw(ArgumentError(
                    "File \$(abspath(config_file)) does not exist"
                ))
            end
        end
        new(
            ["*.mgp"],
            config_file,
            ".my-plugin.toml",
            [
                Badge(
                    "My Plugin",
                    "https://myplugin.com/badge-{{YEAR}}.png",
                    "https://myplugin.com/{{USER}}/{{PKGNAME}}.jl",
                ),
            ],
            Dict{String, Any}("YEAR" => year(today())),
        )
    end
end

interactive(::Type{MyPlugin}) = interactive(MyPlugin; file="my-plugin.toml")
```

The above plugin ignores files ending with `.mgp`, copies `defaults/my-plugin.toml` by
default, and creates a badge that links to the project on its own site, using the default
substitutions with one addition: `{{YEAR}} => year(today())`. Since the default config
template file doesn't follow the generic naming convention, we added another `interactive`
method to correct the assumed filename.
"""
abstract type GenericPlugin <: Plugin end

function Base.show(io::IO, p::GenericPlugin)
    spc = "  "
    println(io, nameof(typeof(p)), ":")

    cfg = if p.src === nothing
        "None"
    else
        dirname(p.src) == DEFAULTS_DIR ? "Default" : p.src
    end
    println(io, spc, "→ Config file: ", cfg)

    n = length(p.gitignore)
    s = n == 1 ? "" : "s"
    print(io, spc, "→ $n gitignore entrie$s")
    n > 0 && print(io, ": ", join(map(repr, p.gitignore), ", "))
end

"""
Custom plugins are plugins whose behaviour does not follow the [`GenericPlugin`](@ref)
pattern. They can implement [`gen_plugin`](@ref), [`badges`](@ref), and
[`interactive`](@ref) in any way they choose, as long as they conform to the usual type
signature.

# Attributes
* `gitignore::Vector{AbstractString}`: Array of patterns to be added to the `.gitignore` of
  generated packages that use this plugin.

# Example
```julia
struct MyPlugin <: CustomPlugin
    gitignore::Vector{AbstractString}
    lucky::Bool

    MyPlugin() = new([], rand() > 0.8)

    function gen_plugin(p::MyPlugin, t::Template, pkg_name::AbstractString)
        return if p.lucky
            text = substitute("You got lucky with {{PKGNAME}}, {{USER}}!", t)
            gen_file(joinpath(t.dir, pkg_name, ".myplugin.yml"), text)
            [".myplugin.yml"]
        else
            println("Maybe next time.")
            String[]
        end
    end

    function badges(p::MyPlugin, user::AbstractString, pkg_name::AbstractString)
        return if p.lucky
            [
                format(Badge(
                    "You got lucky!",
                    "https://myplugin.com/badge.png",
                    "https://myplugin.com/\$user/\$pkg_name.jl",
                )),
            ]
        else
            String[]
        end
    end
end

interactive(:Type{MyPlugin}) = MyPlugin()
```

This plugin doesn't do much, but it demonstrates how [`gen_plugin`](@ref), [`badges`](@ref)
and [`interactive`](@ref) can be implemented using [`substitute`](@ref),
[`gen_file`](@ref), [`Badge`](@ref), and [`format`](@ref).

# Defining Template Files
Often, the contents of the config file that your plugin generates depends on variables like
the package name, the user's username, etc. Template files (which are stored in `defaults`)
can use [here](https://github.com/jverzani/Mustache.jl)'s syntax to define replacements.
"""
abstract type CustomPlugin <: Plugin end

"""
    Badge(hover::AbstractString, image::AbstractString, link::AbstractString) -> Badge

A `Badge` contains the data necessary to generate a Markdown badge.

# Arguments
* `hover::AbstractString`: Text to appear when the mouse is hovered over the badge.
* `image::AbstractString`: URL to the image to display.
* `link::AbstractString`: URL to go to upon clicking the badge.
"""
struct Badge
    hover::String
    image::String
    link::String
end

"""
    format(b::Badge) -> String

Return `badge`'s data formatted as a Markdown string.
"""
format(b::Badge) = "[![$(b.hover)]($(b.image))]($(b.link))"

"""
    gen_plugin(p::Plugin, t::Template, pkg_name::AbstractString) -> Vector{String}

Generate any files associated with a plugin.

# Arguments
* `p::Plugin`: Plugin whose files are being generated.
* `t::Template`: Template configuration.
* `pkg_name::AbstractString`: Name of the package.

Returns an array of generated file/directory names.
"""
gen_plugin(::Plugin, ::Template, ::AbstractString) = String[]

function gen_plugin(p::GenericPlugin, t::Template, pkg_name::AbstractString)
    if p.src === nothing
        return String[]
    end
    text = substitute(
        read(p.src, String),
        t;
        view=merge(Dict("PKGNAME" => pkg_name), p.view),
    )
    gen_file(joinpath(t.dir, pkg_name, p.dest), text)
    return [p.dest]
end

"""
    badges(p::Plugin, user::AbstractString, pkg_name::AbstractString) -> Vector{String}

Generate Markdown badges for the plugin.

# Arguments
* `p::Plugin`: Plugin whose badges we are generating.
* `user::AbstractString`: Username of the package creator.
* `pkg_name::AbstractString`: Name of the package.

Returns an array of Markdown badges.
"""
badges(::Plugin, ::AbstractString, ::AbstractString) = String[]

function badges(p::GenericPlugin, user::AbstractString, pkg_name::AbstractString)
    # Give higher priority to replacements defined in the plugin's view.
    view = merge(Dict("USER" => user, "PKGNAME" => pkg_name), p.view)
    return map(b -> substitute(format(b), view), p.badges)
end

"""
    interactive(T::Type{<:Plugin}; file::Union{AbstractString, Nothing}="") -> Plugin

Interactively create a plugin of type `T`, where `file` is the plugin type's default
config template with a non-standard name (for `MyPlugin`, this is anything but
"myplugin.yml").
"""
function interactive(T::Type{<:GenericPlugin}; file::Union{AbstractString, Nothing}="")
    name = string(nameof(T))
    # By default, we expect the default plugin file template for a plugin called
    # "MyPlugin" to be called "myplugin.yml".
    fn = file != nothing && isempty(file) ? "$(lowercase(name)).yml" : file
    default_config_file = fn == nothing ? fn : joinpath(DEFAULTS_DIR, fn)

    print("$name: Enter the config template filename (\"None\" for no file) ")
    if default_config_file == nothing
        print("[None]: ")
    else
        print("[", replace(default_config_file, homedir() => "~"), "]: ")
    end

    config_file = readline()
    config_file = if uppercase(config_file) == "NONE"
        nothing
    elseif isempty(config_file)
        default_config_file
    else
        config_file
    end

    return T(; config_file=config_file)
end


text:
module JuliaPackageTutorial

using ForwardDiff

include("extra.jl")

end # module

max_stars_repo_name: HazilMohamed/JuliaPackageTutorial

max_issues_repo_name: HazilMohamed/JuliaPackageTutorial

max_forks_repo_name: HazilMohamed/JuliaPackageTutorial


text:
# Julia wrapper for header: /usr/include/libavcodec/ac3_parser.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_ac3_parse_header

function av_ac3_parse_header(buf, size::Csize_t, bitstream_id, frame_size)
    ccall((:av_ac3_parse_header, libavcodec), Cint, (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Ptr{UInt16}), buf, size, bitstream_id, frame_size)
end


text:
usingstatic =   isdefined(TensorCast, :StaticArray)

@testset "scope" begin

    bc = rand(2,3)

    β = 99 # should be ignored
    cb = nothing # should be overwritten
    @cast cb[κ,β] := bc[β,κ]
    @test !isdefined(@__MODULE__,:κ)
    @test β == 99
    @test all(cb .== transpose(bc))

    @cast cb[c,b] = bc[b,c] c:3 # in-place

    β = 2
    @cast f[(c,b)] := bc[b,c] b:β, c:3 # @assert uses β
    @test size(f) == (6,)
    @cast cb2[c,b] := f[(c,b)] b:β # can't not use this
    @test size(cb2) == (3,2)
    @test all(cb2 .== transpose(bc))

    oo = [1,1]
    @cast cb3[c,b] := f[(c,b)] b:sum(oo) # with an expression
    @test size(cb3) == (3,2)

    β = 2
    @cast b2[b] := bc[b,$β] # β should be inserted
    @test b2 == bc[:,2]

end
@testset "permute" begin

    bcde = rand(2,3,4,5);
    bcdef = rand(2,3,4,5,6);

    @cast oo[d,e,b,c] |= bcde[b,c,d,e];
    @test size(oo) == (4,5,2,3)
    @test oo[1,3,1,3] == bcde[1,3,1,3]

    oo[1] = 99
    @test bcde[1] != 99 # copy not a view

    @cast oo[d,e,b,c] = bcde[b,c,d,e]  assert;
    oo[1] = 99
    @test bcde[1] != 99 # still not a view

    @cast oo[e,b,c,f,d] := bcdef[b,c,d,e,f];
    @test size(oo) == (5,2,3,6,4)
    @test oo[1,2,1,2,2] == bcdef[2,1,2,1,2]

    @cast oo[e,b,zc,f,d] := bcdef[b,zc,d,e,f]; # because I alphabetise symbols
    @test size(oo) == (5,2,3,6,4)

    @cast oo[e,zb,c,f,d] := bcdef[zb,c,d,e,f];
    @test size(oo) == (5,2,3,6,4)

end
@testset "slice" begin

    bcd = rand(2,3,4);
    bcde = rand(2,3,4,5);

    @cast Bcd[b][c,d] |= bcd[b,c,d]  assert
    @test size(Bcd) == (2,)
    @test size(first(Bcd)) == (3,4)
    @test Bcd[2][3,1] == bcd[2,3,1]

    Bcd[2][3,1] = 99
    @test bcd[2,3,1] != 99 # copy not a view

    @cast Cbd[c][b,d] := bcd[b,c,d]  assert
    @test size(Cbd) == (3,)
    @test size(first(Cbd)) == (2,4)
    @test Cbd[3][2,1] == bcd[2,3,1]

    @cast BCde[b,c][d,e] := bcde[b,c,d,e]  assert
    @test size(BCde) == (2,3)
    # @test size(first(BCde)) == (4,5)
    @test size(BCde[1,1]) == (4,5) # new julienne

    @cast DBec[d,b][e,c] := bcde[b,c,d,e]  assert
    @test size(DBec) == (4,2)
    # @test size(first(DBec)) == (5,3)
    @test size(DBec[1,1]) == (5,3) # new julienne

    @cast DEbc[d,e][b,c] := bcde[b,c,d,e] # this order can be a view
    @test size(DEbc) == (4,5)
    # @test size(first(DEbc)) == (2,3)
    @test size(DEbc[1,1]) == (2,3)

    DEbc[3,4][1,2] = 99
    @test bcde[1,2,3,4] == 99 # view not a copy

end
@testset "static-slice" begin

    M = rand(Int, 2,3)
    S1 = reinterpret(SVector{2,Int}, vec(M)) # as in readme
    @cast S2[k]{i} := M[i,k]  i:2
    M[end]=42
    @test S1[3][2]==42 # is a view
    @test S2[3][2]==42
    @test first(S1) isa StaticArray
    @test first(S2) isa StaticArray
    @test size(first(S1)) == (2,)
    @test size(first(S2)) == (2,)

    @cast S3[k]{i} := M[k,i]  i:3 # this has a permutedims
    @test first(S3) isa StaticArray
    @test size(first(S3)) == (3,)

    bcd = rand(2,3,4);
    bcde = rand(2,3,4,5);

    @cast DBec[d,b]{e,c} := bcde[b,c,d,e] e:5, c:3
    @test size(DBec) == (4,2)
    @test size(first(DBec)) == (5,3)
    @test first(DBec) isa StaticArray

    @cast DEbc[d,e]{b,c} := bcde[b,c,d,e] b:2, c:3
    @test size(DEbc) == (4,5)
    @test size(first(DEbc)) == (2,3)
    @test first(DEbc) isa StaticArray

    A = rand(3,3,4)
    @cast B[k]{i,j} := A[i,j,k]  i:3, j:3
    @cast C[(j,k),i] := B[k][i,j]
    @test C[1,2] == B[1][2,1] ==  A[2,1,1]

end
@testset "glue" begin

    bcd = rand(2,3,4)
    Bcd = [ rand(3,4) for b=1:2 ]

    @cast bcd[b,c,d] = Bcd[b][c,d]  assert
    @test size(bcd) == (2,3,4)
    @test Bcd[2][3,1] == bcd[2,3,1]

    Cbd = [ rand(2,4) for c=1:3 ]

    @cast oo[b,c,d] := Cbd[c][b,d]
    @test size(oo) == (2,3,4)
    @test Cbd[3][2,1] == oo[2,3,1]

    DEbc = [ rand(2,3) for d=1:4, e=1:5 ]

    @cast bcde[b,c,d,e] := DEbc[d,e][b,c];
    @test size(bcde) == (2,3,4,5)

    BCde = [ rand(4,5) for d=1:2, e=1:3 ]

    @cast bcde[b,c,d,e] = BCde[b,c][d,e]  assert;
    @test size(bcde) == (2,3,4,5)

end
@testset "static-glue" begin

    B = [ SVector{3}(rand(3)) for i=1:2]

    @cast A[j,i] := B[i]{j} i:2, j:3 # view is impossible without Static slices
    @test size(A) == (3,2)
    @test !isa(A, Array)

    A[1,2] = 33
    @test B[2][1] == 33 # mutating B illegally?

    @cast A[i,j] := B[i]{j} j:3
    @test size(A) == (2,3)
    @test !isa(A, StaticArray)

    @test A[2,1] == 33

    C = [ SMatrix{2,3}(rand(2,3)) for i=1:4, j=1:5 ];

    @cast A[l,k,j,i] := C[i,j]{k,l} k:2, l:3
    @test A[1,2,3,4] == C[4,3][2,1]
    @test size(A) == (3,2,5,4)

    @cast A[l,k,j,i] = C[i,j]{k,l} k:2, l:3 # now in-place
    @test A[1,2,3,4] == C[4,3][2,1]

    @cast A[l,k,j,i] = C[i,j]{k,l} # now without static glue
    @test A[1,2,3,4] == C[4,3][2,1]

end
@testset "reshape" begin

    bcd = rand(2,3,4)
    bcde = rand(2,3,4,5)
    bcde2 = similar(bcde)

    @cast v[(b,zc,d)] := bcd[b,zc,d]
    @cast w[d,(zc,b)] := v[(b,zc,d)] b:2, zc:3 # add z sometimes to mess with alphabetisation
    @cast dbc[d,b,c] := w[d,(c,b)] c:3

    @test size(w) == (4,6)
    @test size(dbc) == (4,2,3)
    @test dbc[1,2,3] == bcd[2,3,1]

    @cast g[(b,c),x,y,e] := bcde[b,c,(x,y),e] x:2;
    @test size(g) == (6,2,2,5)

    @cast g[(y,e),(b,c),x] := bcde[b,c,(x,y),e] x:1;
    @test size(g) == (20, 6, 1)
    @cast bcde2[b,c,(x,y),e] = g[(y,e),(b,c),x];
    @test all(bcde2 .== bcde)

end
@testset "combined" begin

    Bc = [ rand(3) for b=1:2 ] # Bc = [(1:3) .+ 10i for i=1:2]
    @cast f[(b,c)] := Bc[b][c]
    @test size(f) == (6,)
    @cast bc[b,c] |= f[(b,c)] b:2 # for in-place to work below, this bc is NOT a view of f
    @test size(bc) == (2,3)
    @test bc[1,2] == Bc[1][2]

    @cast f[(c,b)] = Bc[b][c]  assert  # in-place
    @cast bc[b,c] = f[(c,b)]   assert  # in-place
    @test bc[1,2] == Bc[1][2]

    Cdaeb = [rand(4,1,5,2) for i=1:3]; # sizes match abcde now
    Z = @cast [(c,d),e][a,b] := Cdaeb[c][d,a,e,b] # alphabetical = canonical
    @test size(Z) == (3*4, 5)
    @test size(first(Z)) == (1,2)
    @test Z[2*1,4][1,2] == Cdaeb[2][1,1,4,2]

    @cast Z[e,(c,d)][b,a] := Cdaeb[c][d,a,e,b] # same but with permutedims
    @test size(Z) == (5, 3*4)
    @test size(first(Z)) == (2,1)
    @test Z[4,2*1][2,1] == Cdaeb[2][1,1,4,2]

end
@testset "reverse" begin

    bc = rand(2,3)

    @cast B[i,j] := bc[i,-j]
    D = reverse(bc, dims=2)
    @test B == D

    ccbb = rand(0:99, 3,3,2,2)
    @cast A[b,c,d,e] := ccbb[b,c,-d,e]
    B = similar(A);
    @cast B[b,c,d,e] = ccbb[b,c,-d,e]
    C = reverse(ccbb, dims=3)
    @test A == B == C

    ccbb = rand(0:99, 3,3,2,2)
    @cast A[c,e,b,d] := ccbb[-b,c,d,e] + 100
    @cast B[c,e,b,d] := ccbb[b,c,d,e] + 100
    @test A == reverse(B, dims=3)
    C = similar(A); D = similar(A);
    @cast C[c,e,b,d] = ccbb[-b,c,d,e] + 100
    @test A == C

    @cast A[c,e,b,_] := ccbb[-b,c,2,e] + 100
    @cast B[c,e,b,_] := ccbb[b,c,2,e] + 100
    @test A == reverse(B, dims=3)

end
@testset "fixed indices" begin

    W = rand(2,3,5,7);

    @cast Z[_,i,_,k] := W[2,k,4,i]

    @test size(Z) == (1,7,1,3)

    @cast A[k\i] := W[2,k,4,i] / Z[1,i,_,k]
    @test all(A .== 1)

end
@testset "random" begin

    A = zeros(1,2,3);
    B = randn(1,3*2)

    @cast A[i,j,k] = B[i,(k,j)] # errored on push_checks for a while

    @cast A[i,j,k] = B[i,(j,k)] i:1

    @cast V[(i,j,k)] := B[i,(k,j)]  k:3, assert # was an error for recursive colon reasons

    A = @cast [jk,i] := B[i,jk] # without left-hand name
    @test size(A) == (6,1)

    A = @cast [k,j,i] := B[i,(j,k)] k:3 # without left-hand name
    @test size(A) == (3,2,1)

    C = @cast [k][i,j] := B[i,(j,k)]  k:3 # without left-hand name
    C = @cast [k][i,j] := B[i,(j,k)]  j:2, k:3

    @test size(C) == (3,)
    @test size(first(C)) == (1,2)

    @cast A[k][i,j] := B[i,(j,k)]  k:length(C) # with an expression

    @test size(A) == (3,)

    D = randn(1*2*3, 4)
    # these two both result in a product of sizes like ((sz[2] * :) * s[4]) with : on 2nd level.
    @cast C[i,(j,k,l)] := D[(i,j,k),l] i:1, j:2 # gets sz... from colonise! too
    @cast E[i,(k,j,l)] := D[(i,j,k),l] i:1, j:2 # does not

    @test size(C) == (1, 2*3*4)
    @test size(E) == (1, 2*3*4)

    B = [rand(2) for i=1:3 ]
    A2 = @cast A[(i,j)] := B[i][j]
    A3 = @cast A[(i,j)] = B[i][j] # for a while this in-place thing returned the wrong shape
    @test size(A) == size(A2) == size(A3) == (6,)

    # This was broken above for a while:
    Bc = [(1:3) .+ 10i for i=1:2]
    @cast f[(b,c)] := Bc[b][c]
    @cast bc[b,c] := f[(b,c)] b:2 # this bc is a view of f, NB
    @cast f[(c,b)] = Bc[b][c]
    @cast bc[b,c] = f[(c,b)]  # thus this copyto! is confusing
    @test_broken bc[1,2] != Bc[1][2] # and in fact leads to duplicates here

end
@testset "errors" begin

    bc = rand(2,3)
    cb = rand(3,2)

    ## macro

    # @test_throws ArgumentError @cast a[i] := bc[i,j]
    # @test_throws ArgumentError @cast a[i,k] := bc[(i,j),k]
    # @test_throws ArgumentError @cast a[(i,j),k] := bc[(i,k)]

    ## runtime

    @test_throws DimensionMismatch @cast oo[i,j] := bc[i,j] i:3, assert
    @test_throws DimensionMismatch @cast cb[i,j] = bc[i,j]  assert

end


text:
module Resource
using ..Clarus
import CSV
export read,write

function openFile(filename, mode = "r")
  if isfile(filename)
    return open(filename,mode)
  end
  if isdir(Clarus.Services.credentials.resource_path)
    filename = joinpath(Clarus.Services.credentials.resource_path,filename)
  end
  return open(filename,mode)
end

function read(filenames)
  Clarus.Services.__init_credentials!__(Clarus.Services.credentials_init);
  return readfiles(filenames)
end

function readfiles(filenames)
  strings = String[]
  for filename in split(filenames,","; keep = false)
      f = openFile(strip(filename))
      if f != nothing
        push!(strings,readstring(f))
        close(f)
      end
  end
  return strings
end

function write(filename,r::Clarus.Services.Response)
  if length(dirname(filename)) == 0
    filename = joinpath(Clarus.Services.credentials.savefile_path,filename)
  end
  f = openFile(filename;mode="w")
  Base.write(f,r.httpresponse.data)
  close(f)
  return
end

end #MODULE-END


text:
### NOTE: this script is an addition to KissABC to implement
### multithreading support with in-place operations for the ABCDE method;
### for this we use FLoops and its @floop and @init macros

# TODO/NOTE: how to really check if multithreading runs correctly?
# maybe fix all random seeds

# IMPORTANT NOTE: in this code we use multithreading in the form of "filling pre-allocated
# output" (https://juliafolds.github.io/data-parallelism/tutorials/mutations/#filling_outputs)
# this can be unsafe and cause data races! (dict, sparsearrays, Bit vectors, views)
# but Arrays should be fine (and θs, logπ, Δs, nθs, nlogπ, nΔs are arrays)
# println(typeof(θs)) => Vector{KissABC.Particle{Tuple{Float64, Float64, Float64}}}
# println(typeof(logπ)) => Vector{Float64}
# println(typeof(Δs)) => Vector{Float64}

# NOTE: the op tuple operations (also on a Particle) seem to
# be zero-allocating (both immutable types!, no heap memory needed) and can be
# broadcasted (over multiple elements in a tuple or Particle.x); values inside
# the tuples are all plain data (int, float), so that this should be no problem
# for multithreading (although zero-allocation, nothing is changed in-place,
# just plain data calculation in stack memory); in the end written to nθs
# (which in each generation is created as deepcopy and creates allocations, by identity.());
# see also seems tests in polylox_env_ABC_HSC_test.jl
# so all things at θs[i], logπ[i], Δs[i] are immutable, such as Vector{SomeType}()
# with SomeType immutable as in my question (https://discourse.julialang.org/t/how-to-implement-multi-threading-with-external-in-place-mutable-variables/62610/3)

# NOTE: @reduce needed? if yes
# use solution like this https://discourse.julialang.org/t/using-floops-jl-to-update-array-counters/58805
# or this (histogram)? https://juliafolds.github.io/data-parallelism/tutorials/quick-introduction/

function abcde_init!(prior, dist!, varexternal, θs, logπ, Δs, nparticles, rng, ex)
    # calculate cost/dist for each particle
    # (re-draw parameters if not finite)

    @floop ex for i = 1:nparticles
        # @init allows re-using mutable temporary objects within each base case/thread
        # TODO/NOTE: in ThreadedEx mode one might observe a high % of garbage collection
        # it was hard to check if this @init really works on my varexternal object...
        # (compared with ve = deepcopy(varexternal) alone allocations were similar,
        # but also in sequential mode, suggesting it is just not driving the allocations)
        # NOTE: varexternal is the tuple wrapper of all mutatable external
        # variables as cellstate, stats, ... (same in abcde_swarm!)
        # NOTE: θs, logπ, Δs are read out but never written to, (for this
        # nθs, nlogπ, nΔs are used), so this should be data race free
        @init ve = deepcopy(varexternal)
        trng=rng

        # NOTE: I checked that Threads.threadid() also works in an floop
        # NOTE: seems to have something with thread-safe random numbers, see
        # https://discourse.julialang.org/t/multithreading-and-random-number-generators/49777/8
        # NOTE/TODO: maybe this can be improved performance-wise (see FLoops docs
        # on random numbers)
        ex!=SequentialEx() && (trng=Random.default_rng(Threads.threadid());)

        if isfinite(logπ[i])
            Δs[i] = dist!(θs[i].x, ve)
        end
        while (!isfinite(Δs[i])) || (!isfinite(logπ[i]))
            θs[i] = op(float, Particle(rand(trng, prior)))
            logπ[i] = logpdf(prior, push_p(prior,θs[i].x))
            Δs[i] = dist!(θs[i].x, ve)
        end
    end
end

function abcde_swarm!(prior, dist!, varexternal, θs, logπ, Δs, nθs, nlogπ, nΔs,
                    ϵ_pop, ϵ_target, γ, nparticles, nsims, earlystop, rng, ex)
    @floop ex for i in 1:nparticles
        @init ve = deepcopy(varexternal)

        if earlystop
            Δs[i] <= ϵ_target && continue
        end
        trng=rng
        ex!=SequentialEx() && (trng=Random.default_rng(Threads.threadid());)
        s = i
        ϵ = ifelse(Δs[i] <= ϵ_target, ϵ_target, ϵ_pop)
        if Δs[i] > ϵ
            # NOTE: Δs .<= Δs[i] is this data race safe? note that nΔs and Δs
            # are not referenced through the identity.() broadcast
            s=rand(trng,(1:nparticles)[Δs .<= Δs[i]])
        end
        a = s
        while a == s
            a = rand(trng,1:nparticles)
        end
        b = a
        while b == a || b == s
            b = rand(trng,1:nparticles)
        end
        # θp is a new Particle with new tuple values (.x) [see comment above]
        θp = op(+,θs[s],op(*,op(-,θs[a],θs[b]), γ))
        lπ = logpdf(prior, push_p(prior,θp.x))
        w_prior = lπ - logπ[i]
        log(rand(trng)) > min(0,w_prior) && continue
        nsims[i]+=1
        dp = dist!(θp.x, ve)
        if dp <= max(ϵ, Δs[i])
            nΔs[i] = dp
            nθs[i] = θp
            nlogπ[i] = lπ
        end
    end
end

function abcde!(prior, dist!, ϵ_target, varexternal; nparticles=50, generations=20, α=0, earlystop=false,
                verbose=true, rng=Random.GLOBAL_RNG, proposal_width=1.0,
                ex=ThreadedEx())
    @info("Running on experimental abcde! method")
    @info("Running abcde! with executor ", typeof(ex))

    ### this seems to be initialisation
    @assert 0<=α<1 "α must be in 0 <= α < 1."

    # draw prior parameters for each particle
    θs =[op(float, Particle(rand(rng, prior))) for i = 1:nparticles]
    logπ = [logpdf(prior, push_p(prior,θs[i].x)) for i = 1:nparticles]

    ve = deepcopy(varexternal)
    Δs = fill(dist!(θs[1].x, ve),nparticles)

    abcde_init!(prior, dist!, varexternal, θs, logπ, Δs, nparticles, rng, ex)
    ###

    ### actual ABC run
    nsims = zeros(Int,nparticles)
    γ = proposal_width*2.38/sqrt(2*length(prior))
    iters=0
    complete=1-sum(Δs.>ϵ_target)/nparticles
    while iters<generations
        iters+=1
        # identity.() behaves like deepcopy(), i.e. == is true, === is false
        # so there are n=new object that can be mutated without data races
        nθs = identity.(θs)
        nΔs = identity.(Δs)
        nlogπ=identity.(logπ)

        # returns minimal and maximal distance/cost
        ϵ_l, ϵ_h = extrema(Δs)
        if earlystop
            ϵ_h<=ϵ_target && break
        end
        ϵ_pop = max(ϵ_target,ϵ_l + α * (ϵ_h - ϵ_l))

        abcde_swarm!(prior, dist!, varexternal, θs, logπ, Δs, nθs, nlogπ, nΔs,
                            ϵ_pop, ϵ_target, γ, nparticles, nsims, earlystop, rng, ex)

        θs = nθs
        Δs = nΔs
        logπ = nlogπ
        ncomplete = 1 - sum(Δs .> ϵ_target) / nparticles
        if verbose && (ncomplete != complete || complete >= (nparticles - 1) / nparticles)
            @info "Finished run:" completion=ncomplete nsim = sum(nsims) range_ϵ = extrema(Δs)
        end
        complete=ncomplete
    end
    conv=maximum(Δs) <= ϵ_target
    if verbose
        @info "End:" completion = complete converged = conv nsim = sum(nsims) range_ϵ = extrema(Δs)
    end
    θs = [push_p(prior, θs[i].x) for i = 1:nparticles]
    l = length(prior)
    P = map(x -> Particles(x), getindex.(θs, i) for i = 1:l)
    length(P)==1 && (P=first(P))
    (P=P, C=Particles(Δs), reached_ϵ=conv)
end

export abcde!


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui

function naive_combo(store::AbstractStore, get_state=Redux.get_state)
    if ReduxImGui.Combo(store, s->get_state(s).combos["basic_combo"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
    end
    CImGui.SameLine()
    ShowHelpMarker("Refer to the \"Combo\" section below for an explanation of the full BeginCombo/EndCombo API, and demonstration of various flags.\n")
end


text:
#! /usr/bin/julia

# Rosetta Code, Zero to the zero power

zs = Any[zero(Complex),
         zero(FloatingPoint),
         zero(Rational),
         zero(Integer),
         zero(Bool)]

for i in zs, j in zs
    println(i, "^", j, " = ", i^j, " (", typeof(i^j), ")")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code