text:
facts("Reservoir data structures: reservoir") do
  n_weeks = 52
  n_scenarios = 5
  year = 1999
  envflow = .02m^3/s
  maxflow = 5.m^3/s

  name1 = "RivN1"
  name2 = "RivN2"
  name3 = "RivD3"
  name4 = "RivD4"

  nr1 = gen_river(name1, n_weeks, year, n_scenarios)
  nr2 = gen_river(name2, n_weeks, year, n_scenarios)
  dr1 = gen_river(name3, n_weeks, year, n_scenarios, envflow, maxflow)
  dr2 = gen_river(name4, n_weeks, year, n_scenarios, envflow, maxflow)

  nrs = [nr1, nr2]
  nrs_names = [name1, name2]
  drs = [dr1, dr2]
  drs_names = [name3, name4]
  rs = [nr1, nr2, dr1, dr2]
  rs_names = [name1, name2, name3, name4]

  r_name = "Reservoir: as simple as possible"
  r_min_capacity = 2.e8m^3
  r_max_capacity = 8.e8m^3

  const a = 1.
  const a_discharge = a * 1.m^3/s
  const a_power = a * 1.W
  const a_height = a * 1.m
  const b = 2.
  const b_discharge = b * 1.m^3/s
  const b_power = b * 1.W
  const c = 3.

  function reservoirTestShortVersions(r::Reservoir)
    outputIds = [:spillway, :penstock, :bottomOutlets]
    purposeIds = [:drinkingWater, :environmentalFlow]
    purposeTypes = [WaterWithdrawalPurpose, DeterministicWaterWithdrawalPurpose, DeterministicDrinkingWater,
                    DeterministicEnvironmentalFlow, DeterministicHydropower]

    @fact name(r) --> getName(r)
    @fact variant(r) --> getVariant(r)
    @fact minCapacity(r) --> getMinimumCapacity(r)
    @fact maxCapacity(r) --> getMaximumCapacity(r)
    @fact penstockHydropower(r) --> getPenstockHydropower(r)
    @fact hydropowerUnits(r) --> getHydropowerUnits(r)
    @fact outputs(r) --> getOutputs(r)
    @fact outputsId(r) --> getOutputIds(r)
    @fact bathymetry(r) --> getBathymetry(r)

    # Getter with an ID: return a Nullable (the other option: throw an error when there is nothing; hard to
    # ensure as a generic test).
    for s in outputIds; @fact output(r, s, nullable=true) --> exactly(getOutput(r, s, nullable=true)); end
    for s in purposeIds; @fact purpose(r, s, nullable=true) --> exactly(getPurpose(r, s, nullable=true)); end

    # Getter with a type: returns arrays, hence no exactly().
    for pt in purposeTypes; @fact purpose(r, pt) --> getPurpose(r, pt); end
  end

  function reservoirTestConsistency(r::Reservoir)
    @fact ReservoirManagement.countScenariosWithConsistency(r) --> countScenarios(r)
    @fact ReservoirManagement.countTimeStepsWithConsistency(r) --> countTimeSteps(r)
    @fact ReservoirManagement.getPeriodWithConsistency(r) --> getPeriod(r)
  end

  function reservoirTestRivers(r::Reservoir, all_rivers, all_natural_rivers, all_diverted_rivers,
                               all_rivers_names, all_natural_rivers_names, all_diverted_rivers_names)
    # Iterators.
    coll_all_rivers = collect(rivers(r))
    coll_all_natural_rivers = collect(naturalRivers(r))
    coll_all_diverted_rivers = collect(divertedRivers(r))
    for riv in all_rivers;          @fact in(riv, coll_all_rivers)          --> true; end
    for riv in all_natural_rivers;  @fact in(riv, coll_all_natural_rivers)  --> true; end
    for riv in all_diverted_rivers; @fact in(riv, coll_all_diverted_rivers) --> true; end

    # Nice iterators.
    each_all_rivers = collect(eachriver(r))
    each_all_natural_rivers = collect(eachnaturalriver(r))
    each_all_diverted_rivers = collect(eachdivertedriver(r))
    for riv in all_rivers;          @fact in(riv, coll_all_rivers)          --> true; end
    for riv in all_natural_rivers;  @fact in(riv, coll_all_natural_rivers)  --> true; end
    for riv in all_diverted_rivers; @fact in(riv, coll_all_diverted_rivers) --> true; end

    # Collectives.
    for riv in all_rivers;          @fact in(riv, getRivers(r))         --> true; end
    for riv in all_natural_rivers;  @fact in(riv, getNaturalRivers(r))  --> true; end
    for riv in all_diverted_rivers; @fact in(riv, getDivertedRivers(r)) --> true; end

    # Raw indexing.
    for i in 1:length(all_natural_rivers);  @fact in(getNaturalRiver(r, i),  all_natural_rivers)  --> true; end
    for i in 1:length(all_diverted_rivers); @fact in(getDivertedRiver(r, i), all_diverted_rivers) --> true; end

    # Name indexing.
    @fact getRiver(r, "Nice name for a river") --> isnull
    @fact hasRiver(r, "Nice name for a river") --> false
    for n in all_rivers_names; @fact getRiver(r, n) --> not(isnull); end
    for n in all_rivers_names; @fact hasRiver(r, n) --> true; end
    for n in all_rivers_names; @fact in(get(getRiver(r, n)), all_rivers) --> true; end

    # Name retrieving.
    for n in all_rivers_names;          @fact in(n, getRiverNames(r))         --> true; end
    for n in all_natural_rivers_names;  @fact in(n, getNaturalRiverNames(r))  --> true; end
    for n in all_diverted_rivers_names; @fact in(n, getDivertedRiverNames(r)) --> true; end
  end

  function reservoirTestScenarioIterator(r::Reservoir)
    for s in eachscenario(r)
      @fact countScenarios(s) --> 1
      for riv in rivers(s)
        actual = getScenario(riv, 1)
        expecteds = [getScenario(get(getRiver(r, getName(riv))), sp) for sp in 1:countScenarios(r)]
        @fact in(actual, expecteds) --> true
      end
    end
  end

  context("Error cases") do
    # No input: neither rivers nor extraction.
    @fact_throws Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity))

    # Duplicate IDs.
    p = DeterministicDrinkingWater(a_discharge)
    @fact_throws Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), purposes=[p, p])
  end

  context("Simplest case") do
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2])

    # Basic getters.
    @fact getName(r) --> r_name
    @fact getMinimumCapacity(r) --> ReservoirManagement._from_unitful(r_min_capacity)
    @fact getMaximumCapacity(r) --> ReservoirManagement._from_unitful(r_max_capacity)
    @fact getPenstockHydropower(r) --> 0.
    @fact getHydropowerUnits(r) --> []
    @fact getPurposes(r) --> []
    @fact getPurposeIds(r) --> []
    @fact getOutputs(r) --> []
    @fact getOutputIds(r) --> []
    @fact getBathymetry(r) --> isempty

    reservoirTestShortVersions(r)
    reservoirTestConsistency(r)
    reservoirTestRivers(r, rs, nrs, drs, rs_names, nrs_names, drs_names)
    reservoirTestScenarioIterator(r)

    # Find purposes and outputs.
    @fact_throws hasPurpose(r, Reservoir)
    @fact hasPurpose(r, WaterWithdrawalPurpose) --> false
    @fact hasPurpose(r, DeterministicWaterWithdrawalPurpose) --> false
    @fact hasPurpose(r, DeterministicDrinkingWater) --> false
    @fact hasPurpose(r, DeterministicEnvironmentalFlow) --> false
    @fact hasPurpose(r, DeterministicHydropower) --> false

    @fact_throws getPurpose(r, Reservoir)
    @fact getPurpose(r, WaterWithdrawalPurpose) --> isempty
    @fact getPurpose(r, DeterministicWaterWithdrawalPurpose) --> isempty
    @fact getPurpose(r, DeterministicDrinkingWater) --> isempty
    @fact getPurpose(r, DeterministicEnvironmentalFlow) --> isempty
    @fact getPurpose(r, DeterministicHydropower) --> isempty

    @fact hasPurpose(r, :drinkingWater) --> false
    @fact hasPurpose(r, :environmentalFlow) --> false
    @fact_throws getPurpose(r, :drinkingWater, nullable=false)
    @fact_throws getPurpose(r, :environmentalFlow, nullable=false)
    @fact getPurpose(r, :drinkingWater, nullable=true) --> isnull
    @fact getPurpose(r, :environmentalFlow, nullable=true) --> isnull

    @fact hasOutput(r, :spillway) --> false
    @fact hasOutput(r, :penstock) --> false
    @fact hasOutput(r, :bottomOutlets) --> false
    @fact_throws getOutput(r, :spillway, nullable=false)
    @fact_throws getOutput(r, :penstock, nullable=false)
    @fact_throws getOutput(r, :bottomOutlets, nullable=false)
    @fact getOutput(r, :spillway, nullable=true) --> isnull
    @fact getOutput(r, :penstock, nullable=true) --> isnull
    @fact getOutput(r, :bottomOutlets, nullable=true) --> isnull

    @fact copy(r) --> r

    # Counting things.
    @fact countScenarios(r) --> n_scenarios
    @fact countTimeSteps(r) --> n_weeks
    @fact getPeriod(r) --> Week(1)
  end

  context("Simplest case, agregating") do
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2])

    @fact_throws totalInflow(r) # Multiple scenarios.
    for s in 1:n_scenarios
      # Inflow as arrays.
      actual = totalInflow(r, s, asTimeSeries=false)
      expected = sum([getScenario(riv, s) for riv in rs])
        @fact actual --> expected

      # Inflow as time series.
      actual = totalInflow(r, s, asTimeSeries=true)
      nr1_sc = getScenarios(nr1)[1]
      expected = TimeArray(timestamp(nr1_sc), sum([getScenario(riv, s) for riv in rs]), colnames(nr1_sc), meta(nr1_sc))
      @fact timestamp(actual) --> timestamp(expected)
      @fact values(actual) --> values(expected)
      @fact colnames(actual) --> colnames(expected)
      @fact meta(actual) --> meta(expected)
    end
  end

  context("Simplest case, altering: deleting a scenario") do
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2])
    r2 = removeScenarios(r, IntSet(1))

    @fact getName(r2) --> r_name
    @fact getMinimumCapacity(r2) --> ReservoirManagement._from_unitful(r_min_capacity)
    @fact getMaximumCapacity(r2) --> ReservoirManagement._from_unitful(r_max_capacity)
    @fact getPenstockHydropower(r2) --> 0.
    @fact getHydropowerUnits(r2) --> []
    @fact getPurposes(r2) --> []
    @fact getPurposeIds(r2) --> []
    @fact getOutputs(r2) --> []
    @fact getOutputIds(r2) --> []
    @fact getBathymetry(r2) --> isempty

    reservoirTestShortVersions(r2)
    reservoirTestConsistency(r2)
    # Cannot apply reservoirTestRivers: the rivers have changed, they have fewer scenarios!
    reservoirTestScenarioIterator(r)

    @fact countScenarios(r2) --> n_scenarios - 1
    @fact countTimeSteps(r2) --> n_weeks
    @fact getPeriod(r2) --> Week(1)

    for rn in rs_names
      for s in 1:n_scenarios - 1
        @fact getScenario(get(getRiver(r2, rn)), s) --> getScenario(get(getRiver(r, rn)), s + 1)
      end
    end

    # Other APIs: don't directly give an IntSet.
    r3 = removeScenarios(r, 1)
    r4 = removeScenarios(r, [1])
    for rn in rs_names
      for s in 1:n_scenarios - 1
        @fact getScenario(get(getRiver(r3, rn)), s) --> getScenario(get(getRiver(r, rn)), s + 1)
        @fact getScenario(get(getRiver(r4, rn)), s) --> getScenario(get(getRiver(r, rn)), s + 1)
      end
    end
  end

  context("Simplest case, altering: deleting multiple scenarios") do
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2])
    r2 = removeScenarios(r, IntSet([1, 3, 5]))

    @fact getName(r2) --> r_name
    @fact getMinimumCapacity(r2) --> ReservoirManagement._from_unitful(r_min_capacity)
    @fact getMaximumCapacity(r2) --> ReservoirManagement._from_unitful(r_max_capacity)
    @fact getPenstockHydropower(r2) --> 0.
    @fact getHydropowerUnits(r2) --> []
    @fact getPurposes(r2) --> []
    @fact getPurposeIds(r2) --> []
    @fact getOutputs(r2) --> []
    @fact getOutputIds(r2) --> []
    @fact getBathymetry(r2) --> isempty

    reservoirTestShortVersions(r2)
    reservoirTestConsistency(r2)
    # Cannot apply reservoirTestRivers: the rivers have changed, they have fewer scenarios!
    reservoirTestScenarioIterator(r)

    @fact countScenarios(r2) --> n_scenarios - 3
    @fact countTimeSteps(r2) --> n_weeks
    @fact getPeriod(r2) --> Week(1)

    for rn in rs_names
      @fact getScenario(get(getRiver(r2, rn)), 1) --> getScenario(get(getRiver(r, rn)), 2)
      @fact getScenario(get(getRiver(r2, rn)), 2) --> getScenario(get(getRiver(r, rn)), 4)
    end

    # Other APIs: don't directly give an IntSet.
    r3 = removeScenarios(r, [1, 3, 5])
    for rn in rs_names
      @fact getScenario(get(getRiver(r3, rn)), 1) --> getScenario(get(getRiver(r, rn)), 2)
      @fact getScenario(get(getRiver(r3, rn)), 2) --> getScenario(get(getRiver(r, rn)), 4)
    end
  end

  context("Simplest case, altering: keeping one time step") do
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2])
    chosen_t_step = 8
    r2 = keepTimeSteps(r, chosen_t_step)

    @fact getName(r2) --> r_name
    @fact getMinimumCapacity(r2) --> ReservoirManagement._from_unitful(r_min_capacity)
    @fact getMaximumCapacity(r2) --> ReservoirManagement._from_unitful(r_max_capacity)
    @fact getPenstockHydropower(r2) --> 0.
    @fact getHydropowerUnits(r2) --> []
    @fact getPurposes(r2) --> []
    @fact getPurposeIds(r2) --> []
    @fact getOutputs(r2) --> []
    @fact getOutputIds(r2) --> []
    @fact getBathymetry(r2) --> isempty

    reservoirTestShortVersions(r2)
    reservoirTestConsistency(r2)
    # Cannot apply reservoirTestRivers: the rivers have changed, they have fewer time steps!
    reservoirTestScenarioIterator(r)

    @fact countScenarios(r2) --> n_scenarios
    @fact countTimeSteps(r2) --> 1
    @fact getPeriod(r2) --> Week(1)

    for rn in rs_names
      for s in 1:n_scenarios
        @fact getScenario(get(getRiver(r2, rn)), s, 1) --> getScenario(get(getRiver(r, rn)), s, chosen_t_step)
      end
    end
  end

  context("Simplest case, altering: keeping multiple time steps") do
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity), rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2])
    first_chosen_t_step = 8
    last_chosen_t_step = 42
    n_chosen_t_steps = length(collect(first_chosen_t_step:last_chosen_t_step))
    r2 = keepTimeSteps(r, first_chosen_t_step:last_chosen_t_step)

    @fact getName(r2) --> r_name
    @fact getMinimumCapacity(r2) --> ReservoirManagement._from_unitful(r_min_capacity)
    @fact getMaximumCapacity(r2) --> ReservoirManagement._from_unitful(r_max_capacity)
    @fact getPenstockHydropower(r2) --> 0.
    @fact getHydropowerUnits(r2) --> []
    @fact getPurposes(r2) --> []
    @fact getPurposeIds(r2) --> []
    @fact getOutputs(r2) --> []
    @fact getOutputIds(r2) --> []
    @fact getBathymetry(r2) --> isempty

    reservoirTestShortVersions(r2)
    reservoirTestConsistency(r2)
    # Cannot apply reservoirTestRivers: the rivers have changed, they have fewer time steps!
    reservoirTestScenarioIterator(r)

    @fact countScenarios(r2) --> n_scenarios
    @fact countTimeSteps(r2) --> n_chosen_t_steps
    @fact getPeriod(r2) --> Week(1)

    for rn in rs_names
      for s in 1:n_scenarios
        for t in 1:n_chosen_t_steps
          @fact getScenario(get(getRiver(r2, rn)), s, t) --> getScenario(get(getRiver(r, rn)), s, first_chosen_t_step + t - 1)
        end
      end
    end
  end

  context("Only hydropower") do
    hp = HydropowerUnit(max_discharge=a_discharge, max_power=b_power, efficiency=1./c)
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity),
                  rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2],
                  hydropower=[hp])

    # Basic getters.
    @fact getName(r) --> r_name
    @fact getMinimumCapacity(r) --> ReservoirManagement._from_unitful(r_min_capacity)
    @fact getMaximumCapacity(r) --> ReservoirManagement._from_unitful(r_max_capacity)
    @fact getPenstockHydropower(r) --> ReservoirManagement._from_unitful(a_discharge)
    @fact getHydropowerUnits(r) --> [hp]
    @fact getPurposes(r) --> [DeterministicHydropower(hp)]
    @fact getPurposeIds(r) --> [getPurposeId(hp)]
    @fact getOutputs(r) --> [HydropowerDamOutput(hp)]
    @fact getOutputIds(r) --> [getDamOutputId(hp)]
    @fact getBathymetry(r) --> isempty

    reservoirTestShortVersions(r)
    reservoirTestConsistency(r)
    reservoirTestRivers(r, rs, nrs, drs, rs_names, nrs_names, drs_names)
    reservoirTestScenarioIterator(r)

    # Find purposes and outputs.
    @fact_throws hasPurpose(r, Reservoir)
    @fact hasPurpose(r, WaterWithdrawalPurpose) --> true
    @fact hasPurpose(r, DeterministicWaterWithdrawalPurpose) --> true
    @fact hasPurpose(r, DeterministicDrinkingWater) --> false
    @fact hasPurpose(r, DeterministicEnvironmentalFlow) --> false
    @fact hasPurpose(r, DeterministicHydropower) --> true

    @fact_throws getPurpose(r, Reservoir)
    @fact getPurpose(r, WaterWithdrawalPurpose) --> [DeterministicHydropower(hp)]
    @fact getPurpose(r, DeterministicWaterWithdrawalPurpose) --> [DeterministicHydropower(hp)]
    @fact getPurpose(r, DeterministicDrinkingWater) --> isempty
    @fact getPurpose(r, DeterministicEnvironmentalFlow) --> isempty
    @fact getPurpose(r, DeterministicHydropower) --> [DeterministicHydropower(hp)]

    @fact hasPurpose(r, :drinkingWater) --> false
    @fact hasPurpose(r, :environmentalFlow) --> false
    @fact_throws getPurpose(r, :drinkingWater, nullable=false)
    @fact_throws getPurpose(r, :environmentalFlow, nullable=false)
    @fact getPurpose(r, :drinkingWater, nullable=true) --> isnull
    @fact getPurpose(r, :environmentalFlow, nullable=true) --> isnull

    @fact hasOutput(r, :spillway) --> false
    @fact hasOutput(r, :penstock) --> false
    @fact hasOutput(r, :bottomOutlets) --> false
    @fact_throws getOutput(r, :spillway, nullable=false)
    @fact_throws getOutput(r, :penstock, nullable=false)
    @fact_throws getOutput(r, :bottomOutlets, nullable=false)
    @fact getOutput(r, :spillway, nullable=true) --> isnull
    @fact getOutput(r, :penstock, nullable=true) --> isnull
    @fact getOutput(r, :bottomOutlets, nullable=true) --> isnull

    @fact hasOutput(r, getDamOutputId(hp)) --> true
    @fact getOutput(r, getDamOutputId(hp)) --> HydropowerDamOutput(hp)

    @fact copy(r) --> r

    # Counting things.
    @fact countScenarios(r) --> n_scenarios
    @fact countTimeSteps(r) --> n_weeks
    @fact getPeriod(r) --> Week(1)
  end

  context("Purposes and outputs") do
    p_dw = DeterministicDrinkingWater(a_discharge)
    p_ef = DeterministicEnvironmentalFlow(b_discharge)
    ps = Purpose[p_dw, p_ef]
    o_ep = ConstantDamOutput(:bottomOutlets, a_discharge)
    o_sw = ConstantDamOutput(:spillway, b_discharge)
    os = DamOutput[o_ep, o_sw]
    r = Reservoir(name=r_name, capacity=(r_min_capacity, r_max_capacity),
                  rivers_in=[nr1, nr2], rivers_diverted=[dr1, dr2],
                  purposes=ps, outputs=os)

    # Basic getters.
    @fact getName(r) --> r_name
    @fact getMinimumCapacity(r) --> ReservoirManagement._from_unitful(r_min_capacity)
    @fact getMaximumCapacity(r) --> ReservoirManagement._from_unitful(r_max_capacity)
    @fact getPenstockHydropower(r) --> 0.
    @fact getHydropowerUnits(r) --> isempty
    @fact getPurposes(r) --> ps
    @fact getPurposeIds(r) --> [:drinkingWater, :environmentalFlow]
    @fact getOutputs(r) --> os
    @fact getOutputIds(r) --> [:bottomOutlets, :spillway]
    @fact getBathymetry(r) --> isempty

    reservoirTestShortVersions(r)
    reservoirTestConsistency(r)
    reservoirTestRivers(r, rs, nrs, drs, rs_names, nrs_names, drs_names)
    reservoirTestScenarioIterator(r)

    # Find purposes and outputs.
    @fact_throws hasPurpose(r, Reservoir)
    @fact hasPurpose(r, WaterWithdrawalPurpose) --> true
    @fact hasPurpose(r, DeterministicWaterWithdrawalPurpose) --> true
    @fact hasPurpose(r, DeterministicDrinkingWater) --> true
    @fact hasPurpose(r, DeterministicEnvironmentalFlow) --> true
    @fact hasPurpose(r, DeterministicHydropower) --> false

    @fact_throws getPurpose(r, Reservoir)
    @fact getPurpose(r, WaterWithdrawalPurpose) --> ps
    @fact getPurpose(r, DeterministicWaterWithdrawalPurpose) --> ps
    @fact getPurpose(r, DeterministicDrinkingWater) --> Purpose[p_dw]
    @fact getPurpose(r, DeterministicEnvironmentalFlow) --> Purpose[p_ef]
    @fact getPurpose(r, DeterministicHydropower) --> isempty

    @fact hasPurpose(r, :drinkingWater) --> true
    @fact hasPurpose(r, :environmentalFlow) --> true
    @fact getPurpose(r, :drinkingWater, nullable=false) --> p_dw
    @fact getPurpose(r, :environmentalFlow, nullable=false) --> p_ef
    @fact getPurpose(r, :drinkingWater, nullable=true) --> exactly(Nullable(p_dw))
    @fact getPurpose(r, :environmentalFlow, nullable=true) --> exactly(Nullable(p_ef))

    @fact hasOutput(r, :spillway) --> true
    @fact hasOutput(r, :penstock) --> false
    @fact hasOutput(r, :bottomOutlets) --> true
    @fact getOutput(r, :spillway, nullable=false) --> o_sw
    @fact_throws getOutput(r, :penstock, nullable=false)
    @fact getOutput(r, :bottomOutlets, nullable=false) --> o_ep
    @fact getOutput(r, :spillway, nullable=true) --> exactly(Nullable(o_sw))
    @fact getOutput(r, :penstock, nullable=true) --> isnull
    @fact getOutput(r, :bottomOutlets, nullable=true) --> exactly(Nullable(o_ep))

    @fact copy(r) --> r

    # Counting things.
    @fact countScenarios(r) --> n_scenarios
    @fact countTimeSteps(r) --> n_weeks
    @fact getPeriod(r) --> Week(1)
  end
end


text:
# This file is a part of ArraysOfArrays.jl, licensed under the MIT License (MIT).

using ArraysOfArrays
using StatsBase
using Statistics
using Test

using UnsafeArrays
using Adapt

using ArraysOfArrays: full_consistency_checks, append_elemptr!, element_ptr

@testset "vector_of_arrays" begin
    ref_AoA1(T::Type, n::Integer) = n == 0 ? [Array{T}(undef, 5)][1:0] : [rand(T, rand(1:9)) for i in 1:n]
    ref_AoA2(T::Type, n::Integer) = n == 0 ? [Array{T}(undef, 4, 2)][1:0] : [rand(T, rand(1:4), rand(1:4)) for i in 1:n]
    ref_AoA3(T::Type, n::Integer) = n == 0 ? [Array{T}(undef, 3, 2, 4)][1:0] : [rand(T, rand(1:3), rand(1:3), rand(1:3)) for i in 1:n]

    @testset "element pointer handling" begin
        A = [2, 4, 5, 9]; B = [3, 6, 8]
        @test @inferred(append_elemptr!(deepcopy(A), B)) == [2, 4, 5, 9, 12, 14]
        @test @inferred(append_elemptr!([2], [5])) == [2]
    end

    @testset "ctors" begin
        A1 = ref_AoA1(Float32, 5)
        @test @inferred(VectorOfArrays(deepcopy(A1))) isa VectorOfArrays{Float32,1,0,Array{Float32,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfArrays(deepcopy(A1)) == A1
        @test @inferred(VectorOfArrays{Float64}(deepcopy(A1))) isa VectorOfArrays{Float64,1,0,Array{Float64,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfArrays{Float64}(deepcopy(A1)) == A1
        @test @inferred(VectorOfArrays{Float64,1}(deepcopy(A1))) isa VectorOfArrays{Float64,1,0,Array{Float64,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfArrays{Float64,1}(deepcopy(A1)) == A1

        @test @inferred(VectorOfVectors(deepcopy(A1))) isa VectorOfArrays{Float32,1,0,Array{Float32,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfVectors(deepcopy(A1)) == A1
        @test @inferred(VectorOfVectors{Float64}(deepcopy(A1))) isa VectorOfArrays{Float64,1,0,Array{Float64,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfVectors{Float64}(deepcopy(A1)) == A1

        A1_empty = ref_AoA1(Float32, 0)
        @test @inferred(VectorOfArrays(deepcopy(A1_empty))) isa VectorOfArrays{Float32,1,0,Array{Float32,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfArrays(deepcopy(A1_empty)) == A1_empty
        @test @inferred(VectorOfArrays{Float64}(deepcopy(A1_empty))) isa VectorOfArrays{Float64,1,0,Array{Float64,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfArrays{Float64}(deepcopy(A1_empty)) == A1_empty
        @test @inferred(VectorOfArrays{Float64,1}(deepcopy(A1_empty))) isa VectorOfArrays{Float64,1,0,Array{Float64,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfArrays{Float64,1}(deepcopy(A1_empty)) == A1_empty

        @test @inferred(VectorOfVectors(deepcopy(A1_empty))) isa VectorOfArrays{Float32,1,0,Array{Float32,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfVectors(deepcopy(A1_empty)) == A1_empty
        @test @inferred(VectorOfVectors{Float64}(deepcopy(A1_empty))) isa VectorOfArrays{Float64,1,0,Array{Float64,1},Array{Int,1},Array{Tuple{},1}}
        @test VectorOfVectors{Float64}(deepcopy(A1_empty)) == A1_empty

        A2 = ref_AoA2(Float32, 4)
        @test @inferred(VectorOfArrays(deepcopy(A2))) isa VectorOfArrays{Float32,2,1,Array{Float32,1},Array{Int,1},Array{Tuple{Int},1}}
        @test VectorOfArrays(deepcopy(A2)) == A2
        @test @inferred(VectorOfArrays{Float64}(deepcopy(A2))) isa VectorOfArrays{Float64,2,1,Array{Float64,1},Array{Int,1},Array{Tuple{Int},1}}
        @test VectorOfArrays{Float64}(deepcopy(A2)) == A2
        @test @inferred(VectorOfArrays{Float64,2}(deepcopy(A2))) isa VectorOfArrays{Float64,2,1,Array{Float64,1},Array{Int,1},Array{Tuple{Int},1}}
        @test VectorOfArrays{Float64,2}(deepcopy(A2)) == A2

        A3 = ref_AoA3(Float32, 3)
        @test @inferred(VectorOfArrays(deepcopy(A3))) isa VectorOfArrays{Float32,3,2,Array{Float32,1},Array{Int,1},Array{Tuple{Int,Int},1}}
        @test VectorOfArrays(deepcopy(A3)) == A3
        @test @inferred(VectorOfArrays{Float64}(deepcopy(A3))) isa VectorOfArrays{Float64,3,2,Array{Float64,1},Array{Int,1},Array{Tuple{Int,Int},1}}
        @test VectorOfArrays{Float64}(deepcopy(A3)) == A3
        @test @inferred(VectorOfArrays{Float64,3}(deepcopy(A3))) isa VectorOfArrays{Float64,3,2,Array{Float64,1},Array{Int,1},Array{Tuple{Int,Int},1}}
        @test VectorOfArrays{Float64,3}(deepcopy(A3)) == A3

        A3_empty = ref_AoA3(Float32, 0)
        @test @inferred(VectorOfArrays(deepcopy(A3_empty))) isa VectorOfArrays{Float32,3,2,Array{Float32,1},Array{Int,1},Array{Tuple{Int,Int},1}}
        @test VectorOfArrays(deepcopy(A3_empty)) == A3_empty
        @test @inferred(VectorOfArrays{Float64}(deepcopy(A3_empty))) isa VectorOfArrays{Float64,3,2,Array{Float64,1},Array{Int,1},Array{Tuple{Int,Int},1}}
        @test VectorOfArrays{Float64}(deepcopy(A3_empty)) == A3_empty
        @test @inferred(VectorOfArrays{Float64,3}(deepcopy(A3_empty))) isa VectorOfArrays{Float64,3,2,Array{Float64,1},Array{Int,1},Array{Tuple{Int,Int},1}}
        @test VectorOfArrays{Float64,3}(deepcopy(A3_empty)) == A3_empty
    end

    @testset "append! and vcat" begin
        A1 = ref_AoA3(Float32, 3); A2 = ref_AoA3(Float32, 0)
        A3 = ref_AoA3(Float32, 4); A4 = ref_AoA3(Float64, 2)

        B1 = VectorOfArrays(A1); B2 = VectorOfArrays(A2);
        B3 = VectorOfArrays(A3); B4 = VectorOfArrays(A4);

        @test @inferred(vcat(B1)) === B1
        full_consistency_checks(vcat(B1))

        @test @inferred(vcat(B1, B2)) isa VectorOfArrays
        @test vcat(B1, B2) == vcat(A1, A2)
        @test eltype(vcat(B1, B2)) == Array{Float32,3}
        full_consistency_checks(vcat(B1, B2))

        @test @inferred(vcat(B1, B3)) isa VectorOfArrays
        @test vcat(B1, B3) == vcat(A1, A3)
        @test eltype(vcat(B1, B3)) == Array{Float32,3}
        full_consistency_checks(vcat(B1, B3))

        @test @inferred(vcat(B1, B2, B3, B4)) isa VectorOfArrays
        @test vcat(B1, B2, B3, B4) == vcat(A1, A2, A3, A4)
        @test eltype(vcat(B1, B2, B3, B4)) == Array{Float64,3}
        full_consistency_checks(vcat(B1, B2, B3, B4))

        B1_copy = @inferred(copy(B1)); B3_copy = @inferred(copy(B3))
        append!(B1_copy, B3_copy)
        @test B1_copy.data == vcat(B1.data, B3.data)
    end

    @testset "indexing" begin
        V1 = @inferred(VectorOfArrays(ref_AoA3(Float32, 3)))
        V2 = @inferred(VectorOfArrays(ref_AoA3(Float32, 3)))
        V12 = vcat(V1, V2)
        ind_style = @inferred(IndexStyle(V12))

        @test ind_style == IndexLinear()

        for i in 1:length(V12)
            @test getindex(V12, i) == V12[i]
        end

        @test getindex(V12, 1:length(V12)) == V12

        @test @inferred(element_ptr(V12)) == V12.elem_ptr

        getindex_of_UR = @inferred(Base._getindex(ind_style, V12, 1:length(V12)))
        getindex_of_vector = @inferred(Base._getindex(ind_style, V12, collect(1:length(V12))))
        @test getindex_of_UR == V12
        @test getindex_of_vector == getindex_of_UR

        VV = @inferred(VectorOfVectors{Float64}())
        data = @inferred(rand(5))
        @inferred(push!(VV, data))
        @test @inferred(getindex(VV, 1)) == data
        @test @inferred(size(getindex(VV, 1))) == (5,)

## _view_reshape_spec not yet implemented ##
#       V1_copy = copy(V1)
#       V2_copy = copy(V2)
#       @test setindex!(V1_copy, V1, 1) == V1
#       setindex!(V2_copy, V2, 1)
#       @test V2_copy[1] == V2_copy[2]

## function mul(s) not yet implemented ##
#       sizehint!(v12_copy, 2, (2,2,3))

        zeroed_out = deepmap(x -> 0.0, V12)
        for i in zeroed_out
            @test @inferred(zeros(size(i))) == i
        end

        # Not a good test of variable depth?
        @test innermap(x -> 2*x, V12) == deepmap(x -> 2*x, V12)

    end

    @testset "indexing" begin
        V1 = @inferred(VectorOfArrays(ref_AoA3(Float32, 3)))
        V2 = @inferred(VectorOfArrays(ref_AoA3(Float32, 3)))
        V12 = vcat(V1, V2)
        ind_style = @inferred(IndexStyle(V12))
        @test ind_style == IndexLinear()
        for i in 1:length(V12)
            @test getindex(V12, i) == V12[i]
        end
        @test getindex(V12, 1:length(V12)) == V12

        @test @inferred(element_ptr(V12)) == V12.elem_ptr

## _view_reshape_spec not yet implemented ##
#       V1_copy = copy(V1)
#       V2_copy = copy(V2)
#       @test setindex!(V1_copy, V1, 1) == V1
#       setindex!(V2_copy, V2, 1)
#       @test V2_copy[1] == V2_copy[2]

## function mul(s) not yet implemented ##
#       sizehint!(v12_copy, 2, (2,2,3))

        # -------------------------------------------------------------------

        zeroed_out = deepmap(x -> 0.0, V12)
        for i in zeroed_out
            @test @inferred(zeros(size(i))) == i
        end

        # Not a good test of variable depth?
        @test innermap(x -> 2*x, V12) == deepmap(x -> 2*x, V12)

        f = x -> 0
        A = @inferred(AbstractArray{AbstractArray{Float64, 3},1}([rand(3,3,3), rand(3,3,3), rand(3,3,3), rand(3,3,3)]))
        A_inner = innermap(f,A)
        A_deep = deepmap(f, A)

        @test A_inner == A_deep

        @test @inferred(size(A_inner)) == @inferred(size(A_deep))
        @test innersize(A,1) == innersize(A,2)
        @test innersize(A,1) == innersize(A,3)
        for i in 1:length(A_deep)
            @test A_deep[i] == zeros(3,3,3)
        end
    end

    @testset "copy" begin
        A = ref_AoA3(Float32, 3);
        B = VectorOfArrays(A);

        @test typeof(@inferred copy(B)) == typeof(B)
        @test copy(A) == A
        @test copy(A) == B
    end

    @testset "empty" begin
        A = ref_AoA3(Float32, 3);
        B = VectorOfArrays(A);

        @test typeof(@inferred empty(B)) == typeof(B)
        @test empty(A) == empty(B)
    end

    @testset "adapt" begin
        A1 = VectorOfArrays(ref_AoA1(Float32, 3))
        @test @inferred(adapt(identity, A1)) == A1
        @test typeof(adapt(identity, A1)) == typeof(A1)

        A3 = VectorOfArrays(ref_AoA3(Float32, 3))
        @test @inferred(adapt(identity, A3)) == A3
        @test typeof(adapt(identity, A3)) == typeof(A3)
    end

    @testset "examples" begin
        VA = @inferred(VectorOfArrays{Float64, 2}())

        @inferred(push!(VA, rand(2, 3)))
        @inferred(push!(VA, rand(4, 2)))

        @test @inferred(size(VA[1]) == (2,3))
        @test @inferred(size(VA[2]) == (4,2))

        # -------------------------------------------------------------------

        VV = @inferred(VectorOfVectors{Float64}())
        d1 = @inferred(rand(5))
        d2 = @inferred(rand(4))

        @inferred(push!(VV, d1))
        @inferred(push!(VV, d2))

        @test VV[1] == d1
        @test VV[2] == d2

        @test @inferred(size(VV[1])) == (5,)
        @test @inferred(size(VV[2])) == (4,)

        # -------------------------------------------------------------------

        VA_flat = flatview(VA)
        @test VA_flat isa Vector{Float64}

        @test @inferred(uview(VA)) == VA
        @test @inferred(Base.unsafe_view(VA, 1:size(VA)[1])) == VA

        # -------------------------------------------------------------------

        VA_flat = flatview(VA)
        @test view(VA_flat, 7:14) == vec(VA[2])

        fill!(view(VA_flat, 7:14), 2.4)
        @test all(x -> x == 2.4, VA[2])

        fill!(view(VA_flat, 7:14), 4.2)
        @test all(x -> x == 4.2, VA[2])

        # -------------------------------------------------------------------

        @test length(@inferred resize!(VA, 1)) == 1
        @test_throws ArgumentError resize!(VA, 4)

        # -------------------------------------------------------------------

        using UnsafeArrays

        A = nestedview(rand(2,3,4,5), 2)

        @test isbits(A[2,2]) == false

        @test @uviews A begin
            isbits(A[2,2]) == true
        end

        # -------------------------------------------------------------------

        A = [1, 1, 2, 3, 3, 2, 2, 2]
        A_grouped_ref = [[1, 1], [2], [3, 3], [2, 2, 2]]
        elem_ptr = consgrouped_ptrs(A)
        elem_ptr32 = Int32.(consgrouped_ptrs(A))
        @test first.(@inferred(VectorOfVectors(A, elem_ptr))) == [1, 2, 3, 2]
        @test first.(@inferred(VectorOfVectors(A, elem_ptr32))) == [1, 2, 3, 2]

        B = [1, 2, 3, 4, 5, 6, 7, 8]
        B_grouped_ref = [[1, 2], [3], [4, 5], [6, 7, 8]]
        @test @inferred(VectorOfVectors(B, elem_ptr)) == B_grouped_ref
        @test @inferred(VectorOfVectors(B, elem_ptr32)) == B_grouped_ref

        C = [1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8]
        C_grouped_ref = [[1.1, 2.2], [3.3], [4.4, 5.5], [6.6, 7.7, 8.8]]

        @test @inferred(consgroupedview(A, B)) isa VectorOfVectors
        @test consgroupedview(A, B) == B_grouped_ref

        @test @inferred(consgroupedview(A, (B, C))) isa NTuple{2, VectorOfVectors}
        @test consgroupedview(A, (B, C)) == (B_grouped_ref, C_grouped_ref)

        data = (a = A, b = B, c = C)
        @test @inferred(consgroupedview(A, data)) isa NamedTuple{(:a, :b, :c),<:NTuple{3,AbstractVector}}
        result = consgroupedview(A, data)
        @test result == (a = A_grouped_ref, b = B_grouped_ref, c = C_grouped_ref)
        @test flatview(result.a) === data.a
        @test flatview(result.b) === data.b
        @test flatview(result.c) === data.c

        nestedV = @inferred(AbstractVector{AbstractArray{Float64, 4}}([rand(4,2,3,1), rand(5,3,1,3), rand(6,4,3,1), rand(9,2,1,2)]))
        VoA1 = @inferred(convert(VectorOfArrays, nestedV))
        @test @inferred(flatview(VoA1)) == VoA1.data
        VoA2 = @inferred(convert(VectorOfArrays{Float32, 4}, nestedV))
        @test @inferred(map(Float32, VoA1.data)) == VoA2.data

    end
end

ext: jl


text:
module StructJuMPwithMPI

using StructJuMP
using MPI

export mpi_wrapper

type MPIWrapper
    comm::MPI.Comm
    init::Function

    function MPIWrapper()
        instance = new(MPI.Comm(-1))
        finalizer(instance, freeMPIWrapper)

        instance.init = function(ucomm::MPI.Comm)
            if MPI.Initialized() && ucomm.val == -1
                instance.comm = MPI.COMM_WORLD
            elseif !MPI.Initialized()
                MPI.Init()
                instance.comm = MPI.COMM_WORLD
            elseif MPI.Initialized() && ucomm.val != -1
                instance.comm = ucomm
            elseif MPI.Finalized()
                error("MPI is already finalized!")
            else
                #doing nothing
            end
        end
        return instance
    end
end
function freeMPIWrapper(instance::MPIWrapper)
    if isdefined(:MPI) && MPI.Initialized() && !MPI.Finalized()
        MPI.Finalize()
    end
end

const mpi_wrapper = MPIWrapper()

function getMyRank()
    myrank = 0
    mysize = 1
    if isdefined(:MPI) && MPI.Initialized() && !MPI.Finalized()
        comm = MPI.COMM_WORLD
        mysize = MPI.Comm_size(comm)
        myrank = MPI.Comm_rank(comm)
    end
    return myrank, mysize
end

function StructJuMP.getProcIdxSet(mpi_wrapper::MPIWrapper, numScens::Integer)
    mysize = 1;
    myrank = 0;
    if isdefined(:MPI) == true && MPI.Initialized() == true
        comm = MPI.COMM_WORLD
        mysize = MPI.Comm_size(comm)
        myrank = MPI.Comm_rank(comm)
    end
    # Why don't we just take a round-and-robin?
    proc_idx_set = Int[];
    for s = myrank:mysize:(numScens-1)
        push!(proc_idx_set, s+1);
    end
    return proc_idx_set;
end

end # module


text:
using QuantumOpticsBase
using Test
using LinearAlgebra
import LinearAlgebra: mul!
using Random
import Base: ==, -

# Implement custom type with AbstractArray interface
mutable struct TestData{T,N,X} <: AbstractArray{T,N}
    x::X
    function TestData(x::X) where X
        x_ = convert(Matrix{ComplexF64}, x)
        new{ComplexF64,length(axes(x)),typeof(x_)}(x_)
    end
end
Base.size(A::TestData) = size(A.x)
Base.getindex(A::TestData, inds...) = getindex(A.x, inds...)
Base.setindex!(A::TestData, val, inds...) = setindex!(A.x, val, inds...)
Base.isapprox(A::TestData,B::TestData,args...;kwargs...) = isapprox(A.x,B.x,args...;kwargs...)

# Additional methods
Base.copy(A::TestData) = TestData(copy(A.x))
Base.kron(A::TestData,B::TestData) = TestData(kron(A.x,B.x))
# -(A::TestData) = TestData(-A.x)

# Auxiliary functions
function randtestoperator(args...)
    op_ = randoperator(args...)
    return Operator(op_.basis_l,op_.basis_r,TestData(copy(op_.data)))
end

D(a::DataOperator,b::DataOperator,tol=1e-14) = isapprox(a.data,b.data,atol=tol)
D(a::StateVector,b::StateVector,tol=1e-14) = isapprox(a.data,b.data,atol=tol)
D(a::AbstractOperator,b::AbstractOperator,tol=1e-12) = (tol > abs(tracedistance_nh(dense(a), dense(b))))

@testset "abstract-data" begin

Random.seed!(0)

b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(1)
b3b = GenericBasis(5)

b_l = b1a⊗b2a⊗b3a
b_r = b1b⊗b2b⊗b3b

# Test creation
@test_throws DimensionMismatch Operator(b1a, TestData([1 1 1; 1 1 1]))
@test_throws DimensionMismatch Operator(b1a, b1b, TestData([1 1; 1 1; 1 1]))
op1 = Operator(b1a, b1b, TestData([1 1 1; 1 1 1]))
op2 = Operator(b1b, b1a, TestData([1 1; 1 1; 1 1]))
op3 = op1*op2
@test op1 == dagger(op2)

# Test ' shorthand
@test dagger(op2) == op2'
@test transpose(op2) == conj(op2')

# Test copy
op1 = randtestoperator(b1a)
op2 = copy(op1)
@test op1.data == op2.data
@test !(op1.data === op2.data)
op2.data[1,1] = complex(10.)
@test op1.data[1,1] != op2.data[1,1]

# Arithmetic operations
# =====================
op_zero = DenseOperator(b_l, b_r)
op1 = randtestoperator(b_l, b_r)
op2 = randtestoperator(b_l, b_r)
op3 = randtestoperator(b_l, b_r)

x1 = Ket(b_r, rand(ComplexF64, length(b_r)))
x2 = Ket(b_r, rand(ComplexF64, length(b_r)))

xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))

# Addition
@test_throws DimensionMismatch op1 + dagger(op2)
@test D(op1 + op_zero, op1)
@test D(op1 + op2, op2 + op1)
@test D(op1 + (op2 + op3), (op1 + op2) + op3)

# Subtraction
@test_throws DimensionMismatch op1 - dagger(op2)
@test D(op1-op_zero, op1)
@test D(op1-op2, op1 + (-op2))
@test D(op1-op2, op1 + (-1*op2))
@test D(op1-op2-op3, op1-(op2+op3))

# Test multiplication
@test_throws DimensionMismatch op1*op2
@test D(op1*(x1 + 0.3*x2), op1*x1 + 0.3*op1*x2, 1e-12)
@test D((op1+op2)*(x1+0.3*x2), op1*x1 + 0.3*op1*x2 + op2*x1 + 0.3*op2*x2, 1e-12)

@test D((xbra1+0.3*xbra2)*op1, xbra1*op1 + 0.3*xbra2*op1)
@test D((xbra1+0.3*xbra2)*(op1+op2), xbra1*op1 + 0.3*xbra2*op1 + xbra1*op2 + 0.3*xbra2*op2)

@test D(op1*dagger(0.3*op2), 0.3*dagger(op2*dagger(op1)))
@test D((op1 + op2)*dagger(0.3*op3), 0.3*op1*dagger(op3) + 0.3*op2*dagger(op3), 1e-12)
@test D(0.3*(op1*dagger(op2)), op1*(0.3*dagger(op2)))

tmp = copy(op1)
conj!(tmp)
@test tmp == conj(op1) && conj(tmp.data) == op1.data

# Internal layout
b1 = GenericBasis(2)
b2 = GenericBasis(3)
b3 = GenericBasis(4)
op1 = randtestoperator(b1, b2)
op2 = randtestoperator(b2, b3)
x1 = randstate(b2)
d1 = op1.data
d2 = op2.data
v = x1.data
@test (op1*x1).data ≈ [d1[1,1]*v[1] + d1[1,2]*v[2] + d1[1,3]*v[3], d1[2,1]*v[1] + d1[2,2]*v[2] + d1[2,3]*v[3]]
@test (op1*op2).data[2,3] ≈ d1[2,1]*d2[1,3] + d1[2,2]*d2[2,3] + d1[2,3]*d2[3,3]

# Test division
@test D(op1/7, (1/7)*op1)

# Tensor product
# ==============
op1a = randtestoperator(b1a, b1b)
op1b = randtestoperator(b1a, b1b)
op2a = randtestoperator(b2a, b2b)
op2b = randtestoperator(b2a, b2b)
op3a = randtestoperator(b3a, b3b)
op123 = op1a ⊗ op2a ⊗ op3a
@test isa(op123.data,TestData)
@test op123.basis_l == b_l
@test op123.basis_r == b_r

# Associativity
@test D((op1a ⊗ op2a) ⊗ op3a, op1a ⊗ (op2a ⊗ op3a))

# Linearity
@test D(op1a ⊗ (0.3*op2a), 0.3*(op1a ⊗ op2a))
@test D((0.3*op1a) ⊗ op2a, 0.3*(op1a ⊗ op2a))

# Distributivity
@test D(op1a ⊗ (op2a + op2b), op1a ⊗ op2a + op1a ⊗ op2b)
@test D((op2a + op2b) ⊗ op3a, op2a ⊗ op3a + op2b ⊗ op3a)

# Mixed-product property
@test D((op1a ⊗ op2a) * dagger(op1b ⊗ op2b), (op1a*dagger(op1b)) ⊗ (op2a*dagger(op2b)))

# Transpose
@test D(dagger(op1a ⊗ op2a), dagger(op1a) ⊗ dagger(op2a))
@test D(dagger(op1a ⊗ op2a), dagger(op1a) ⊗ dagger(op2a))

# Internal layout
a = Ket(b1a, rand(ComplexF64, length(b1a)))
b = Ket(b2b, rand(ComplexF64, length(b2b)))
ab = a ⊗ dagger(b)
@test ab.data[2,3] == a.data[2]*conj(b.data[3])
@test ab.data[2,1] == a.data[2]*conj(b.data[1])

shape = tuple(op123.basis_l.shape..., op123.basis_r.shape...)
idx = LinearIndices(shape)[2, 1, 1, 3, 4, 5]
@test op123.data[idx] == op1a.data[2,3]*op2a.data[1,4]*op3a.data[1,5]
@test reshape(op123.data, shape...)[2, 1, 1, 3, 4, 5] == op1a.data[2,3]*op2a.data[1,4]*op3a.data[1,5]

idx = LinearIndices(shape)[2, 1, 1, 1, 3, 4]
@test op123.data[idx] == op1a.data[2,1]*op2a.data[1,3]*op3a.data[1,4]
@test reshape(op123.data, shape...)[2, 1, 1, 1, 3, 4] == op1a.data[2,1]*op2a.data[1,3]*op3a.data[1,4]

# Test tr and normalize
op = Operator(GenericBasis(3), TestData([1 3 2;5 2 2;-1 2 5]))
@test 8 == tr(op)
op_normalized = normalize(op)
@test 8 == tr(op)
@test 1 == tr(op_normalized)
op_copy = deepcopy(op)
normalize!(op_copy)
@test tr(op) != tr(op_copy)
@test 1 ≈ tr(op_copy)
@test op === normalize!(op)

# Test partial tr of operators
b1 = GenericBasis(3)
b2 = GenericBasis(5)
b3 = GenericBasis(7)
op1 = randtestoperator(b1)
op2 = randtestoperator(b2)
op3 = randtestoperator(b3)
op123 = op1 ⊗ op2 ⊗ op3

@test D(op1⊗op2*tr(op3), ptrace(op123, 3))
@test D(op1⊗op3*tr(op2), ptrace(op123, 2))
@test D(op2⊗op3*tr(op1), ptrace(op123, 1))

@test D(op1*tr(op2)*tr(op3), ptrace(op123, [2,3]))
@test D(op2*tr(op1)*tr(op3), ptrace(op123, [1,3]))
@test D(op3*tr(op1)*tr(op2), ptrace(op123, [1,2]))

@test_throws ArgumentError ptrace(op123, [1,2,3])
x = randtestoperator(b1, b1⊗b2)
@test_throws ArgumentError ptrace(x, [1])
x = randtestoperator(b1⊗b1⊗b2, b1⊗b2)
@test_throws ArgumentError ptrace(x, [1, 2])
x = randtestoperator(b1⊗b2)
@test_throws ArgumentError ptrace(x, [1, 2])
x = randtestoperator(b1⊗b2, b2⊗b1)
@test_throws ArgumentError ptrace(x, [1])

op1 = randtestoperator(b1, b2)
op2 = randtestoperator(b3)

@test D(op1*tr(op2), ptrace(op1⊗op2, 2))

# Test expect
b1 = GenericBasis(3)
b2 = GenericBasis(5)
b3 = GenericBasis(7)
op1 = randtestoperator(b1)
op2 = randtestoperator(b2)
op3 = randtestoperator(b3)
op123 = op1 ⊗ op2 ⊗ op3
b_l = b1 ⊗ b2 ⊗ b3

state = randstate(b_l)
@test expect(op123, state) ≈ dagger(state)*op123*state
@test expect(1, op1, state) ≈ expect(op1, ptrace(state, [2, 3]))
@test expect(2, op2, state) ≈ expect(op2, ptrace(state, [1, 3]))
@test expect(3, op3, state) ≈ expect(op3, ptrace(state, [1, 2]))

state = randtestoperator(b_l)
@test expect(op123, state) ≈ tr(op123*state)
@test expect(1, op1, state) ≈ expect(op1, ptrace(state, [2, 3]))
@test expect(2, op2, state) ≈ expect(op2, ptrace(state, [1, 3]))
@test expect(3, op3, state) ≈ expect(op3, ptrace(state, [1, 2]))

# Permute systems
op1 = randtestoperator(b1a, b1b)
op2 = randtestoperator(b2a, b2b)
op3 = randtestoperator(b3a, b3b)
op123 = op1⊗op2⊗op3

op132 = op1⊗op3⊗op2
@test D(permutesystems(op123, [1, 3, 2]), op132)

op213 = op2⊗op1⊗op3
@test D(permutesystems(op123, [2, 1, 3]), op213)

op231 = op2⊗op3⊗op1
@test D(permutesystems(op123, [2, 3, 1]), op231)

op312 = op3⊗op1⊗op2
@test D(permutesystems(op123, [3, 1, 2]), op312)

op321 = op3⊗op2⊗op1
@test D(permutesystems(op123, [3, 2, 1]), op321)

# Test projector
xket = normalize(Ket(b_l, rand(ComplexF64, length(b_l))))
yket = normalize(Ket(b_l, rand(ComplexF64, length(b_l))))
xbra = dagger(xket)
ybra = dagger(yket)

@test D(projector(xket)*xket, xket)
@test D(xbra*projector(xket), xbra)
@test D(projector(xbra)*xket, xket)
@test D(xbra*projector(xbra), xbra)
@test D(ybra*projector(yket, xbra), xbra)
@test D(projector(yket, xbra)*xket, yket)

# Test gemv
b1 = GenericBasis(3)
b2 = GenericBasis(5)
op = randtestoperator(b1, b2)
xket = randstate(b2)
xbra = dagger(randstate(b1))
rket = randstate(b1)
rbra = dagger(randstate(b2))
alpha = complex(0.7, 1.5)
beta = complex(0.3, 2.1)
@test isa(op.data, TestData)

rket_ = deepcopy(rket)
QuantumOpticsBase.mul!(rket_,op,xket,complex(1.0),complex(0.))
@test D(rket_, op*xket)

@test isa(op.data, TestData)
rket_ = deepcopy(rket)
QuantumOpticsBase.mul!(rket_,op,xket,alpha,beta)
@test D(rket_, alpha*op*xket + beta*rket)

rbra_ = deepcopy(rbra)
QuantumOpticsBase.mul!(rbra_,xbra,op,complex(1.0),complex(0.))
@test D(rbra_, xbra*op)

rbra_ = deepcopy(rbra)
QuantumOpticsBase.mul!(rbra_,xbra,op,alpha,beta)
@test D(rbra_, alpha*xbra*op + beta*rbra)

# Test gemm
b1 = GenericBasis(37)
b2 = GenericBasis(53)
b3 = GenericBasis(41)
op1 = randtestoperator(b1, b2)
op2 = randtestoperator(b2, b3)
r = randtestoperator(b1, b3)
alpha = complex(0.7, 1.5)
beta = complex(0.3, 2.1)

r_ = deepcopy(r)
QuantumOpticsBase.mul!(r_,op1,op2,complex(1.0),complex(0.))
@test D(r_, op1*op2)

r_ = deepcopy(r)
QuantumOpticsBase.mul!(r_,op1,op2,alpha,beta)
@test D(r_, alpha*op1*op2 + beta*r, 1e-12)

# Test Hermitian
bspin = SpinBasis(1//2)
bnlevel = NLevelBasis(2)
@test ishermitian(Operator(bspin, bspin, TestData([1.0 im; -im 2.0]))) == true
@test ishermitian(Operator(bspin, bnlevel, TestData([1.0 im; -im 2.0]))) == false

# Test broadcasting
op1_ = copy(op1)
@test isa(op1_.data, TestData)
op1 .= 2*op1
@test op1 == op1_ .+ op1_
op1 .= op1_
@test op1 == op1_
op1 .= op1_ .+ 3 * op1_
@test op1 == 4*op1_
@test_throws DimensionMismatch op1 .= op2
bf = FockBasis(3)
op3 = randtestoperator(bf)
@test_throws QuantumOpticsBase.IncompatibleBases op1 .+ op3
@test_throws ErrorException cos.(op1)

####################
# Test lazy tensor #
b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(6)
b3b = GenericBasis(5)

b_l = b1a⊗b2a⊗b3a
b_r = b1b⊗b2b⊗b3b

op1 = randtestoperator(b1a, b1b)
op2 = randtestoperator(b2a, b2b)
op3 = randtestoperator(b3a, b3b)

# Test creation
@test_throws AssertionError LazyTensor(b_l, b_r, [1], (randtestoperator(b1a),))
@test_throws AssertionError LazyTensor(b_l, b_r, [1, 2], (op1,))
@test_throws AssertionError LazyTensor(b_l, b_r, [1, 2], (op1, sparse(randtestoperator(b_l, b_l))))
@test_throws AssertionError LazyTensor(b_l, b_r, [1, 2], (randtestoperator(b_r, b_r), sparse(op2)))

# @test LazyTensor(b_l, b_r, [2, 1], [op2, op1]) == LazyTensor(b_l, b_r, [1, 2], [op1, op2])
x = randtestoperator(b2a)
@test LazyTensor(b_l, 2, x) == LazyTensor(b_l, b_l, [2], (x,))

# Test copy
x = 2*LazyTensor(b_l, b_r, [1,2], (randtestoperator(b1a, b1b), sparse(randtestoperator(b2a, b2b))))
x_ = copy(x)
@test x == x_
@test !(x === x_)
x_.operators[1].data[1,1] = complex(10.)
@test x.operators[1].data[1,1] != x_.operators[1].data[1,1]
x_.factor = 3.
@test x_.factor != x.factor
x_.indices[2] = 100
@test x_.indices != x.indices

# Test dense & sparse
I2 = identityoperator(b2a, b2b)
x = LazyTensor(b_l, b_r, [1, 3], (op1, sparse(op3)), 0.3)
@test D(0.3*op1⊗dense(I2)⊗op3, dense(x))
@test D(0.3*sparse(op1)⊗I2⊗sparse(op3), sparse(x))

# Test suboperators
@test QuantumOpticsBase.suboperator(x, 1) == op1
@test QuantumOpticsBase.suboperator(x, 3) == sparse(op3)
@test QuantumOpticsBase.suboperators(x, [1, 3]) == [op1, sparse(op3)]

# Arithmetic operations
subop1 = randtestoperator(b1a, b1b)
subop2 = randtestoperator(b2a, b2b)
subop3 = randtestoperator(b3a, b3b)
I1 = dense(identityoperator(b1a, b1b))
I2 = dense(identityoperator(b2a, b2b))
I3 = dense(identityoperator(b3a, b3b))
op1 = LazyTensor(b_l, b_r, [1, 3], (subop1, sparse(subop3)), 0.1)
op1_ = 0.1*subop1 ⊗ I2 ⊗ subop3
op2 = LazyTensor(b_l, b_r, [2, 3], (sparse(subop2), subop3), 0.7)
op2_ = 0.7*I1 ⊗ subop2 ⊗ subop3
op3 = 0.3*LazyTensor(b_l, b_r, 3, subop3)
op3_ = 0.3*I1 ⊗ I2 ⊗ subop3

x1 = Ket(b_r, rand(ComplexF64, length(b_r)))
x2 = Ket(b_r, rand(ComplexF64, length(b_r)))
xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))

# Addition
@test_throws ArgumentError op1 + op2
@test_throws ArgumentError op1 - op2
@test D(-op1_, -op1, 1e-12)

# Test multiplication
@test_throws DimensionMismatch op1*op2
@test D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2))
@test D((xbra1 + 0.3*xbra2)*op1, (xbra1 + 0.3*xbra2)*op1_)
@test D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2)
@test D(dagger(x1)*dagger(0.3*op2), dagger(x1)*dagger(0.3*op2_))
@test D(op1_*dagger(0.3*op2), op1_*dagger(0.3*op2_))
@test D(dagger(0.3*op2)*op1_, dagger(0.3*op2_)*op1_)
@test D(dagger(0.3*op2)*op1, dagger(0.3*op2_)*op1_)

#####################
# Test lazy product #
op1a = randtestoperator(b_l, b_r)
op1b = randtestoperator(b_r, b_l)
op2a = randtestoperator(b_l, b_r)
op2b = randtestoperator(b_r, b_l)
op3a = randtestoperator(b_l, b_l)
op1 = LazyProduct([op1a, sparse(op1b)])*0.1
op1_ = 0.1*(op1a*op1b)
op2 = LazyProduct([sparse(op2a), op2b], 0.3)
op2_ = 0.3*(op2a*op2b)
op3 = LazyProduct(op3a)
op3_ = op3a

x1 = Ket(b_l, rand(ComplexF64, length(b_l)))
x2 = Ket(b_l, rand(ComplexF64, length(b_l)))
xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))

# Addition
@test_throws ArgumentError op1 + op2
@test D(-op1_, -op1)

# Test multiplication
@test_throws DimensionMismatch op1a*op1a
@test D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2), 1e-11)
@test D((xbra1 + 0.3*xbra2)*op1, (xbra1 + 0.3*xbra2)*op1_, 1e-11)
@test D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2, 1e-11)
@test D(dagger(x1)*dagger(0.3*op2), dagger(x1)*dagger(0.3*op2_), 1e-11)
@test 0.3*LazyProduct(op1, sparse(op2)) == LazyProduct([op1, sparse(op2)], 0.3)
@test 0.3*LazyProduct(op1)*LazyProduct(sparse(op2)) == LazyProduct([op1, sparse(op2)], 0.3)

# Test gemv
op1 = randtestoperator(b_l, b_r)
op2 = randtestoperator(b_r, b_l)
op3 = randtestoperator(b_l, b_r)
op = LazyProduct([op1, sparse(op2), op3], 0.2)
op_ = 0.2*op1*op2*op3
tmp = 0.2*prod(dense.([op1*op2*op3]))

state = Ket(b_r, rand(ComplexF64, length(b_r)))
result_ = Ket(b_l, rand(ComplexF64, length(b_l)))
result = copy(result_)
QuantumOpticsBase.mul!(result,op,state,complex(1.),complex(0.))
@test D(result, op_*state, 1e-11)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op,state,alpha,beta)
@test D(result, alpha*op_*state + beta*result_, 1e-9)

state = Bra(b_l, rand(ComplexF64, length(b_l)))
result_ = Bra(b_r, rand(ComplexF64, length(b_r)))
result = copy(result_)
QuantumOpticsBase.mul!(result,state,op,complex(1.),complex(0.))
@test D(result, state*op_, 1e-9)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op,alpha,beta)
@test D(result, alpha*state*op_ + beta*result_, 1e-9)

# Test gemm
op1 = randtestoperator(b_l, b_r)
op2 = randtestoperator(b_r, b_l)
op3 = randtestoperator(b_l, b_r)
op = LazyProduct([op1, sparse(op2), op3], 0.2)
op_ = 0.2*op1*op2*op3

state = randtestoperator(b_r, b_r)
result_ = randtestoperator(b_l, b_r)
result = copy(result_)
result.data[1] = 1
result.data[1] != result_.data[1] || error("")
QuantumOpticsBase.mul!(result,op,state,complex(1.),complex(0.))
@test D(result, op_*state, 1e-9)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op,state,alpha,beta)
@test D(result, alpha*op_*state + beta*result_, 1e-9)

state = randtestoperator(b_l, b_l)
result_ = randtestoperator(b_l, b_r)
result = copy(result_)
QuantumOpticsBase.mul!(result,state,op,complex(1.),complex(0.))
@test D(result, state*op_, 1e-9)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op,alpha,beta)
@test D(result, alpha*state*op_ + beta*result_, 1e-9)

end # testset


text:
"""
    lanczos(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, numepoch::Number=10, perm::Bool=false, cper::Number=1f0)

Lanczos method.

Input Arguments
---------
- `input` : Julia Binary file generated by `OnlinePCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `rowmeanlist` : The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `rowvarlist` : The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `colsumlist` : The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `dim` : The number of dimension of PCA.
- `perm` : Whether the data matrix is shuffled at random.
- `cper` : Count per X (e.g. CPM: Count per million <1e+6>)

Output Arguments
---------
- `W` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `V` : Loading vectors of covariance matrix (No. rows of the data matrix × dim)
- `Scores` : Principal component scores
- `ExpVar` : Explained variance by the eigenvectors
- `TotalVar` : Total variance of the data matrix
- stop : Whether the calculation is converged
"""
function lanczos(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, numepoch::Number=10, perm::Bool=false, cper::Number=1f0)
    # Initial Setting
    pca = LANCZOS()
    N, M = nm(input)
    pseudocount, W, X, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar = init(input, pseudocount, dim, rowmeanlist, rowvarlist, colsumlist, pca, numepoch, cper, scale)
    T = zeros(Float32, numepoch, numepoch)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    # If not 0 the calculation is converged
    stop = 0
    k = 1
    n = 1
    s = 0
    # Each epoch s
    progress = Progress(numepoch)
    while(stop == 0 && k <= numepoch)
        v = zeros(Float32, M)
        open(input) do file
            stream = ZstdDecompressorStream(file)
            read!(stream, tmpN)
            read!(stream, tmpM)
            # Each step n
            while(stop == 0 && n <= N)
                next!(progress)
                # Row vector of data matrix
                read!(stream, x)
                normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
                if perm
                    normx .= normx[randperm(length(normx))]
                end
                # Power iteration
                v .+= normx * (normx' * W[:,k])
                n += 1
            end
            close(stream)
        end
        # NaN
        checkNaN(W, GD())
        # Check Float32
        @assert W[1,1] isa Float32
        # Ortho-normalization in Lanczos process
        T[k,k] = W[:,k]' * v
        if k == 1
            v = v .- T[k,k] * W[:,k]
        else
            v = v .- s * W[:,k-1] .- T[k,k] * W[:,k]
        end
        s = norm(v)
        if k != numepoch
            T[k+1,k] = s
            T[k,k+1] = s
        end
        W[:,k+1] = v / s
        k += 1
        if n == N + 1
            n = 1
        end
    end

    # Return, W, λ, V
    Y = eigvecs(T[1:numepoch, :])
    X = (W[:,1:numepoch] * Y[:,1:dim])
    out = WλV(X, input, dim, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, TotalVar, cper)
    out = (out[1], out[2], out[3], out[4], out[5], out[6], stop)
    if outdir isa String
        output(outdir, out, expvar)
    end
    return out
end


text:
# Build a decoding graph from a set of models, a lexicon and a language
# model (LM).
#

using ArgParse
using BSON
using MarkovModels

function getargs()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--add-sil", "-s"
            help = "add silence add the beginning/end of decoding graph (used only if lm is not specified)"
            action = :store_true
        "--disambig-sym", "-d"
            help = "Disambiguation symbol to distinguish word state from phone state"
            default = "wrd"
        "--lm", "-l"
            help = "Language Model graph to use (if not provided, build a flat LM)"
            default = ""
        "lexicon"
            help = "mapping word -> pronunciation"
            arg_type = String
            required = true
        "hmms"
            help = "unit hmms"
            arg_type = String
            required = true
        "decgraph"
            help = "output decoding graph"
            arg_type = String
            required = true
    end
    args = parse_args(s)
end

const args = getargs()
const dis_sym = args["disambig-sym"]
const lmfile = args["lm"]
const lexiconfile = args["lexicon"]
const hmmsfile = args["hmms"]
const out = args["decgraph"]

run(`mkdir -p $(dirname(out))`)

@info "loading the hmms"
const hmms = BSON.load(hmmsfile)[:hmms]

@info "loading the pronunciations"
const pronuns = Dict{String, Vector{FSM}}()
open(lexiconfile, "r") do f
    for line in eachline(f)
        tokens = split(line)
        word = "$dis_sym:$(tokens[1])"
        pronun = tokens[2:end]
        ps = get(pronuns, word, FSM[])
        push!(ps, LinearFSM(pronun))
        pronuns[word] = ps
    end
end

@info "compiling pronunciation graphs"
const lexicon = Dict{String, FSM}()
for (k, v) in pronuns
    lexicon[k] = union(v...) |> minimize!
end

function build_flat_lm(lexicon, addsil::Bool)
    lm = FSM()
    init = initstate(lm)
    final = finalstate(lm)
    if addsil
        init = addstate!(lm, label = "sil")
        link!(lm, initstate(lm), init)
        final = addstate!(lm, label = "sil")
        link!(lm, final, finalstate(lm))
    end
    anchor = addstate!(lm)
    link!(lm, init, anchor)
    link!(lm, anchor, final)
    for w in keys(lexicon)
        s = addstate!(lm, label = w)
        link!(lm, anchor, s)
        link!(lm, s, anchor)
    end
    lm |> weightnormalize!
end

const lm = begin
    if lmfile == ""
        @info "no LM provided, building a flat LM"
        build_flat_lm(lexicon, args["add-sil"])
    else
        @info "using a flat LM"
        lm = BSON.load(lmfile)[:lm]
    end
end

@info "compiling the decoding graph"
decodegraph = compose!(compose!(lm, lexicon), hmms) |> removenilstates!
bson(out, Dict(:decodegraph => decodegraph))

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
## Using HomotopyContinuation.jl for bifurcation diagrams

# function starting_solution(system, ksvalues; xs=HC.variables(system), start_system=:total_degree)
#     @debug "Solving system with :total_degree"
#     result = HC.solve(system(xs, ksvalues);
#                       start_system=start_system,
#                       show_progress=false,
#                       )
#     # @debug "Solving system with :polyhedral"
#     # result = HomotopyContinuation.solve(system(xs, ksvalues), show_progress=false)
#     return realpositve_solutions(result)
# end

# function track_solution(system, solution, start, target)
#     result = HC.solve(system, solution; start_parameters=start, target_parameters=target, show_progress=false)
#     return realpositve_solutions(result)
# end

function realpositve_solutions(result;
                               positive_tol = -1e-6,
                               only_nonsingular = false,
                               kwargs...)
    solutions = HC.real_solutions(result; only_nonsingular=only_nonsingular, kwargs...)
    return filter(s -> minimum(s) > positive_tol, solutions)
end

function solve_all!(out, hc_system, K_span;
                    variables = HC.variables(hc_system),
                    start_system = :total_degree,
                    show_progress = false,
                    kwargs...)
    for j in 1:length(K_span)
        @debug "Iteration $j, point and Ks $(Ks[j])"
        result = HC.solve(hc_system(variables, [K_span[j]]);
                          start_system = start_system,
                          show_progress = show_progress,
                          kwargs...)
        # @debug "Solving system with :polyhedral"
        # result = HomotopyContinuation.solve(system(xs, ksvalues), show_progress=false)
        out[j] = realpositve_solutions(result)
        # out[j] = starting_solution(hc_sys, [Ks[j]]; start_system=start_system)
        # solution = track_solution(hc_sys, out[j], [Ks[j]], [Ks[j + 1]])
        @debug "Done."
    end
end

function solve_all(hc_system, K_span; kwargs...)
    out = Vector{Any}(undef, length(K_span))
    solve_all!(out, hc_system, K_span; kwargs...)
    return out
end

# function solve_all(hc_system, Ks;
#                    variables = HC.variables(hc_system),
#                    start_system = :total_degree,
#                    show_progress = false,
#                    kwargs...)
#     return solve_all!(Vector{Any}(undef, length(Ks)), hc_system, Ks;
#                       variables = variables,
#                       start_system = start_system,
#                       show_progress = show_progress,
#                       kwargs...)
# end

## Defining hc_system and hc_homo
# hc_eqs_smooth = F_smooth(hc_vars[new_ind], hc_params)
# hc_system = HC.System(hc_eqs_smooth; variables=hc_vars[new_ind], parameters=[hc_param]);
# hc_homo = HC.Homotopy(hc_eqs_smooth, hc_vars[new_ind], hc_param)

function solve_branch!(out, hc_homo, Xs_init, K_init;
                       K_min = 0.0, K_max = 1.0,
                      dsmin = 1e-8, dsmax = 1e-2, ds = 1e-3,
                      variables = HC.variables(hc_system),
                      show_progress = false,
                      kwargs...)
    tracker = HC.Tracker(hc_homo)
    # track(tracker::Tracker, x::AbstractVector, t₁ = 1.0, t₀ = 0.0; debug::Bool = false)
    k = K_init
    x = X_init
    push!(out, (k, x))
    dds = ds
    k_new = k + dds
    while k <= K_max
        result = HC.track(tracker, x, t₁ = k, t₀ = k_new)
        if HC.is_success(result)
            k = k_new
            x = HC.solution(result)
            push!(out, (k, x))
            dds = min(dsmax, 1.1*dds)
            k_new += dds
        else
            dds = max(dsmin, dds/10)
            k_new = k + dds
        end
    end
    return out
end

# function generatedata_tracking(hc_sys, Ks, start_system=:total_degree)
#     @debug "Preallocate output... "
#     out = Vector{Any}(undef, length(Ks))
#     @debug "Done."
#     ## Compute initial solution to track along all the values of K
#     @debug "Compute initial solutions"
#     out[1] = starting_solution(hc_sys, [first(Ks)]; start_system=start_system)[1]
#     @debug "Done."
#     for j in 1:(length(Ks) - 1)
#         @debug "Iteration $j, point $(out[j]) and Ks $(Ks[j])"
#         solution = track_solution(hc_sys, out[j], [Ks[j]], [Ks[j + 1]])
#         @debug "Done."
#         @debug "Solution computed $(solution)"
#         if length(solution) < 1
#             @debug "Recomputing Initial point"
#             solution = starting_solution(hc_sys, [Ks[j + 1]])[1]
#         end
#         out[j + 1] = solution[1]
#     end
#     return out
# end


text:
export DummyDAQ, DummyDAQParams

Base.@kwdef mutable struct DummyDAQParams <: DeviceParams
  # Here "only" for user/toml parameter
  samplesPerPeriod::Int
  amplitude::Float32 = 1.0
  frequency::Float32 = 1.0
end
DummyDAQParams(dict::Dict) = params_from_dict(DummyDAQParams, dict)

Base.@kwdef mutable struct DummyDAQ <: AbstractDAQ
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::DummyDAQParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  # Here for any other internal parameters
end

function init(daq::DummyDAQ)
  @debug "Initializing dummy DAQ with ID `$(daq.deviceID)`."

  daq.present = true
end

neededDependencies(::DummyDAQ) = []
optionalDependencies(::DummyDAQ) = [TxDAQController, SurveillanceUnit]

Base.close(daq::DummyDAQ) = nothing

function startTx(daq::DummyDAQ)
end

function stopTx(daq::DummyDAQ)
end

# Version 1
# Parameter list can grow very large, but all are "known" in function signature
#=function setTxParams(daq::DummyDAQ, amplitude, frequency)
  daq.params.amplitude = amplitude # Only set them if they need to be remembered/used obvs
  daq.params.frequency = frequency
  doSomething(daq.params.amplitude)
  doSomething(daq.params.amplitude)
end
function setTxParams(daq::DummyDAQ, sequence::Sequence)
  temp = electricalTxChannels(sequence)
  channels = [channel for channel in temp if channel isa periodicElectricalTxChannels]
  component = components(channels[1])
  setTxParams(daq, amplitude(component), divider(component))
end=#

# Version 2
# Paraemter list is hidden, but if someone wants to call function directly they need to know which parameters to set
function setTxParams(daq::DummyDAQ, sequence::Sequence)
  temp = electricalTxChannels(sequence)
  channels = [channel for channel in temp if channel isa periodicElectricalTxChannels]
  component = components(channels[1])
  daq.params.amplitude = amplitude(component)
  daq.params.frequency = divider(component)
  daq.state = 1
  setTxParams(daq)
end
function setTxParams(daq::DummyDAQ)
  doSomething(daq.params.amplitude)
  doSomething(daq.params.amplitude)
end
# Or maybe an additional helper function, could also be used in f(daq, seq) then
function setTxParams(daq::DummyDAQ, amplitude, frequency)
  daq.params.amplitude = amplitude
  daq.params.frequency = frequency
  setTxParams(daq)
end

function currentFrame(daq::DummyDAQ)
  return 1;
end

function currentPeriod(daq::DummyDAQ)
	return 1;
end

function disconnect(daq::DummyDAQ)
end

enableSequence(daq::DummyDAQ, enable::Bool, numFrames=0,
              ffRampUpTime=0.4, ffRampUpFraction=0.8) = 1

function readData(daq::DummyDAQ, startFrame, numFrames)
    uMeas = zeros(daq.params.samplesPerPeriod,1,1,1)
    uRef = zeros(daq.params.samplesPerPeriod,1,1,1)

    uMeas[:,1,1,1] = sin.(range(0,2*pi, length=daq.params.samplesPerPeriod))
    uRef[:,1,1,1] = sin.(range(0,2*pi, length=daq.params.samplesPerPeriod))

    return uMeas, uRef
end

function readDataPeriods(daq::DummyDAQ, startPeriod, numPeriods)
    uMeas = zeros(daq.params.samplesPerPeriod,1,1,1)
    uRef = zeros(daq.params.samplesPerPeriod,1,1,1)

    uMeas[:,1,1,1] = sin.(range(0,2*pi, length=daq.params.samplesPerPeriod))
    uRef[:,1,1,1] = sin.(range(0,2*pi, length=daq.params.samplesPerPeriod))

    return uMeas, uRef
end
refToField(daq::DummyDAQ, d::Int64) = 0.0

mutable struct DummyAsyncBuffer <: AsyncBuffer
    samples::Union{Matrix{Float32}, Nothing}
end

function channelType(daq::DummyDAQ)
    return Matrix{Float32}
end

function AsyncBuffer(daq::DummyDAQ)
    return DummyAsyncBuffer(nothing)
end

function frameAverageBufferSize(daq::DummyDAQ, frameAverages)
    # 1 Rx channel and 1 period per frame
    return daq.params.rxNumSamplingPoints, 1, 1, frameAverages
end

function updateAsyncBuffer!(buffer::DummyAsyncBuffer, chunk)
    samples = chunk
    if !isnothing(buffer.samples)
        buffer.samples = hcat(buffer.samples, samples)
    else
        buffer.samples = samples
    end
end

function retrieveMeasAndRef!(buffer::DummyAsyncBuffer, daq::DummyDAQ)
    unusedSamples = buffer.samples
    samples = unusedSamples
    frames = nothing
    samplesInBuffer = size(samples)[2]
    framesInBuffer = div(samplesInBuffer, daq.params.samplesPerPeriod)

    if framesInBuffer > 0
        samplesToConvert = samples[:, 1:(daq.params.samplesPerPeriod * framesInBuffer)]
        temp = reshape(samplesToConvert, 2, daq.params.samplesPerPeriod, 1, 1)
        frames = zeros(Float32, daq.params.samplesPerPeriod, 2, 1, 1)
        frames[:, 1, :, :] = temp[1, :, :, :]
        frames[:, 2, :, :] = temp[2, :, :, :]

        if (daq.params.samplesPerPeriod * framesInBuffer) + 1 <= samplesInBuffer
            unusedSamples = samples[:, (daq.params.samplesPerPeriod * framesInBuffer) + 1:samplesInBuffer]
        else
          unusedSamples = nothing
        end

    end

    buffer.samples = unusedSamples

    uMeas = nothing
    uRef = nothing
    if !isnothing(frames)
        uMeas = frames[:, [1], :, :]
        uRef = frames[:, [2], :, :]
    end

    return uMeas, uRef
end

function startProducer(channel::Channel, daq::DummyDAQ)
    startTx(daq)
    startFrame = 1
    endFrame = numFrames + 1
    currentFrame = startFrame
    try
        while currentFrame < endFrame
            samples = zeros(Float32, 2, daq.params.samplesPerPeriod)
            samples[1, :] = daq.params.amplitude * sin.(daq.params.frequency .* range(0,2*pi, length=daq.params.samplesPerPeriod))
            samples[2, :] = daq.params.amplitude * sin.(daq.params.frequency .* range(0,2*pi, length=daq.params.samplesPerPeriod))
            put!(channel, samples)
            currentFrame += 1
        end
    catch e
        @error e
    end
    return endFrame
end

function prepareTx(daq::DummyDAQ, sequence::Sequence; allowControlLoop = true)
    # NOP
end

function setSequenceParams(daq::DummyDAQ, sequence::Sequence)
    # NOP
end

function prepareSequence(daq::DummyDAQ, sequence::Sequence)
    # NOP
end

function endSequence(daq::DummyDAQ, endFrame)
    # NOP
end

function readData(daq::DummyDAQ, startFrame, numFrames, numBlockAverages)
  uMeas=zeros(2,2,2,2)
  uRef=zeros(2,2,2,2)
  return uMeas, uRef
end

function readDataPeriods(daq::DummyDAQ, startPeriod, numPeriods, numBlockAverages)
  uMeas=zeros(2,2,2,2)
  uRef=zeros(2,2,2,2)
  return uMeas, uRef
end

numTxChannelsTotal(daq::DummyDAQ) = 1
numRxChannelsTotal(daq::DummyDAQ) = 1
numTxChannelsActive(daq::DummyDAQ) = 1
numRxChannelsActive(daq::DummyDAQ) = 1
numRxChannelsReference(daq::DummyDAQ) = 0
numRxChannelsMeasurement(daq::DummyDAQ) = 1

canPostpone(daq::DummyDAQ) = false
canConvolute(daq::DummyDAQ) = false


text:
abstract type AbstractCostFunctional end

"""
Default weight: 1 for position, (1/π^2) for angular velocity (1 m ≈ 180 deg)
# Notes
- Q_p: weight matrix for running cost (position)
- Q_ω: weight matrix for running cost (angular velocity)
- F_p: weight matrix for terminal cost (position)
- F_ω: weight matrix for terminal cost (angular velocity)
"""
struct PositionAngularVelocityCostFunctional <: AbstractCostFunctional
    Q_p::AbstractMatrix
    Q_ω::AbstractMatrix
    F_p::AbstractMatrix
    F_ω::AbstractMatrix
    function PositionAngularVelocityCostFunctional(
            Q_p=Matrix(I, 3, 3), Q_ω=Matrix(I, 3, 3),
            F_p=(1/π^2)*Matrix(I, 3, 3), F_ω=(1/π^2)*Matrix(I, 3, 3),
        )
        new(Q_p, Q_ω, F_p, F_ω)
    end
end

"""
ts: an array of time
e_ps: an array of position error
e_ωs: an array of angular velocity error
"""
function cost(cf::PositionAngularVelocityCostFunctional,
        ts::AbstractVector, e_ps::AbstractVector, e_ωs::AbstractVector,
    )
    @unpack Q_p, Q_ω, F_p, F_ω = cf
    running_costs = zip(e_ps, e_ωs) |> MapSplat((e_p, e_ω) -> e_p'*Q_p*e_p + e_ω'*Q_ω*e_ω) |> collect
    e_p_f = e_ps[end]
    e_ω_f = e_ωs[end]
    terminal_cost = e_p_f'*F_p*e_p_f + e_ω_f'*F_ω*e_ω_f
    integrate(ts, running_costs) + terminal_cost
end


text:
"""
    cartesian_tensor_component(t::Tensor{1}, c)

This returns the Cartesian tensor component `c` (valid choices are
`:x`, `:y`, or `:z`) of the rank-1 tensor `t`, as a linear combination
of its "natural" [`TensorComponent`](@ref)s `1`, `0`, and `-1`. The
transform matrix is `M(+1, 0, -1 ← x, y, z)` given in Table 1.2 on
p. 14 of Varshalovich (1988).

# Examples

```jldoctest
julia> cartesian_tensor_component(Gradient(), :x)
- 0.707107 𝛁̂⁽¹⁾₁ + 0.707107 𝛁̂⁽¹⁾₋₁
```
"""
function cartesian_tensor_component(t::Tensor{1}, c::Symbol)
    if c == :x
        (-TensorComponent(t, 1) + TensorComponent(t, -1))/√2
    elseif c == :y
        (im*TensorComponent(t, 1) + im*TensorComponent(t, -1))/√2
    elseif c == :z
        TensorComponent(t, 0)
    else
        throw(ArgumentError("Unknown Cartesian tensor component $(c)"))
    end
end

export cartesian_tensor_component


text:
"""
GRADIENT_DESCENT_1 A simple implementation of pure gradient descent

x = gradient_descent_1(f,x0) computes x, the argmin of f using a pure
gradient descent algorithm.  This routine SHOULD NOT be used for anything
other than pedagogical purposes as it does not include any technique to
make gradient descent a practical algorithm.

The function f(x0) must return both the function value and the gradient
at x0, i.e. [fx,gx] = f(x0);  See FCOMBINE for a helper function or
ROSENBROCK for an example.

[x,f,g,hist,histx] = gradient_descent_1(...) returns the final function
value, final gradient, a history of iteration data, and a history of
all the iterates x.  Memory for histx is only allocated if this output is
requested as it could grow rather large.

The function call
[...] = gradient_descent_1(f,x0;key=value)
provide optional arguments to the function that all have reasonable
default values set.

The options are:

   maxiter (Int) : the maximum number of iterations/function evaluations
   tol (Float32) : the stopping tolerance in terms of the infinity norm of the
     gradient
   quiet (Bool): do not display iteration output
   histx (Array) : an array that will be updated with the history of iterates,
     each is pushed on the back

Example:
using Optim
   x = gradient_descent_1(@rosenbrock,[1+0.1;1+0.1]); % does not converge
"""
function gradient_descent_1(fg,x0;
    maxiter=10000,tol=1.0e-8,quiet=false,histx=[],gamma=0.0001)

    x = copy(x0)
    n = length(x)

    hist = zeros(2,maxiter)
    savehistx = eltype(histx) == Vector{Float64} ? true : false

    f = Inf
    normg = Inf
    lastiter = 0
    g = Vector{Float64}()

    if !quiet
        @printf("  %6s  %9s %9s %9s\n", "iter",
            "val", "normg", "fdiff");
    end

    for iter=1:maxiter
        if savehistx
            push!(histx, x)
        end

        if iter>1
            x = x - gamma*g
        end

        flast = f
        f,g = fg(x)
        normg = norm(g,Inf)

        fdiff = flast - f

        if !quiet
            @printf("  %6i  %9.2e %9.2e %9.2e\n",
                iter, f, normg, fdiff)
        end

        hist[:,iter] = [f; normg]
        lastiter = iter

        if normg <= tol
            break
        end
        if !isfinite(normg)
            break
        end
    end

    if lastiter < maxiter
        hist = hist[:,1:lastiter]
    end

    if normg > tol
        warn("Did not converge")
    end

    return x,f,g,hist
end

"""
GRADIENT_DESCENT_2 A simple implementation of pure gradient descent
with normalized gradients

x = gradient_descent_1(f,x0) computes x, the argmin of f using a pure
gradient descent algorithm.  This routine SHOULD NOT be used for anything
other than pedagogical purposes as it does not include any technique to
make gradient descent a practical algorithm.

The function f should return both the function and gradient at the
provided point.

x,f,g,hist = gradient_descent_1(...) returns the final function
value, final gradient, a history of iteration data, and a history of
all the iterates x.

The function call
[...] = gradient_descent_1(f,x0;key=value)
provide optional arguments to the function that all have reasonable
default values set.

The options are:

   maxiter (Int) : the maximum number of iterations/function evaluations
   tol (Float32) : the stopping tolerance in terms of the infinity norm of the
     gradient
   quiet (Bool): do not display iteration output
   histx (Array) : an array that will be updated with the history of iterates,
     each is pushed on the back
   gamma (Float64) : A scaling of the _normalized_ gradient

Example:
using Optim
   x = gradient_descent_1(@rosenbrock,[1+0.1;1+0.1]); % does not converge
"""
function gradient_descent_2(fg,x0;
    maxiter=10000,tol=1.0e-8,quiet=false,histx=[],gamma=1.0)

    x = copy(x0)
    n = length(x)

    hist = zeros(2,maxiter)
    savehistx = eltype(histx) == Vector{Float64} ? true : false

    f = Inf
    normg = Inf
    lastiter = 0
    g = Vector{Float64}()

    if !quiet
        @printf("  %6s  %9s %9s %9s\n", "iter",
            "val", "normg", "fdiff");
    end

    for iter=1:maxiter
        if savehistx
            push!(histx, x)
        end

        if iter>1
            x = x - gamma*g/norm(g)
        end

        flast = f
        f,g = fg(x)
        normg = norm(g,Inf)

        fdiff = flast - f

        if !quiet
            @printf("  %6i  %9.2e %9.2e %9.2e\n",
                iter, f, normg, fdiff)
        end

        hist[:,iter] = [f; normg]
        lastiter = iter

        if normg <= tol
            break
        end
        if !isfinite(normg)
            break
        end
    end

    if lastiter < maxiter
        hist = hist[:,1:lastiter]
    end

    if normg > tol
        warn("Did not converge")
    end

    return x,f,g,hist
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CglBuilder"
version = v"0.59.10"

# Collection of sources required to build CglBuilder
sources = [
   "https://github.com/coin-or/Cgl/archive/releases/0.59.10.tar.gz" =>
   "2a4038dfeb00b4f4084bbe00b144b6bae4ddfe87873ac5913603af248cb3f1d4",
   "https://github.com/coin-or/Clp/archive/releases/1.16.11.tar.gz" =>
    "ac42c00ba95e1e034ae75ba0e3a5ff03b452191e0c9b2f5e2d5e65bf652fb0a1",

]

# Bash recipe for building across all platforms
script = raw"""

## STATIC BUILD START

# For static build, build and embedd a static version of Clp

if [ $target = "x86_64-apple-darwin14" ]; then
    export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

cd $WORKSPACE/srcdir
cd Clp-releases-1.16.11/
update_configure_scripts
# temporary fix
for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done
mkdir build
cd build/
export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H"

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
  --with-blas="-L${prefix}/lib -lcoinblas" \
  --with-lapack="-L${prefix}/lib -lcoinlapack" \
  --with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
  --with-mumps-lib="-L${prefix}/lib -lcoinmumps" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin"
else
  ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
  --with-blas="-L${prefix}/lib -lcoinblas" \
  --with-lapack="-L${prefix}/lib -lcoinlapack" \
  --with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
  --with-mumps-lib="-L${prefix}/lib -lcoinmumps" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
  LDFLAGS=-ldl;
fi

make -j${nproc}
make install

cd $WORKSPACE/srcdir
cd Cgl-releases-0.59.10/
update_configure_scripts
# temporary fix
for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done

mkdir build
cd build/

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  ../configure --prefix=$prefix --disable-pkg-config --with-pic --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
  --with-osiclp-lib="-L${prefix}/lib -lOsiClp" --with-osiclp-incdir="$prefix/include/coin"
else
  ../configure --prefix=$prefix --disable-pkg-config --with-pic --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
  --with-osiclp-lib="-L${prefix}/lib -lOsiClp" --with-osiclp-incdir="$prefix/include/coin" \
  LDFLAGS=-ldl;
fi

## STATIC BUILD END

## DYNAMIC BUILD START
#cd $WORKSPACE/srcdir
#cd Cgl-releases-0.59.10/
#update_configure_scripts
## temporary fix
#for path in ${LD_LIBRARY_PATH//:/ }; do
#    for file in $(ls $path/*.la); do
#        echo "$file"
#        baddir=$(sed -n "s|libdir=||p" $file)
#        sed -i~ -e "s|$baddir|'$path'|g" $file
#    done
#done
#
#mkdir build
#cd build/
#
# To fix an OSX toolchain bug. See https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/388
#if [[ ${target} == *darwin* ]]; then
#    export AR=/opt/${target}/bin/${target}-ar
#fi
#if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
#   export LDFLAGS="-L${prefix}/lib ${prefix}/lib/libClp.a"
#fi
#../configure --prefix=$prefix --disable-pkg-config --with-pic --host=${target} --enable-shared --disable-static \
#--enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
#--with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
#--with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
#--with-osiclp-lib="-L${prefix}/lib -lOsiClp" --with-osiclp-incdir="$prefix/include/coin"
## DYNAMIC BUILD END

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)
# To fix gcc4 bug in Windows
platforms = setdiff(platforms, [Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)), Windows(:i686, compiler_abi=CompilerABI(:gcc4))])
push!(platforms, Windows(:i686,compiler_abi=CompilerABI(:gcc6)))
push!(platforms, Windows(:x86_64,compiler_abi=CompilerABI(:gcc6)))

# The products that we will ensure are always built
products(prefix) = [
   LibraryProduct(prefix, "libCgl", :libCgl),
   LibraryProduct(prefix, "libOsiClp", :libOsiClp),
   LibraryProduct(prefix, "libClp", :libClp),
   LibraryProduct(prefix, "libClpSolver", :libClpSolver)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaOpt/OsiBuilder/releases/download/v0.107.9-1-static/build_OsiBuilder.v0.107.9.jl",
    "https://github.com/JuliaOpt/CoinUtilsBuilder/releases/download/v2.10.14-1-static/build_CoinUtilsBuilder.v2.10.14.jl",
    "https://github.com/JuliaOpt/COINMumpsBuilder/releases/download/v1.6.0-1-static/build_COINMumpsBuilder.v1.6.0.jl",
    "https://github.com/JuliaOpt/COINMetisBuilder/releases/download/v1.3.5-1-static/build_COINMetisBuilder.v1.3.5.jl",
    "https://github.com/JuliaOpt/COINLapackBuilder/releases/download/v1.5.6-1-static/build_COINLapackBuilder.v1.5.6.jl",
    "https://github.com/JuliaOpt/COINBLASBuilder/releases/download/v1.4.6-1-static/build_COINBLASBuilder.v1.4.6.jl",
    "https://github.com/JuliaOpt/ASLBuilder/releases/download/v3.1.0-1-static/build_ASLBuilder.v3.1.0.jl"#,
 #  "https://github.com/juan-pablo-vielma/ClpBuilder/releases/download/v1.16.11-1-static/build_ClpBuilder.v1.16.11.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: juan-pablo-vielma/CglBuilder

max_issues_repo_name: juan-pablo-vielma/CglBuilder

max_forks_repo_name: juan-pablo-vielma/CglBuilder


text:
function emit(io::IO, inst::Instruction, id_bound = 999; pad_assignment = false)
    if !isnothing(inst.result_id)
        printstyled(io, '%', inst.result_id; color=:yellow)
        print(io, " = ")
    end

    printstyled(io, inst.opcode; color=:light_cyan)

    !isempty(inst.arguments) && print(io, '(')

    arginfo = collect(zip(inst.arguments, operand_kinds(inst)))
    !isempty(arginfo) && emit_argument(io, popfirst!(arginfo)...)
    while !isempty(arginfo)
        print(io, ", ")
        emit_argument(io, popfirst!(arginfo)...)
    end

    !isempty(inst.arguments) && print(io, ')')

    if !isnothing(inst.type_id)
        printstyled(io, "::"; color=:light_green, bold=true)
        printstyled(io, '%', inst.type_id; color=:light_green, bold=true)
    end
end

function emit_argument(io, arg::AbstractVector, kind)
    category = kind_to_category[kind]
    args = copy(arg)
    emit_argument(io, popfirst!(args), kind, category)
    while !isempty(args)
        print(io, ", ")
        emit_argument(io, popfirst!(args), kind, category)
    end
end

function emit_argument(io, arg, kind)
    category = kind_to_category[kind]
    emit_argument(io, arg, kind, category)
end

function emit_argument(io, arg, kind, category)
    @match category begin
        "ValueEnum" => printstyled(io, replace(string(arg), Regex("^$(nameof(kind))") => ""); color=:208)
        "BitEnum" => printstyled(io, arg; color=:light_magenta)
        "Literal" => @match kind begin
            &LiteralString => printstyled(io, '"', arg, '"'; color=150)
            _ => printstyled(io, arg; color=153)
        end
        "Id" => printstyled(io, '%', arg; color=:yellow)
    end
end

show(io::IO, ::MIME"text/plain", inst::Instruction) = emit(io, inst)

"""
    disassemble(io, spir_module)

Transform the content of `spir_module` into a human-readable format and prints it to `io`.
"""
function disassemble(io::IO, mod::Module)
    if mod.magic_number == magic_number
        println(io, "SPIR-V")
    else
        println(io, "Magic number: ", mod.magic_number)
    end
    println(io, "Version: ", join([mod.version.major, mod.version.minor], "."))
    println(io, "Generator: ", hex(mod.generator_magic_number))
    println(io, "Bound: ", mod.bound)
    println(io, "Schema: ", mod.schema)
    println(io)

    padding(id) = length(string(mod.bound)) - (isnothing(id) ? -4 : length(string(id)))
    for inst ∈ mod.instructions
        print(io, ' '^padding(inst.result_id))
        emit(io, inst, mod.bound)
        println(io)
    end
end

hex(x) = "0x" * lpad(string(x, base=16), sizeof(x) * 2, '0')

disassemble(obj) = disassemble(stdout, obj)
disassemble(io::IO, mod::PhysicalModule) = disassemble(io, Module(mod))

show(io::IO, ::MIME"text/plain", mod::Module) = disassemble(io, mod)


text:
# This script solves the NKCapital model with a risk-adjusted linearization.
using RiskAdjustedLinearizations, JLD2, LinearAlgebra, Test, SparseArrays

# Settings
define_functions      = true
testing               = true         # check model's solution under default parameters against saved output
nlsolve_sparsity      = true         # exploit sparsity in calls to nlsolve
autodiff              = false
algorithm             = :relaxation
euler_equation_errors = false
test_price_dispersion = false        # check if price dispersion in steady state is always bounded below by 1
plot_irfs             = false
horizon               = 40    # horizon for IRFs
N_approx              = 5     # Number of periods ahead used for forward-difference equations
n_GH                  = 5     # number of nodes for Gauss-Hermite quadrature

if define_functions
    include("nk_with_capital.jl")
end

# Set up
m_nk = NKCapital(; N_approx = N_approx) # create parameters
m = nk_capital(m_nk) # instantiate risk-adjusted linearization
autodiff_method = autodiff ? :forward : :central

if nlsolve_sparsity
    solve!(m; algorithm = :deterministic, autodiff = autodiff_method) # ensure Ψ is nonzero
    jac_cache = preallocate_jac_cache(m, algorithm; sparsity_detection = false)
else
    jac_cache = nothing
end

# Solve!
solve!(m; algorithm = algorithm, autodiff = autodiff_method,
       jac_cache = jac_cache, sparse_jacobian = nlsolve_sparsity)

if testing
    out = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "nk_with_capital_output.jld2"), "r")

    test_m_nk = NKCapital(; N_approx = 1) # create parameters
    test_m = nk_capital(test_m_nk) # instantiate risk-adjusted linearization

    solve!(test_m; algorithm = :deterministic, verbose = :none)
    zdet = copy(test_m.z)
    ydet = copy(test_m.y)
    Psidet = copy(test_m.Ψ)
    @test test_m.z ≈ out["z_det"]
    @test test_m.y ≈ out["y_det"]
    @test test_m.Ψ ≈ out["Psi_det"]

    solve!(test_m; algorithm = :relaxation, verbose = :none)
    @test test_m.z ≈ out["z"]
    @test test_m.y ≈ out["y"]
    @test test_m.Ψ ≈ out["Psi"]

    test_5_m_nk = NKCapital(; N_approx = 5) # create parameters
    test_5_m = nk_capital(test_5_m_nk) # instantiate risk-adjusted linearization
    solve!(test_5_m; algorithm = :relaxation, verbose = :none)
    @test test_5_m.z ≈ out["z_5"]
    @test test_5_m.y ≈ out["y_5"]
    @test test_5_m.Ψ ≈ out["Psi_5"]
end

if test_price_dispersion
    π_ss_vec = log.(range(1 - .005, stop = 1 + .005, length = 10)) # On annualized basis, range from -2% to 2% target inflation
    det_soln = Dict()
    sss_soln = Vector{RiskAdjustedLinearization}(undef, length(π̃_ss_vec))

    for (i, π_ss) in enumerate(π_ss_vec)
        local m_nk = NKCapital(; π_ss = π_ss)
        local m = nk_capital(m_nk)
        solve!(m; algorithm = :deterministic, verbose = :none)
        det_soln[i] = Dict()
        det_soln[i][:z] = copy(m.z)
        det_soln[i][:y] = copy(m.y)
        det_soln[i][:Ψ] = copy(m.Ψ)
        solve!(m; algorithm = algorithm, verbose = :none)
        sss_soln[i] = m
    end

    det_v = exp.([det_soln[i][:z][3] for i in 1:length(det_soln)])
    sss_v = exp.([sss_soln[i].z[3] for i in 1:length(sss_soln)])
    @test all(det_v .> 1.)
    @test all(sss_v .> 1.)
end

if euler_equation_errors

    # Load shocks. Using CRW ones b/c that model also has 2 standard normal random variables
    shocks = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "nk_with_capital_shocks.jld2"), "r")["shocks"][:, 1:50]

    # Calculate Euler equation for bonds
    @test abs(euler_equation_error(m, nk_cₜ, (a, b, c, d) -> nk_log_euler(a, b, c, d; β = m_nk.β, γ = m_nk.γ, J = m_nk.J),
                                   nk_𝔼_quadrature, shocks, summary_statistic = x -> norm(x, Inf))) ≈ 0.

    # Can calculate the Euler equation error for q, s₁, and s₂ as well by treating these variables as "consumption variables"
    # but need to do the Euler equation error calculation "semi-manually" b/c of the forward difference equations
    impl_output = Dict()
    for k in [:dq, :pq, :ds₁, :ps₁, :ds₂, :ps₂]
        impl_output[k] = Dict()
    end

    _states, _jumps = simulate(m, size(shocks, 2), shocks, m.z)
    q_ral  = _jumps[m_nk.J[:q], :]
    s₁_ral = _jumps[m_nk.J[:s₁], :]
    s₂_ral = _jumps[m_nk.J[:s₂], :]

    for i in 1:m_nk.N_approx
        impl_output[:dq][i] = log.(euler_equation_error(m, (m, zₜ) -> nk_dqₜ(m, zₜ, i, m_nk.J),
                                                        (a, b, c, d) -> nk_log_dq(a, b, c, d; β = m_nk.β,
                                                                                γ = m_nk.β, i = i, J = m_nk.J, S = m_nk.S),
                                                        nk_𝔼_quadrature, shocks, return_soln = true))
        impl_output[:pq][i] = log.(euler_equation_error(m, (m, zₜ) -> nk_pqₜ(m, zₜ, i, m_nk.J),
                                                        (a, b, c, d) -> nk_log_pq(a, b, c, d; β = m_nk.β,
                                                                                γ = m_nk.β, i = i, J = m_nk.J, S = m_nk.S),
                                                        nk_𝔼_quadrature, shocks, return_soln = true))
        impl_output[:ds₁][i - 1] = log.((i == 1) ? [nk_ds₁ₜ(m, _states[:, t], i - 1, m_nk.J) for t in 1:size(_states, 2)] :
                                        euler_equation_error(m, (m, zₜ) -> nk_ds₁ₜ(m, zₜ, i - 1, m_nk.J),
                                                             (a, b, c, d) -> nk_log_ds₁(a, b, c, d; β = m_nk.β,
                                                                                γ = m_nk.β, θ = m_nk.θ, ϵ = m_nk.ϵ,
                                                                                i = i - 1, J = m_nk.J, S = m_nk.S),
                                                             nk_𝔼_quadrature, shocks, return_soln = true))
        impl_output[:ps₁][i] = log.(euler_equation_error(m, (m, zₜ) -> nk_ps₁ₜ(m, zₜ, i, m_nk.J),
                                                         (a, b, c, d) -> nk_log_ps₁(a, b, c, d; β = m_nk.β,
                                                                                γ = m_nk.β, θ = m_nk.θ, ϵ = m_nk.ϵ,
                                                                                i = i, J = m_nk.J, S = m_nk.S),
                                                         nk_𝔼_quadrature, shocks, return_soln = true))
        impl_output[:ds₂][i - 1] = log.((i == 1) ? [nk_ds₂ₜ(m, _states[:, t], i - 1, m_nk.J) for t in 1:size(_states, 2)] :
                                        euler_equation_error(m, (m, zₜ) -> nk_ds₂ₜ(m, zₜ, i - 1, m_nk.J),
                                                             (a, b, c, d) -> nk_log_ds₂(a, b, c, d; β = m_nk.β,
                                                                                γ = m_nk.β, θ = m_nk.θ, ϵ = m_nk.ϵ,
                                                                                i = i - 1, J = m_nk.J, S = m_nk.S),
                                                             nk_𝔼_quadrature, shocks, return_soln = true))
        impl_output[:ps₂][i] = log.(euler_equation_error(m, (m, zₜ) -> nk_ps₂ₜ(m, zₜ, i, m_nk.J),
                                                         (a, b, c, d) -> nk_log_ps₂(a, b, c, d; β = m_nk.β,
                                                                                γ = m_nk.β, θ = m_nk.θ, ϵ = m_nk.ϵ,
                                                                                i = i, J = m_nk.J, S = m_nk.S),
                                                         nk_𝔼_quadrature, shocks, return_soln = true))
    end

    q_impl = log.(sum([exp.(x) for x in collect(values(impl_output[:dq]))]) + exp.(impl_output[:pq][m_nk.N_approx]))
    s₁_impl = log.(sum([exp.(x) for x in collect(values(impl_output[:ds₁]))]) + exp.(impl_output[:ps₁][m_nk.N_approx]))
    s₂_impl = log.(sum([exp.(x) for x in collect(values(impl_output[:ds₂]))]) + exp.(impl_output[:ps₂][m_nk.N_approx]))

    @test maximum(abs.((exp.(q_impl) - exp.(q_ral)) ./ exp.(q_ral))) < .1
    @test maximum(abs.((exp.(s₁_impl) - exp.(s₁_ral)) ./ exp.(s₁_ral))) < .1
    @test maximum(abs.((exp.(s₂_impl) - exp.(s₂_ral)) ./ exp.(s₂_ral))) < .1
end

if plot_irfs
    # Show IRFs of interesting variables (discount rate, labor supply, productivity, and MP shocks)
    m_nk = NKCapital()
    m = nk_capital(m_nk)

    solve!(m; algorithm = algorithm, autodiff = autodiff_method)

    z_irfs = Dict()
    y_irfs = Dict()

    for k in keys(m_nk.SH)
        z_irfs[k], y_irfs[k] = impulse_responses(m, horizon, m_nk.SH[k], 1.) # 1 positive standard deviation shock
    end

    using Plots
    plot_dicts = Dict()

    for k in keys(m_nk.SH)
        plot_dicts[k] = Dict()

        plot_dicts[k][:output] = plot(1:horizon, y_irfs[k][m_nk.J[:output], :], label = "Output",
                                      linewidth = 3, color = :black)
        plot_dicts[k][:l] = plot(1:horizon, y_irfs[k][m_nk.J[:l], :], label = "Hours",
                                 linewidth = 3, color = :black)
        plot_dicts[k][:w] = plot(1:horizon, y_irfs[k][m_nk.J[:w], :], label = "Real Wage",
                                 linewidth = 3, color = :black)
        plot_dicts[k][:rk] = plot(1:horizon, y_irfs[k][m_nk.J[:rk], :], label = "Rental Rate of Capital",
                                  linewidth = 3, color = :black)
        plot_dicts[k][:k] = plot(1:horizon, z_irfs[k][m_nk.S[:k₋₁], :], label = "Capital Stock",
                                 linewidth = 3, color = :black)
        plot_dicts[k][:π] = plot(1:horizon, y_irfs[k][m_nk.J[:π], :], label = "Inflation",
                                 linewidth = 3, color = :black)
        plot_dicts[k][:q] = plot(1:horizon, y_irfs[k][m_nk.J[:q], :], label = "Price of Capital",
                                 linewidth = 3, color = :black)
        plot_dicts[k][:x] = plot(1:horizon, y_irfs[k][m_nk.J[:x], :], label = "Investment",
                                 linewidth = 3, color = :black)
        plot_dicts[k][:r] = plot(1:horizon, y_irfs[k][m_nk.J[:r], :], label = "Nominal Interest Rate",
                                 linewidth = 3, color = :black)

        # excess returns on capital (exploits properties of IRFs)
        EₜRₖₜ₊₁ = exp.(y_irfs[k][m_nk.J[:rk], 2:end] .+ m.y[m_nk.J[:rk]])
        EₜQₜ₊₁ = exp.(y_irfs[k][m_nk.J[:q], 2:end] .+ m.y[m_nk.J[:q]])
        EₜΩₜ₊₁ = exp.(y_irfs[k][m_nk.J[:ω], 2:end] .+ m.y[m_nk.J[:ω]])
        exc_ret = ((EₜRₖₜ₊₁ + EₜQₜ₊₁ .* EₜΩₜ₊₁) ./ exp.(y_irfs[k][m_nk.J[:q], 1:end - 1] .+ m.y[m_nk.J[:q]])) -
            (exp.(y_irfs[k][m_nk.J[:r], 1:end - 1] - y_irfs[k][m_nk.J[:π], 1:end - 1] .+
                  (m.y[m_nk.J[:r]] - m.y[m_nk.J[:π]]))) .-
            (((exp(m.y[m_nk.J[:rk]]) + exp.(m.y[m_nk.J[:q]] + m.y[m_nk.J[:ω]])) / exp.(m.y[m_nk.J[:q]])) -
             exp(m.y[m_nk.J[:r]] - m.y[m_nk.J[:π]]))
        plot_dicts[k][:real_excess_ret] = plot(1:(horizon - 1), exc_ret, label = "Real Excess Returns",
                                               linewidth = 3, color = :black)
    end
end

nothing