text:
module TestDefault
  using Test
  using BioEnergeticFoodWebs

  food_chain = [0 1 0 ; 0 0 1; 0 0 0]
  metab_status = [true, false, false]
  p = model_parameters(food_chain, vertebrates = metab_status)

  # DEFAULT ARG.
  # growth rates
  @test p[:r] == [1, 1, 1]
  # metabolic rates
  a_invertebrate = 0.3141
  a_vertebrate = 0.88
  a_producer = 0.138
  @test p[:x] == [a_vertebrate, a_invertebrate, a_producer]
  # maximum consumption rate
  y_vertebrate = 4.0
  y_invertebrate = 8.0
  y_producer = 0.0
  @test p[:y] == [y_vertebrate, y_invertebrate, y_producer]
  # handling time
  @test p[:ht] == 1 ./ p[:y]
  # attack rate
  @test p[:ar] == 1 ./ (0.5 * p[:ht])
  # half saturation constant
  hsc = 1 ./ (p[:ar] .* p[:ht])
  hsc[isnan.(hsc)] .= 0.0
  @test p[:Γ] == hsc

  # PASSED ARG.
  p = model_parameters(food_chain, vertebrates = metab_status,
                       handlingtime = NoEffectTemperature(:handlingtime, parameters_tuple = (y_vertebrate = 3.0, y_invertebrate = 7.0)),
                       attackrate = NoEffectTemperature(:attackrate, parameters_tuple = (Γ = 0.8,)),
                       metabolicrate = NoEffectTemperature(:metabolism, parameters_tuple = (a_vertebrate = 0.8, a_invertebrate = 0.3, a_producer = 0.1)),
                       growthrate = NoEffectTemperature(:growth, parameters_tuple = (ar = 2.0,)))
  @test p[:r] == [2.0, 2.0, 2.0]
  # metabolic rates
  a_invertebrate = 0.3
  a_vertebrate = 0.8
  a_producer = 0.1
  @test p[:x] == [a_vertebrate, a_invertebrate, a_producer]
  # maximum consumption rate
  y_vertebrate = 3.0
  y_invertebrate = 7.0
  y_producer = 0.0
  @test p[:y] == [y_vertebrate, y_invertebrate, y_producer]
  # handling time
  @test p[:ht] == 1 ./ p[:y]
  # attack rate
  @test p[:ar] == 1 ./ (0.8 * p[:ht])
  # half saturation constant
  hsc = 1 ./ (p[:ar] .* p[:ht])
  hsc[isnan.(hsc)] .= 0.0
  @test p[:Γ] == hsc
  @test p[:Γ] == [0.8, 0.8, 0.0]

  # NORMALISATION
  p_norm = model_parameters(food_chain, vertebrates = metab_status,
                       handlingtime = NoEffectTemperature(:handlingtime, parameters_tuple = (y_vertebrate = 3.0, y_invertebrate = 7.0)),
                       attackrate = NoEffectTemperature(:attackrate, parameters_tuple = (Γ = 0.8,)),
                       metabolicrate = NoEffectTemperature(:metabolism, parameters_tuple = (a_vertebrate = 0.8, a_invertebrate = 0.3, a_producer = 0.1)),
                       growthrate = NoEffectTemperature(:growth, parameters_tuple = (ar = 2.0,)),
                       scale_growth = true,
                       scale_maxcons = true,
                       scale_metabolism = true)
  rspp = 2.0
  @test p_norm[:r] == [2.0, 2.0, 2.0] ./ rspp
  # metabolic rates
  @test p_norm[:x] == [a_vertebrate, a_invertebrate, a_producer] ./ rspp
  # maximum consumption rate
  @test p_norm[:y] == [y_vertebrate, y_invertebrate, y_producer] ./ [a_vertebrate, a_invertebrate, a_producer]
  # handling time
  # attack rate
  @test p_norm[:ar] == 1 ./ (0.8 * p[:ht])
  # half saturation constant
  @test p_norm[:Γ] == hsc
  @test p_norm[:Γ] == [0.8, 0.8, 0.0]

  # Different temperatures, same rates
  temp2 = 293.15
  p2 = model_parameters(food_chain, T = temp2, vertebrates = metab_status,
                       handlingtime = NoEffectTemperature(:handlingtime, parameters_tuple = (y_vertebrate = 3.0, y_invertebrate = 7.0)),
                       attackrate = NoEffectTemperature(:attackrate, parameters_tuple = (Γ = 0.8,)),
                       metabolicrate = NoEffectTemperature(:metabolism, parameters_tuple = (a_vertebrate = 0.8, a_invertebrate = 0.3, a_producer = 0.1)),
                       growthrate = NoEffectTemperature(:growth, parameters_tuple = (ar = 2.0,)))
  @test p[:r] == p2[:r]
  @test p[:x] == p2[:x]
  @test p[:ht] == p2[:ht]
  @test p[:ar] == p2[:ar]

  # ERRORS

  @test_throws Exception model_parameters(food_chain, vertebrates = metab_status, handlingtime = NoEffectTemperature(:y))
  @test_throws Exception model_parameters(food_chain, vertebrates = metab_status, attackrate = NoEffectTemperature(:y))
  @test_throws Exception model_parameters(food_chain, vertebrates = metab_status, metabolicrate = NoEffectTemperature(:y))
  @test_throws Exception model_parameters(food_chain, vertebrates = metab_status, growthrate = NoEffectTemperature(:y))

end

module TestEppley
  using Test
  using BioEnergeticFoodWebs

  omnivory = [0 1 1 ; 0 0 1 ; 0 0 0]
  metabolic_status = [:true, :false, :false]
  bmass = [100.0, 10.0, 1.0]
  temp = 270.0

  # GROWTH - defaults
  p_g_def = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedEppley(:r))
  r_g_def = bmass .^ -0.25 .* 0.81 .* exp(0.0631 .* (temp.-273.15)) * (1 .- (((temp.-273.15) .- (298.15-273.15)) ./ (35/2)).^2)
  @test p_g_def[:r] == r_g_def

  # GROWTH - changing temperature
  temp_2 = 250.0
  p_g_temp = model_parameters(omnivory, T = temp_2, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedEppley(:r))
  r_g_temp = bmass .^ -0.25 .* 0.81 .* exp(0.0631 .* (temp_2.-273.15)) * (1 .- (((temp_2.-273.15) .- (298.15-273.15)) ./ (35/2)).^2)
  @test p_g_temp[:r] == r_g_temp

  # GROWTH - negative outside of range
  mintemp = 298.15-17.5
  maxtemp = 298.15+17.5
  p_r_min = model_parameters(omnivory, T = mintemp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedEppley(:r))
  p_r_max = model_parameters(omnivory, T = maxtemp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedEppley(:r))
  @test p_r_min[:r] == p_r_max[:r] == [0.0, 0.0, 0.0]
  inftemp = mintemp - 1
  suptemp = maxtemp + 1
  p_r_inf = model_parameters(omnivory, T = inftemp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedEppley(:r))
  p_r_sup = model_parameters(omnivory, T = suptemp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedEppley(:r))
  @test !any(x -> x > 0, p_r_inf[:r])

  # GROWTH - passed
  pt = (maxrate_0 = 0.6, eppley_exponent = 0.1, z = 215, β = -0.2, range = 20)
  p_g = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedEppley(:r, parameters_tuple = pt))
  r_g = bmass .^ -0.2 .* 0.6 .* exp(0.1 .* (temp.-273.15)) * (1 .- (((temp.-273.15) .- (215-273.15)) ./ (20/2)).^2)
  @test p_g[:r] == r_g

  # METABOLISM - default
  p_x_def = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedEppley(:x))
  @test p_x_def[:x] == p_g_def[:r]

  # METABOLISM - changing temperature
  p_x_temp = model_parameters(omnivory, T = temp_2, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedEppley(:x))
  @test p_x_temp[:x] == p_g_temp[:r]

  # METABOLISM - negative outside of range
  mintemp = 298.15-17.5
  maxtemp = 298.15+17.5
  p_x_min = model_parameters(omnivory, T = mintemp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedEppley(:metabolicrate))
  p_x_max = model_parameters(omnivory, T = maxtemp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedEppley(:metabolicrate))
  @test p_x_min[:x] == p_x_max[:x] == [0.0, 0.0, 0.0]
  inftemp = mintemp - 1
  suptemp = maxtemp + 1
  p_x_inf = model_parameters(omnivory, T = inftemp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedEppley(:metabolicrate))
  p_x_sup = model_parameters(omnivory, T = suptemp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedEppley(:metabolicrate))
  @test !any(x -> x > 0, p_x_inf[:x])

  # METABOLISM - passed arguments
  pt_x = (maxrate_0_producer = 0.8, maxrate_0_invertebrate = 0.6, maxrate_0_vertebrate = 5,
      eppley_exponent_producer = 0.1, eppley_exponent_invertebrate = 0.06, eppley_exponent_vertebrate = 0.09,
      z_producer = 298, z_invertebrate = 250, z_vertebrate = 270,
      range_producer = 30, range_invertebrate = 32, range_vertebrate = 33,
      β_producer = -0.2, β_invertebrate = -0.3, β_vertebrate = -0.27)
  p_x = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedEppley(:x, parameters_tuple = pt_x))
  maxrate_0_all = [5, 0.6, 0.8]
  eppley_exponent_all = [0.09, 0.06, 0.1]
  z_all = [270, 250, 298]
  z_all = z_all .- 273.15
  range_all = [33, 32, 30]
  β_all = [-0.27, -0.3, -0.2]
  x_x = bmass.^β_all .* maxrate_0_all .* exp.(eppley_exponent_all .* (temp.-273.15)) .* (1 .- (((temp.-273.15) .- z_all) ./ (range_all./2)).^2)
  @test x_x == p_x[:x]

  # ERRORS
  @test_throws Exception model_parameters(omnivory, metabolicrate = ExtendedEppley(:handlingtime))
  @test_throws Exception model_parameters(omnivory, metabolicrate = ExtendedEppley(:attackrate))
  @test_throws Exception model_parameters(omnivory, metabolicrate = ExtendedEppley(:y))

end

module TestExponentialBA
  using Test
  using BioEnergeticFoodWebs

  omnivory = [0 1 1 ; 0 0 1 ; 0 0 0]
  metabolic_status = [:true, :false, :false]
  bmass = [100.0, 10.0, 1.0]
  temp = 270.0

  #GROWTH
  #defaults
  p_r_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExponentialBA(:r))
  k = 8.617e-5
  r_d = exp(-15.68) .* 4e6 .* (bmass .^-0.25) .* exp.(-0.84 .* (293.15 .- temp) ./ (k .* temp .* 293.15))
  @test p_r_d[:r] == r_d
  #change temperature
  temp2 = 250.0
  p_r_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExponentialBA(:r))
  r_t = exp(-15.68) .* 4e6 .* (bmass .^-0.25) .* exp.(-0.84 .* (293.15 .- temp2) ./ (k .* temp2 .* 293.15))
  @test p_r_t[:r] == r_t
  #rate increases with temperature
  @test !any(p_r_t[:r] .> p_r_d[:r])
  #passed arguments
  pt_r = (norm_constant = exp(-18), activation_energy = -0.8, T0 = 290, β = -0.31)
  p_r_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExponentialBA(:r, parameters_tuple = pt_r))
  r_2 = exp(-18) .* (bmass .^-0.31) .* exp.(-0.8 .* (290 .- temp) ./ (k .* temp .* 290))
  @test p_r_2[:r] == r_2

  #METABOLISM
  #defaults
  p_x_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExponentialBA(:x))
  x_d = exp(-16.54) .* 4e6 .* (bmass .^-0.31) .* exp.(-0.69 .* (293.15 .- temp) ./ (k .* temp .* 293.15))
  @test p_x_d[:x] == x_d
  #change temperature
  p_x_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExponentialBA(:x))
  x_t = exp(-16.54) .* 4e6 .* (bmass .^-0.31) .* exp.(-0.69 .* (293.15 .- temp2) ./ (k .* temp2 .* 293.15))
  @test p_x_t[:x] == x_t
  #rate increases with temperature
  @test !any(p_x_t[:x] .> p_x_d[:x])
  #passed arguments
  pt_x = (norm_constant_producer = exp(-16), norm_constant_invertebrate = exp(-17), norm_constant_vertebrate = exp(-18),
             activation_energy_producer = -0.6, activation_energy_invertebrate = -0.7, activation_energy_vertebrate = -0.8,
             T0_producer = 270, T0_invertebrate = 280, T0_vertebrate = 290,
             β_producer = -0.2, β_invertebrate = -0.3, β_vertebrate = -0.4)
  norm_constant_all = [-18, -17, -16]
  activation_energy_all = [-0.8, -0.7, -0.6]
  T0_all = [290, 280, 270]
  β_all = [-0.4, -0.3, -0.2]
  p_x_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExponentialBA(:x, parameters_tuple = pt_x))
  x_2 = exp.(norm_constant_all) .* (bmass .^β_all) .* exp.(activation_energy_all .* (T0_all .- temp) ./ (k .* temp .* T0_all))
  @test p_x_2[:x] == x_2

  #ATTACK
  #defaults
  p_ar_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExponentialBA(:attackrate))
  ar_d = exp.([-13.1, -13.1, 0.0]) .* 4e6 .* (bmass .^[-0.8, -0.8, 0.0]) .* (bmass' .^[0 -0.8 0.25; 0 0 0.25 ; 0 0 0]) .* exp.([-0.38, -0.38, 0.0] .* ([293.15, 293.15, 0.0] .- temp) ./ (k .* temp .* [293.15, 293.15, 0.0]))
  ar_d[isnan.(ar_d)] .= 0
  @test p_ar_d[:ar] == ar_d
  #change temperature
  p_ar_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExponentialBA(:attackrate))
  ar_t = exp.([-13.1, -13.1, 0.0]) .* 4e6 .* (bmass .^[-0.8, -0.8, 0.0]) .* (bmass' .^[0 -0.8 0.25; 0 0 0.25 ; 0 0 0]) .* exp.([-0.38, -0.38, 0.0] .* ([293.15, 293.15, 0.0] .- temp2) ./ (k .* temp2 .* [293.15, 293.15, 0.0]))
  ar_t[isnan.(ar_t)] .= 0
  @test p_ar_t[:ar] == ar_t
  #rate increases with temperature
  @test !any(p_ar_t[:ar] .> p_ar_d[:ar])
  #passed arguments
  pt_ar = (norm_constant_vertebrate = exp(-12), norm_constant_invertebrate = exp(-14),
  						activation_energy_vertebrate = -0.3, activation_energy_invertebrate = -0.4,
  						T0_vertebrate = 290, T0_invertebrate = 270,
  						β_producer = 0.2, β_vertebrate = -0.9, β_invertebrate = -0.7)
  p_ar_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExponentialBA(:attackrate, parameters_tuple = pt_ar))
  ar_2 = exp.([-12, -14, 0.0]) .* (bmass .^[-0.9, -0.7, 0.0]) .* (bmass' .^[0 -0.7 0.2; 0 0 0.2 ; 0 0 0]) .* exp.([-0.3, -0.4, 0.0] .* ([290, 270, 0.0] .- temp) ./ (k .* temp .* [290, 270, 0.0]))
  ar_2[isnan.(ar_2)] .= 0
  @test p_ar_2[:ar] == ar_2

  #HANDLING
  #defaults
  p_ht_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, handlingtime = ExponentialBA(:handlingtime))
  ht_d = exp.([9.66, 9.66, 0.0]) .* 4e6 .* (bmass .^[0.47, 0.47, 0.0]) .* (bmass' .^[0.0 0.47 -0.45 ; 0 0 -0.45; 0 0 0]) .* exp.([0.26, 0.26, 0.0] .* ([293.15, 293.15, 0.0] .- temp) ./ (k .* temp .* [293.15, 293.15, 0.0]))
  ht_d[isnan.(ht_d)] .= 0
  @test p_ht_d[:ht] == ht_d
  #change temperature
  p_ht_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, handlingtime = ExponentialBA(:handlingtime))
  ht_t = exp.([9.66, 9.66, 0.0]) .* 4e6 .* (bmass .^[0.47, 0.47, 0.0]) .* (bmass' .^[0.0 0.47 -0.45 ; 0 0 -0.45; 0 0 0]) .* exp.([0.26, 0.26, 0.0] .* ([293.15, 293.15, 0.0] .- temp2) ./ (k .* temp2 .* [293.15, 293.15, 0.0]))
  ht_t[isnan.(ht_t)] .= 0
  @test p_ht_t[:ht] == ht_t
  #rate decreases with temperature
  @test !any(p_ht_t[:ht] .< p_ht_d[:ht])
  #passed arguments
  pt_ht = (norm_constant_vertebrate = exp(9), norm_constant_invertebrate = exp(10),
  						activation_energy_vertebrate = 0.2, activation_energy_invertebrate = 0.3,
  						T0_vertebrate = 290, T0_invertebrate = 270,
  						β_producer = -0.4, β_vertebrate = 0.3, β_invertebrate = 0.5)
  p_ht_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, handlingtime = ExponentialBA(:handlingtime, parameters_tuple = pt_ht))
  ht_2 = exp.([9, 10, 0.0]) .* (bmass .^[0.3, 0.5, 0.0]) .* (bmass' .^[0.0 0.5 -0.4 ; 0 0 -0.4; 0 0 0]) .* exp.([0.2, 0.3, 0.0] .* ([290, 270, 0.0] .- temp) ./ (k .* temp .* [290, 270, 0.0]))
  ht_2[isnan.(ht_2)] .= 0
  @test p_ht_2[:ht] == ht_2

  #ERRORS
  @test_throws Exception model_parameters(omnivory, metabolicrate = ExponentialBA(:y))

end

module TestExtendedBA
  using Test
  using BioEnergeticFoodWebs

  omnivory = [0 1 1 ; 0 0 1 ; 0 0 0]
  metabolic_status = [:true, :false, :false]
  bmass = [100.0, 10.0, 1.0]
  temp = 270.0
  temp2 = 250.0
  k = 8.617e-5 # Boltzmann constant

  #GROWTH
  #defaults
  p_r_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedBA(:r))
  Δenergy = 1.15 - 0.53
  r_d = 1.8e9 .* bmass .^(-0.25) .* exp.(.-0.53 ./ (k * temp)) .* (1 ./ (1 + exp.(-1 / (k * temp) .* (1.15 .- (1.15 ./ 298.15 .+ k .* log(0.53 ./ Δenergy)).*temp))))
  @test p_r_d[:r] == r_d
  #change temperature
  p_r_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedBA(:r))
  @test p_r_d[:r] != p_r_t[:r]
  #passed arguments
  pt_r = (norm_constant = 4000,)
  p_r_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedBA(:r, parameters_tuple = pt_r))
  r_2 = (r_d ./ 1.8e9) .* 4000
  @test p_r_2[:r] ≈ r_2 atol=1e-6
  #max rate at Topt
  p_r_inf = model_parameters(omnivory, T = 297.0, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedBA(:r))
  p_r_sup = model_parameters(omnivory, T = 299.0, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedBA(:r))
  p_r_opt = model_parameters(omnivory, T = 298.15, bodymass = bmass, vertebrates = metabolic_status, growthrate = ExtendedBA(:r))
  @test !any(p_r_opt[:r] .< p_r_sup[:r])
  @test !any(p_r_opt[:r] .< p_r_inf[:r])

  #METABOLISM
  #defaults
  p_x_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedBA(:x))
  x_d = 3e8 .* bmass .^(-0.25) .* exp.(.-[0.53, 0.53, 0.53] ./ (k * temp)) .* (1 ./ (1 .+ exp.(-1 ./ (k * temp) .* (1.15 .- (1.15 ./ 298.15 .+ k .* log.(0.53 ./ Δenergy)).* temp))))
  @test p_x_d[:x] == x_d
  #change temperature
  p_x_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedBA(:x))
  @test p_x_t[:x] != p_x_d[:x]
  #passed arguments
  pt_x = (β_vertebrate = -0.5,)
  p_x_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedBA(:x, parameters_tuple = pt_x))
  p_x_2[:x] != p_x_d[:x]
  x_2 = 3e8 .* bmass .^([-0.5, -0.25, -0.25]) .* exp.(.-[0.53, 0.53, 0.53] ./ (k * temp)) .* (1 ./ (1 .+ exp.(-1 ./ (k * temp) .* (1.15 .- (1.15 ./ 298.15 .+ k .* log.(0.53 ./ Δenergy)).* temp))))
  @test p_x_2[:x] == x_2
  #max rate at Topt
  p_x_inf = model_parameters(omnivory, T = 297.0, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedBA(:x))
  p_x_sup = model_parameters(omnivory, T = 299.0, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedBA(:x))
  p_x_opt = model_parameters(omnivory, T = 298.15, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = ExtendedBA(:x))
  @test !any(p_x_opt[:x] .< p_x_sup[:x])
  @test !any(p_x_opt[:x] .< p_x_inf[:x])

  #ATTACK
  #defaults
  p_ar_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:attackrate))
  Δenergy = 1.15 - 0.8
  ar_d = [5e13, 5e13, 5e13] .* bmass .^([0.25, 0.25, 0.0]) .* bmass' .^([0.0 0.25 0.25 ; 0.0 0.0 0.25 ; 0.0 0.0 0.0]) .* exp.(.-[0.8, 0.8, 0.8] ./ (k * temp)) .* (1 ./ (1 .+ exp.(-1 / (k * temp) .* ([1.15, 1.15, 0.0] .- ([1.15, 1.15, 0.0] ./ [298.15, 298.15, 0.0] .+ k .* log.([0.8, 0.8, 0.0] ./ Δenergy)) .* temp))))
  ar_d[isnan.(ar_d)] .= 0
  @test p_ar_d[:ar] == ar_d
  #change temperature
  p_ar_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:attackrate))
  @test p_ar_d[:ar] != p_ar_t[:ar]
  #passed arguments
  pt_ar = (norm_constant_invertebrate = 3e7,)
  p_ar_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:attackrate, parameters_tuple = pt_ar))
  @test p_ar_2[:ar] != p_ar_d[:ar]
  ar_2 = [5e13, 3e7, 5e13] .* bmass .^([0.25, 0.25, 0.0]) .* bmass' .^([0.0 0.25 0.25 ; 0.0 0.0 0.25 ; 0.0 0.0 0.0]) .* exp.(.-[0.8, 0.8, 0.8] ./ (k * temp)) .* (1 ./ (1 .+ exp.(-1 / (k * temp) .* ([1.15, 1.15, 1.15] .- ([1.15, 1.15, 1.15] ./ [298.15, 298.15, 0.0] .+ k .* log.([0.8, 0.8, 0.8] ./ Δenergy)) .* temp))))
  @test p_ar_2[:ar] == ar_2
  #max rate at Topt
  p_ar_inf = model_parameters(omnivory, T = 297.0, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:attackrate))
  p_ar_sup = model_parameters(omnivory, T = 299.0, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:attackrate))
  p_ar_opt = model_parameters(omnivory, T = 298.15, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:attackrate))
  @test !any(p_ar_opt[:ar] .< p_ar_sup[:ar])
  @test !any(p_ar_opt[:ar] .< p_ar_inf[:ar])

  #ERRORS
  @test_throws Exception model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:handlingtime))
  @test_throws Exception model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = ExtendedBA(:y))

end

module TestGaussian
  using Test
  using BioEnergeticFoodWebs

  omnivory = [0 1 1 ; 0 0 1 ; 0 0 0]
  metabolic_status = [:true, :false, :false]
  bmass = [100.0, 10.0, 1.0]
  temp = 270.0
  temp2 = 250.0

  #GROWTH
  #defaults
  p_r_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = Gaussian(:r))
  r_d = bmass .^ -0.25 .* 1.0 .* exp( .- (temp .- 298.15) .^ 2 ./ (2 .* 20 .^ 2))
  @test p_r_d[:r] == r_d
  #change temperature
  p_r_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, growthrate = Gaussian(:r))
  @test p_r_t[:r] != p_r_d[:r]
  #passed arguments
  pt_r = (norm_constant = 0.7,)
  p_r_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, growthrate = Gaussian(:r, parameters_tuple = pt_r))
  @test p_r_2[:r] != p_r_d[:r]
  r_2 = bmass .^ -0.25 .* 0.7 .* exp( .- (temp .- 298.15) .^ 2 ./ (2 .* 20 .^ 2))
  @test p_r_2[:r] == r_2
  #max rate at Topt
  p_r_inf = model_parameters(omnivory, T = 297.0, bodymass = bmass, vertebrates = metabolic_status, growthrate = Gaussian(:r))
  p_r_sup = model_parameters(omnivory, T = 299.0, bodymass = bmass, vertebrates = metabolic_status, growthrate = Gaussian(:r))
  p_r_opt = model_parameters(omnivory, T = 298.15, bodymass = bmass, vertebrates = metabolic_status, growthrate = Gaussian(:r))
  @test !any(p_r_opt[:r] .< p_r_sup[:r])
  @test !any(p_r_opt[:r] .< p_r_inf[:r])

  #METABOLISM
  #defaults
  p_x_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = Gaussian(:x))
  x_d = bmass .^ -0.25 .* [0.9, 0.35, 0.2] .* exp.( .- (temp .- 298.15) .^ 2 ./ (2 .* 20 .^ 2))
  @test p_x_d[:x] == x_d
  #change temperature
  p_x_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = Gaussian(:x))
  @test p_x_t[:x] != p_x_d[:x]
  #passed arguments
  pt_x = (range_producer = 30,)
  p_x_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = Gaussian(:x, parameters_tuple = pt_x))
  @test p_x_2[:x] != p_x_d[:x]
  x_2 = bmass .^ -0.25 .* [0.9, 0.35, 0.2] .* exp.( .- (temp .- 298.15) .^ 2 ./ (2 .* [20,20,30] .^ 2))
  @test p_x_2[:x] == x_2
  #max rate at Topt
  p_x_inf = model_parameters(omnivory, T = 297.0, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = Gaussian(:x))
  p_x_sup = model_parameters(omnivory, T = 299.0, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = Gaussian(:x))
  p_x_opt = model_parameters(omnivory, T = 298.15, bodymass = bmass, vertebrates = metabolic_status, metabolicrate = Gaussian(:x))
  @test !any(p_x_opt[:x] .< p_x_sup[:x])
  @test !any(p_x_opt[:x] .< p_x_inf[:x])

  #ATTACK
  #defaults
  p_ar_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = Gaussian(:attackrate))
  ar_d = bmass .^ [-0.25,-0.25,0.0] .* bmass' .^ [0.0 -0.25 -0.25 ; 0.0 0.0 -0.25 ; 0.0 0.0 0.0] .* [16, 16, 0.0] .* exp.( .- (temp .- [295,295,0]) .^ 2 ./ (2 .* [20,20,0] .^ 2))
  ar_d[isnan.(ar_d)] .= 0
  @test p_ar_d[:ar] == ar_d
  #change temperature
  p_ar_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, attackrate = Gaussian(:attackrate))
  @test p_ar_t[:ar] != p_ar_d[:ar]
  #passed arguments
  pt_ar = (T_opt_invertebrate = 270,)
  p_ar_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, attackrate = Gaussian(:attackrate, parameters_tuple = pt_ar))
  @test p_ar_2[:ar] != p_ar_d[:ar]
  ar_2 = bmass .^ [-0.25,-0.25,0.0] .* bmass' .^ [0.0 -0.25 -0.25 ; 0.0 0.0 -0.25 ; 0.0 0.0 0.0] .* [16, 16, 0.0] .* exp.( .- (temp .- [295,270,0]) .^ 2 ./ (2 .* [20,20,0] .^ 2))
  ar_2[isnan.(ar_2)] .= 0
  @test p_ar_2[:ar] == ar_2
  #max rate at Topt
  p_ar_inf = model_parameters(omnivory, T = 293.0, bodymass = bmass, vertebrates = metabolic_status, attackrate = Gaussian(:attackrate))
  p_ar_sup = model_parameters(omnivory, T = 297.0, bodymass = bmass, vertebrates = metabolic_status, attackrate = Gaussian(:attackrate))
  p_ar_opt = model_parameters(omnivory, T = 295.0, bodymass = bmass, vertebrates = metabolic_status, attackrate = Gaussian(:attackrate))
  @test !any(p_ar_opt[:ar] .< p_ar_sup[:ar])
  @test !any(p_ar_opt[:ar] .< p_ar_inf[:ar])

  #HANDLING
  #defaults
  p_ht_d = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, handlingtime = Gaussian(:handlingtime))
  ht_d = bmass .^ [-0.25,-0.25,0.0] .* bmass' .^ [0.0 -0.25 -0.25 ; 0.0 0.0 -0.25 ; 0.0 0.0 0.0] .* [0.5, 0.5, 0.0] .* exp.((temp .- [295,295,0]) .^ 2 ./ (2 .* [20,20,0] .^ 2))
  ht_d[isnan.(ht_d)] .= 0
  @test p_ht_d[:ht] == ht_d
  #change temperature
  p_ht_t = model_parameters(omnivory, T = temp2, bodymass = bmass, vertebrates = metabolic_status, handlingtime = Gaussian(:handlingtime))
  @test p_ht_t[:ht] != p_ht_d[:ht]
  #passed arguments
  pt_ht = (T_opt_invertebrate = 270,)
  p_ht_2 = model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, handlingtime = Gaussian(:handlingtime, parameters_tuple = pt_ht))
  @test p_ht_2[:ht] != p_ht_d[:ht]
  ht_2 = bmass .^ [-0.25,-0.25,0.0] .* bmass' .^ [0.0 -0.25 -0.25 ; 0.0 0.0 -0.25 ; 0.0 0.0 0.0] .* [0.5, 0.5, 0.0] .* exp.((temp .- [295,270,0]) .^ 2 ./ (2 .* [20,20,0] .^ 2))
  ht_2[isnan.(ht_2)] .= 0
  @test p_ht_2[:ht] == ht_2
  #min rate at Topt
  p_ht_inf = model_parameters(omnivory, T = 290.0, bodymass = bmass, vertebrates = metabolic_status, handlingtime = Gaussian(:handlingtime))
  p_ht_sup = model_parameters(omnivory, T = 305.0, bodymass = bmass, vertebrates = metabolic_status, handlingtime = Gaussian(:handlingtime))
  p_ht_opt = model_parameters(omnivory, T = 298.15, bodymass = bmass, vertebrates = metabolic_status, handlingtime = Gaussian(:handlingtime))
  @test !any(p_ht_opt[:ht] .> p_ht_sup[:ht])
  @test !any(p_ht_opt[:ht] .> p_ht_inf[:ht])

  #ERRORS
  @test_throws Exception model_parameters(omnivory, T = temp, bodymass = bmass, vertebrates = metabolic_status, handlingtime = Gaussian(:y))

end

ext: jl


text:
mutable struct BehaviourActionValues <: LoggerKeyData
    max_q::Array{Float64}
    function BehaviourActionValues(logger_init_info)
        new([])
    end
end

function lg_step!(self::BehaviourActionValues, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    base_s = proc_input(agent, s)
    q = maximum(agent.behaviour_learner(base_s))
    push!(self.max_q,q)
end

function lg_episode_end!(self::BehaviourActionValues, cur_step_in_episode, cur_step_total)
end

function save_log(self::BehaviourActionValues, save_dict::Dict)
    save_dict[:max_behaviour_q] = self.max_q
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using SegyIO
using Test

include("test_read.jl")
include("test_SeisBlock.jl")
include("test_SeisCon.jl")
include("test_write.jl")


text:
immutable dampCoef
    dVxx  :: SparseMatrixCSC{Float64, Int64}
    dVxz  :: SparseMatrixCSC{Float64, Int64}
    dVzx  :: SparseMatrixCSC{Float64, Int64}
    dVzz  :: SparseMatrixCSC{Float64, Int64}
    dTxxx :: SparseMatrixCSC{Float64, Int64}
    dTxxz :: SparseMatrixCSC{Float64, Int64}
    dTzzx :: SparseMatrixCSC{Float64, Int64}
    dTzzz :: SparseMatrixCSC{Float64, Int64}
    dTxzx :: SparseMatrixCSC{Float64, Int64}
    dTxzz :: SparseMatrixCSC{Float64, Int64}
end

function InitDampCoef(nz::Int64, nx::Int64, ext::Int64, iflag::Int64)
  if iflag == 1
     Nz = nz + 2*ext
  elseif iflag == 2
     Nz = nz +   ext
  end
  Nx = nx + 2*ext
  dpCoef = dampCoef(spzeros(Nz,Nx), spzeros(Nz,Nx), spzeros(Nz,Nx), spzeros(Nz,Nx), spzeros(Nz,Nx),
                    spzeros(Nz,Nx), spzeros(Nz,Nx), spzeros(Nz,Nx), spzeros(Nz,Nx), spzeros(Nz,Nx))
  return dpCoef
end

function DampBound(nz::Int64, nx::Int64, maxv::Float64, ext::Int64, rho::Array{Float64,2}, iflag::Int64, dz::Float64)
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    if ext == 5
       R = 0.01
    elseif ext == 10
       R = 0.001
    elseif ext >= 20
       R = 0.0001
    end
    ll = 2
    dVxx  = spzeros(Nz, Nx)
    dVxz  = spzeros(Nz, Nx)
    dVzx  = spzeros(Nz, Nx)
    dVzz  = spzeros(Nz, Nx)
    dTxxx = spzeros(Nz, Nx)
    dTxxz = spzeros(Nz, Nx)
    dTzzx = spzeros(Nz, Nx)
    dTzzz = spzeros(Nz, Nx)
    dTxzx = spzeros(Nz, Nx)
    dTxzz = spzeros(Nz, Nx)
    for i = 1 : ext
        tmph = -1.5 * maxv/(ext*dz) * log(R) * (((ext+0.5-i)/ext)^ll)
        tmpo = -1.5 * maxv/(ext*dz) * log(R) * (((ext+1.0-i)/ext)^ll)
        # upper
        if iflag == 1
           dVxz[i,:]  = tmpo * rho[i,:]
           dVzz[i,:]  = tmph * rho[i,:]
           dTxxz[i,:] = tmpo * ones(Nx)
           dTzzz[i,:] = tmpo * ones(Nx)
           dTxzz[i,:] = tmph * ones(Nx)
        end
        # left
        dVxx[:,i]  = tmph * rho[:,i]
        dVzx[:,i]  = tmpo * rho[:,i]
        dTxxx[:,i] = tmpo * ones(Nz)
        dTzzx[:,i] = tmpo * ones(Nz)
        dTxzx[:,i] = tmph * ones(Nz)
    end
    # lower side
    if iflag == 1
       for i = nz+ext+1 : Nz
           tmph = -1.5 * maxv/(ext*dz) * log(R) * (((i+0.5-nz-ext)/ext)^ll)
           tmpo = -1.5 * maxv/(ext*dz) * log(R) * (((i    -nz-ext)/ext)^ll)
           dVxz[i,:]  = tmpo * rho[i,:]
           dVzz[i,:]  = tmph * rho[i,:]
           dTxxz[i,:] = tmpo * ones(Nx)
           dTzzz[i,:] = tmpo * ones(Nx)
           dTxzz[i,:] = tmph * ones(Nx)
       end
    elseif iflag == 2
        for i = nz+1 : Nz
            tmph = -1.5 * maxv/(ext*dz) * log(R) * (((i+0.5-nz)/ext)^ll)
            tmpo = -1.5 * maxv/(ext*dz) * log(R) * (((i    -nz)/ext)^ll)
            dVxz[i,:]  = tmpo * rho[i,:]
            dVzz[i,:]  = tmph * rho[i,:]
            dTxxz[i,:] = tmpo * ones(Nx)
            dTzzz[i,:] = tmpo * ones(Nx)
            dTxzz[i,:] = tmph * ones(Nx)
        end
    end
    # right side
    for i = nx+ext+1 : Nx
        tmph = -1.5 * maxv/(ext*dz) * log(R) * (((i+0.5-nx-ext)/ext)^ll)
        tmpo = -1.5 * maxv/(ext*dz) * log(R) * (((i    -nx-ext)/ext)^ll)
        dVxx[:,i]  = tmph * rho[:,i]
        dVzx[:,i]  = tmpo * rho[:,i]
        dTxxx[:,i] = tmpo * ones(Nz)
        dTzzx[:,i] = tmpo * ones(Nz)
        dTxzx[:,i] = tmph * ones(Nz)
    end
    dpCoef = dampCoef(dVxx, dVxz, dVzx, dVzz, dTxxx, dTxxz, dTzzx, dTzzz, dTxzx, dTxzz)
    return dpCoef
end

function modExpand(par::Array{Float64,2}, ext::Int64, iflag::Int64)
    temp1 = repmat(par[:,  1], 1, ext)
    temp2 = repmat(par[:,end], 1, ext)
    par   = hcat(temp1, par, temp2)
    if iflag == 1
       temp1 = repmat((par[1,  :])', ext, 1)
       temp2 = repmat((par[end,:])', ext, 1)
       par   = vcat(temp1, par, temp2)
    elseif iflag == 2
       temp2 = repmat((par[end,:])', ext, 1)
       par   = vcat(par, temp2)
    end
    return par
end

function modelPadding(lambda::Array{Float64,2}, mu::Array{Float64,2}, rho::Array{Float64,2}, ext::Int64, iflag::Int64)
    lambda = modExpand(lambda, ext, iflag)
    mu     = modExpand(mu    , ext, iflag)
    rho    = modExpand(rho   , ext, iflag)
    return lambda, mu, rho
end

function modSmooth(par::Array{Float64,2}, L::Int64)
    n = 2*L + 1
    par1 = modExpand(par, L, 1)
    s = convert(Array{typeof(par1[1]), 1}, hanning(n))
    s = s / sum(s)
    par1 = conv2(s,s,par1)
    par1 = par1[2*L+1:end-2*L, 2*L+1:end-2*L]
    return par1
end

# function DampBound(nz::Int64, nx::Int64, maxv::Float64, ext::Int64, rho::Array{Float64,2}, iflag::Int64, dz::Float64)
#     if iflag == 1
#        Nz = nz + 2*ext
#     elseif iflag == 2
#        Nz = nz +   ext
#     end
#     Nx = nx + 2*ext
#     if ext == 5
#        R = 0.01
#     elseif ext == 10
#        R = 0.001
#     elseif ext >= 20
#        R = 0.0001
#     end
#     m =0.25; n = 0.75; ll = 4;
#     dVxx  = spzeros(Nz, Nx)
#     dVxz  = spzeros(Nz, Nx)
#     dVzx  = spzeros(Nz, Nx)
#     dVzz  = spzeros(Nz, Nx)
#     dTxxx = spzeros(Nz, Nx)
#     dTxxz = spzeros(Nz, Nx)
#     dTzzx = spzeros(Nz, Nx)
#     dTzzz = spzeros(Nz, Nx)
#     dTxzx = spzeros(Nz, Nx)
#     dTxzz = spzeros(Nz, Nx)
#     #  upper and left sides
#     for i = 1 : ext
#         tmph = -maxv/ext * log(R) * (m*(ext+0.5-i)/ext + n*((ext+0.5-i)/ext)^ll)
#         tmpo = -maxv/ext * log(R) * (m*(ext+1.0-i)/ext + n*((ext+1.0-i)/ext)^ll)
#         if iflag == 1
#            dVxz[i,:]  = tmpo * rho[i,:]
#            dVzz[i,:]  = tmph * rho[i,:]
#            dTxxz[i,:] = tmpo * ones(Nx)
#            dTzzz[i,:] = tmpo * ones(Nx)
#            dTxzz[i,:] = tmph * ones(Nx)
#         end
#         dVxx[:,i]  = tmph * rho[:,i]
#         dVzx[:,i]  = tmpo * rho[:,i]
#         dTxxx[:,i] = tmpo * ones(Nz)
#         dTzzx[:,i] = tmpo * ones(Nz)
#         dTxzx[:,i] = tmph * ones(Nz)
#     end
#     # lower side
#     if iflag == 1
#        for i = nz+ext+1 : Nz
#            tmph = -maxv/ext * log(R) * (m*(i+0.5-nz-ext)/ext + n*((i+0.5-nz-ext)/ext)^ll)
#            tmpo = -maxv/ext * log(R) * (m*(i    -nz-ext)/ext + n*((i    -nz-ext)/ext)^ll)
#            dVxz[i,:]  = tmpo * rho[i,:]
#            dVzz[i,:]  = tmph * rho[i,:]
#            dTxxz[i,:] = tmpo * ones(Nx)
#            dTzzz[i,:] = tmpo * ones(Nx)
#            dTxzz[i,:] = tmph * ones(Nx)
#        end
#     elseif iflag == 2
#        for i = nz+1 : Nz
#            tmph = -maxv/ext * log(R) * (m*(i+0.5-nz-ext)/ext + n*((i+0.5-nz-ext)/ext)^ll)
#            tmpo = -maxv/ext * log(R) * (m*(i    -nz-ext)/ext + n*((i    -nz-ext)/ext)^ll)
#            dVxz[i,:]  = tmpo * rho[i,:]
#            dVzz[i,:]  = tmph * rho[i,:]
#            dTxxz[i,:] = tmpo * ones(Nx)
#            dTzzz[i,:] = tmpo * ones(Nx)
#            dTxzz[i,:] = tmph * ones(Nx)
#        end
#     end
#     # right side
#     for i = nx+ext+1 : Nx
#         tmph = -maxv/ext * log(R) * (m*(i+0.5-nx-ext)/ext + n*((i+0.5-nx-ext)/ext)^ll)
#         tmpo = -maxv/ext * log(R) * (m*(i    -nx-ext)/ext + n*((i    -nx-ext)/ext)^ll)
#         dVxx[:,i]  = tmph * rho[:,i]
#         dVzx[:,i]  = tmpo * rho[:,i]
#         dTxxx[:,i] = tmpo * ones(Nz)
#         dTzzx[:,i] = tmpo * ones(Nz)
#         dTxzx[:,i] = tmph * ones(Nz)
#     end
#     dpCoef = dampCoef(dVxx, dVxz, dVzx, dVzz, dTxxx, dTxxz, dTzzx, dTzzz, dTxzx, dTxzz)
#     return dpCoef
# end

max_stars_repo_name: WneleiGao/ElasticWave


text:
# Arvind Prasadan
# 2019 January
# DMF Package
# SOBI Function (R) Wrapper

"""
	 SOBI_Wrapper(X, lags)

A wrapper function for the Second Order Blind Identification (SOBI) Source Separation algorithm. This is a wrapper that calls the function from R, and depends on the presence of an R install as well as the jointDiag package being installed.

# Arguments
- `X`: Data matrix, p dimensions/sensors x n samples
- `lags`: List of (Integral) lags for the algorithm to operate at

# Outputs
- `Q_hat`: Eigenvector Estimates
- `S_hat: The recovered latent time series
"""
function SOBI_Wrapper(X, lags)

    n = size(X, 2)
    p = size(X, 1)

    # Sanity check lags
    lags = abs.(Int.(round.(lags)))
    if 1 < length(lags)
        lags = [min(x, n - 1) for x in lags]
        lags = sort(unique(lags))
    end

    # Whiten X
    mu = mean(X; dims = 2)
    X = X .- mu
    cov_X = X * X'
    cov_X = (cov_X + cov_X') / 2.0
    cW, cV = eigen(cov_X)
    cov_invsqrt = cV * Diagonal(sqrt.(abs.(pinv(cW)))) * cV'
    cov_sqrt = cV * Diagonal(sqrt.(abs.(cW))) * cV'
    Xw = cov_invsqrt * X

    # Store matrices at different lags
    A = zeros(p, p, length(lags))
    for ll = 1:1:length(lags)
	    A[:, :, ll] = Xw[:, (lags[ll] + 1):end] * pinv(Xw[:, 1:(end - lags[ll])]) + Diagonal(ones(p)) # Regularizer
    end

    # Put Variables in R
    R"library(jointDiag)"
    R"A = $(A)"

    # Run SOBI: Joint Diagonalization
    Q_hat = reval("ajd(A, method = \"jedi\")\$A")
    Q_hat = convert(Array{Float64, 2}, Q_hat)
    Q_hat = cov_sqrt * Q_hat

    S_hat = (pinv(Q_hat) * (X .+ mu))'

    return (Q_hat, S_hat)
end


text:
@testset "Transfer operator from triangulations" begin
    # Embeddings
    E = StateSpaceReconstruction.cembed([diff(rand(15)) for i = 1:3])
    E_invariant = invariantize(E)

    # Triangulations
    triang = triangulate(E)
    triang_inv = triangulate(E_invariant)

    # Transfer operators from *invariant* triangulations
    TO = transferoperator_triang(triang_inv)
    TO_approx = transferoperator_triang(triang_inv, exact = false, parallel = false)
    TO_approx_rand = transferoperator_triang(triang_inv, exact = false, parallel = false, sample_randomly = true)
    #TO_exact = transferoperator(triang_inv, exact = true, parallel = false)
    #TO_exact_p = transferoperator(triang_inv, exact = true, parallel = true)

    @test typeof(TO) == ApproxSimplexTransferOperator #default = approx
    @test typeof(TO_approx) == ApproxSimplexTransferOperator
    @test typeof(TO_approx_rand) == ApproxSimplexTransferOperator
    #@test typeof(TO_exact) == ExactSimplexTransferOperator
    #@test typeof(TO_exact_p) == ExactSimplexTransferOperator

    #@test all(TO_exact.transfermatrix .== TO_exact_p.transfermatrix)

    # Approximations to intersections
    @test is_markov(TO)
    @test is_markov(TO_approx)
    @test is_markov(TO_approx_rand)
    @test is_almost_markov(TO)
    @test is_almost_markov(TO_approx)
    @test is_almost_markov(TO_approx_rand)

    # Exact intersections

    #@test is_almost_markov(TO_exact)
    #@test is_almost_markov(TO_exact_p)
    #@test is_markov(TO_exact)
    #@test is_markov(TO_exact_p)

    # Transfer operators from regular triangulations *not guaranteed to be invariant*
    TO = transferoperator_triang(triang)
    TO_approx = transferoperator_triang(triang, exact = false, parallel = false)
    TO_approx_rand = transferoperator_triang(triang, exact = false,
					parallel = false, sample_randomly = true)
    #TO_exact = transferoperator(triang, exact = true, parallel = false)
    #TO_exact_p = transferoperator(triang, exact = true, parallel = true)

    @test typeof(TO) == ApproxSimplexTransferOperator #default = approx
    @test typeof(TO_approx) == ApproxSimplexTransferOperator
    @test typeof(TO_approx_rand) == ApproxSimplexTransferOperator
    #@test typeof(TO_exact) == ExactSimplexTransferOperator
    #@test typeof(TO_exact_p) == ExactSimplexTransferOperator
    #@test all(TO_exact.transfermatrix .== TO_exact_p.transfermatrix)

end


text:
"""
    levene(observations::Array{Union{Number, Vector{Number}}})
    levene(observations::Vector{Number}, factorassignments::Vector{Vector{Any}})
    levene(df::DataFrame, observationscolumn::Symbol, factorcolumns::Vector{Symbol})

Performs Levene's homogeniety of variance test.

Operates on at least 3 crossed factors (fixed or random) and arbitrarily many random nested factors on balanced data. Levene's test is identical to performing an ANOVA on the absolute deviations of each cell.

If Levene's test is significant, the data is heteroscedastic (variances differ). If it is not, variances could differ or the sample size could be too small to tell.

# Arguments
- `observations`: Array containing the values to test. For the array, each dimension is a factor level, such that observations[2,5,3] indicates the 2nd level of the first factor, the 5th level of the second factor, and the 3rd level of the third factor. May contain values or vectors of values, where the vector contains replicates. Factors should be ordered with least significant first. For the vector, must provide `factorassignments` to specify factor levels.
- `factorassignments`: Vector of vectors of integers specifying how each observation is assigned to a factor level. Provide this when `observations` is given as a vector. Factor levels do not have to be consecutive or ordered. Nested factors must reuse factor levels currently.

Notes: Throws an exception if it can tell there are no replicates. Do not use with on any data structure you would pass into `anova()` using `hasrepliates=false`.

Output: `AnovaData` structure containing the test results for each factor.

# Examples
```julia
levene(observations)
```
"""
function levene(observations::AbstractArray{T}) where {T <: Union{Number, AbstractVector{<:Number}}}
    length(observations) > 0 || return

    firstlevelreplicates = eltype(observations) <: Number ? true : false
    nfactors = ndims(observations) - (firstlevelreplicates ? 1 : 0)

    nreplicates = firstlevelreplicates ? size(observations, 1) : length(observations[1])
    nreplicates > 1 || throw(ErrorException("Levene's test is invalid if there are no replicates."))
    firstlevelreplicates || all(c -> length(c) == nreplicates, observations) || throw(ErrorException("All cells must have the same number of replicates."))

    ncells = Int.(length(observations) / (firstlevelreplicates ? nreplicates : 1))
    nfactorlevels = firstlevelreplicates ? [size(observations)...][Not(1)] : [size(observations)...]

    observations = firstlevelreplicates ? reshape(observations, nreplicates, ncells) : reshape(observations, ncells)
    observations = upcat(observations)
    levenekernel(observations, nreplicates)
end

function levene(observations::AbstractVector{T}, factorassignments::AbstractVector{<:AbstractVector}) where {T <: Number}
    length(observations) > 0 || return
    nfactors = length(factorassignments)
    N = length(observations)
    N % nfactors == 0 || error("Design is unbalanced.")
    all(length.(factorassignments) .== N) || error("Each observation must have an assignment for each factor.")

    factorlevels = factorassignments .|> unique .|> sort
    nfactorlevels = length.(factorlevels)
    all(N .& nfactorlevels .== 0) || error("Design is unbalanced.")
    factorlevelcounts = [[count(l -> l == factorlevels[i][j], factorassignments[i]) for j ∈ 1:nfactorlevels[i]] for i ∈ 1:nfactors]
    nperfactorlevel = factorlevelcounts .|> unique
    all(nperfactorlevel .|> length .== 1) || error("Design is unbalanced.")
    nperfactorlevel = nperfactorlevel .|> first

    if !(isa(factorassignments, Number)) || any(maximum.(factorlevels) .> nfactorlevels)
        compressedfactorlevels = [1:i for i ∈ nfactorlevels]
        factorlevelremapping = [factorlevels[i] .=> compressedfactorlevels[i] for i ∈ 1:nfactors]
        factorassignments = [replace(factorassignments[i], factorlevelremapping[i]...) for i ∈ 1:nfactors]
    end

    #deduplicate above

    ncells = prod(nfactorlevels)
    nreplicates = Int(N / ncells)

    nlevels = [nreplicates; nfactorlevels]
    sortorder = sortperm(repeat(1:nreplicates, Int(N / nreplicates)) .+
                         sum([factorassignments[i] .* prod(nlevels[1:i]) for i ∈ 1:nfactors]))
    observationsmatrix = reshape(observations[sortorder], nreplicates, ncells)

    levenekernel(observationsmatrix, nreplicates)
end

function levenekernel(observations, nreplicates)
    cellmeans = meanfirstdim(observations)
    if eltype(observations) <: Number
        absdevs = abs.(observations .- cellmeans')
    else
        absdevs = [abs.(observations[i] .- cellmeans[i]) for i ∈ eachindex(observations)]
    end
    anovakernel(absdevs, ["Groups"], [fixed], nreplicates > 1)
end

export levene


text:
export LazyObservable, lazy, set!

import Base: show

# TODO: Untested and undocumented
lazy(::Type{T} = Any) where T = LazyObservable{T}(Vector{Tuple{LazySubscription, Any}}(), nothing)

# Lazy subscription represents a reference to an original subscription or a lazy observable
mutable struct LazySubscription <: Teardown
    subscription_or_lazy :: Any
end

@subscribable mutable struct LazyObservable{D} <: Subscribable{D}
    pending :: Vector{Tuple{LazySubscription, Any}}
    stream  :: Any
end

getpending(lazy::LazyObservable) = lazy.pending
pushpending!(lazy::LazyObservable, subscription::LazySubscription, actor) = push!(lazy.pending, (subscription, actor))

getstream(lazy::LazyObservable)          = lazy.stream
setstream!(lazy::LazyObservable, stream) = lazy.stream = stream

@inline set!(lazy::LazyObservable, observable::S) where S = on_lazy_set!(lazy, as_subscribable(S), observable)

@inline on_lazy_set!(lazy::LazyObservable{D},  ::InvalidSubscribableTrait,       observable) where D                = throw(InvalidSubscribableTraitUsageError(observable))
@inline on_lazy_set!(lazy::LazyObservable{D1}, ::SimpleSubscribableTrait{D2},    observable) where { D1, D2 <: D1 } = _on_lazy_set!(lazy, observable)
@inline on_lazy_set!(lazy::LazyObservable{D1}, ::ScheduledSubscribableTrait{D2}, observable) where { D1, D2 <: D1 } = _on_lazy_set!(lazy, observable)

@inline function _on_lazy_set!(lazy, observable)

    if getstream(lazy) !== nothing
        error("Lazy stream cannot be set twice")
    end

    setstream!(lazy, observable)

    foreach(lazy.pending) do (subscription, actor)
        subscription.subscription_or_lazy = subscribe!(observable, actor)
    end

    empty!(lazy.pending)
    resize!(lazy.pending, 0)

    return nothing
end

function on_subscribe!(observable::LazyObservable{D}, actor) where D
    stream = getstream(observable)
    if stream !== nothing
        return LazySubscription(subscribe!(stream, actor))
    else
        subscription = LazySubscription(observable)
        pushpending!(observable, subscription, actor)
        return subscription
    end
end

Base.show(io::IO, ::LazyObservable{D}) where D = print(io, "LazyObservable($D)")

as_teardown(::Type{ <: LazySubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::LazySubscription)
    return unsubscribe_lazy!(subscription, subscription.subscription_or_lazy)
end

function unsubscribe_lazy!(subscription::LazySubscription, observable::LazyObservable)
    filter!(p -> first(p) !== subscription, observable.pending)
    return nothing
end

function unsubscribe_lazy!(subscription::LazySubscription, teardown::Teardown)
    return unsubscribe!(teardown)
end