text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures
using Printf

@testset verbose=true "test_bileveltso_market" begin

    function create_instance(ECH, S, TS,
                            fo::Minute,
                            load_1::Vector{Vector{Float64}}, load_2::Vector{Vector{Float64}}, wind_1::Vector{Vector{Float64}},
                            preceding_decision_prod_1::Vector{Vector{PSCOPF.GeneratorState}},
                            preceding_decision_prod_2::Vector{Vector{PSCOPF.GeneratorState}},
                            limit::Float64=35.,
                            logs=nothing)
        network = PSCOPFFixtures.network_2buses(limit=limit)
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 200.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 200.,
                                                10000., 10.,
                                                Dates.Second(4*60*60), Dates.Second(15*60))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "prod_2_1", PSCOPF.Networks.PILOTABLE,
                                                20., 200.,
                                                50000., 50.,
                                                Dates.Second(1*60*60), Dates.Second(15*60))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()

        for (s_index, s) in enumerate(S)
            for (ts_index, ts) in enumerate(TS)
                PSCOPF.add_uncertainty!(uncertainties, ECH[1], "bus_1",  ts, s, load_1[s_index][ts_index]  )
                PSCOPF.add_uncertainty!(uncertainties, ECH[1], "bus_2",  ts, s, load_2[s_index][ts_index])
                PSCOPF.add_uncertainty!(uncertainties, ECH[1], "wind_1_1", ts, s, wind_1[s_index][ts_index])
            end
        end

        gen_initial_state = SortedDict{String,PSCOPF.GeneratorState}(
                                "prod_1_1" => PSCOPF.OFF,
                                "prod_2_1" => PSCOPF.OFF,
                            )

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.ManagementMode("test_mode", fo),
                                    gen_initial_state,
                                    uncertainties, nothing, logs)

        # c.f. tso_bilevel_does_not_start_units_for_eod_reasons_from_missing
        # important in context.market_schedule for TSOBilevel, because tso bilevel won't start units for EOD reasons
        for (s_index, s) in enumerate(S)
            for (ts_index, ts) in enumerate(TS)
                PSCOPF.set_commitment_value!(context.market_schedule, "prod_1_1", ts, s, preceding_decision_prod_1[s_index][ts_index])
                PSCOPF.set_commitment_value!(context.tso_schedule, "prod_1_1", ts, s, preceding_decision_prod_1[s_index][ts_index])
                PSCOPF.set_commitment_value!(context.market_schedule, "prod_2_1", ts, s, preceding_decision_prod_2[s_index][ts_index])
                PSCOPF.set_commitment_value!(context.tso_schedule, "prod_2_1", ts, s, preceding_decision_prod_2[s_index][ts_index])
            end
        end

        PSCOPF.set_horizon_timepoints(context, ECH)
        PSCOPF.check(context)

        return context
    end

################################################################################################
#           TEST CASE
################################################################################################

    #=
        TS: [11h]
        S: [S1]
                        bus 1                   bus 2
                            |                      |
        (limitable) wind_1_1|       "1_2"          |
        Pmin=0, Pmax=200    |                      |
        Csta=0, Cprop=1     |                      |
        S1: 85              |----------------------|
                            |         35           |
                            |                      |
        (pilotable) prod_1_1|                      |(pilotable) prod_2_1
        Pmin=20, Pmax=200   |                      | Pmin=20, Pmax=200
        Csta=10k, Cprop=10  |                      | Csta=50k, Cprop=50
        OFF->ON             |                      | OFF->OFF
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        S1: 50              |                      | S1: 70
    =#
    @testset "test_two_modes_on_a_single_ech_a_single_ts" begin
        ECH = [DateTime("2015-01-01T07:00:00"), DateTime("2015-01-01T10:00:00")]
        S = ["S1"]
        TS = [DateTime("2015-01-01T11:00:00")]
        future_ech = ECH[2]

        #[[S1_TS1, S1_TS_2, ...], [S2_TS1, S2_TS_2, ...]]
        load_1 = [[50.]]
        load_2 = [[70.]]
        wind_1 = [[85.]]
        test2_preceding_decision_prod_1 = [[PSCOPF.ON]]
        test2_preceding_decision_prod_2 = [[PSCOPF.OFF]]

        #=
        TS: [11h]
        S: [S1]
                        bus 1                   bus 2
                            |                      |
        (limitable) wind_1_1|       "1_2"          |
        Pmin=0, Pmax=200    |                      |
        Csta=0, Cprop=1     |                      |
        S1: 85              |----------------------|
                            |         35           |
                            |                      |
        (pilotable) prod_1_1|                      |(pilotable) prod_2_1
        Pmin=20, Pmax=200   |                      | Pmin=20, Pmax=200
        Csta=10k, Cprop=10  |                      | Csta=50k, Cprop=50
        OFF->ON             |                      | OFF->OFF
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        S1: 50              |                      | S1: 70

        MARKET:
        demand = 120
        market chooses economically,
        It uses the limitable prod => P_wind_1_1 = 85.
        remaining demand : 120-85 = 35MW
        Then the cheaper unit (start and prop cost) => prod_1_1 => P_prod_1_1 = 35.

        This causes a flow of 70 (>35.) on branch_1_2, the TSO needs to react

        TSO:
        The tso solves the RSO constraint while trying to stick to the market levels

        option 1 :
            -35 on prod_1_1, +35 on prod_2_1
            => Delta = 70 => cost_1 = 70.
            Pwind_1_1 = 85, P_prod_2_1 = 35
            cost = 85 + 35. * 50 + 50000 = 51835
        option 2 :
            -20 on wind_1_1, -15 on prod_1_1, +35 on prod_2_1
            => Delta = 70 => cost_1 = 70. + 1e-3 (cause one limitation)
            cost_step1 = 70.001 => option 1 is better
            Pwind_1_1 = 65, P_prod_2_1 = 20, P_prod_2_1 = 35
            cost = 65 + 20. * 10 + 35. * 50 + 50000 = 52015

        Market Decision :
            P_wind_1_1 = 85
            P_prod_1_1 = 35
            P_prod_2_1 = 0
        TSO reaction :
            P_wind_1_1 = 85 (unchanged)
            P_prod_1_1 = 0 (-35)
            P_prod_2_1 = 35 (+35)

        =#
        @testset "energymarket_tso" begin

            context = create_instance(ECH, S, TS, Dates.Minute(1),
                                    load_1,
                                    load_2,
                                    wind_1,
                                    test2_preceding_decision_prod_1,
                                    test2_preceding_decision_prod_2,
                                    35.)

            @testset "energymarket_before_tso" begin
                PSCOPF.run_step!(context, PSCOPF.EnergyMarket(), ECH[1], ECH[2])

                @test 85. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
                @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1")
                @test PSCOPF.get_prod_value(context.market_schedule, "prod_2_1", TS[1], "S1") < 1e-09
            end

            @testset "tso" begin
                PSCOPF.run_step!(context, PSCOPF.TSOOutFO(), ECH[1], ECH[2])

                # market schedule is not changed
                @test 85. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
                @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1")
                @test PSCOPF.get_prod_value(context.market_schedule, "prod_2_1", TS[1], "S1") < 1e-09

                # TSO reaction
                @test 85. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1")
                @test PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1") < 1e-09
                @test 35. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_2_1", TS[1], "S1")

                #no limitation #Note : even if there was, they would only appear for units that reached the DP
                @test isempty(PSCOPF.get_limitations(context.tso_actions))
                #impositions only appear at DP
                @test isempty(PSCOPF.get_impositions(context.tso_actions))
                @test ismissing(PSCOPF.get_imposition(context.tso_actions, "prod_1_1", TS[1], "S1"))
                @test ismissing(PSCOPF.get_imposition(context.tso_actions, "prod_2_1", TS[1], "S1"))

            end

        end

        #=
        TS: [11h]
        S: [S1]
                        bus 1                   bus 2
                            |                      |
        (limitable) wind_1_1|       "1_2"          |
        Pmin=0, Pmax=200    |                      |
        Csta=0, Cprop=1     |                      |
        S1: 85              |----------------------|
                            |         35           |
                            |                      |
        (pilotable) prod_1_1|                      |(pilotable) prod_2_1
        Pmin=20, Pmax=200   |                      | Pmin=20, Pmax=200
        Csta=10k, Cprop=10  |                      | Csta=50k, Cprop=50
        OFF->ON             |                      | OFF->OFF
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        S1: 50              |                      | S1: 70

        MARKET:
        demand = 120
        market chooses economically,
        It uses the limitable prod => P_wind_1_1 = 85.
        remaining demand : 120-85 = 35MW
        Then the cheaper unit => prod_1_1 => P_prod_1_1 = 35.
        This would cause a flow of 70 (>35.) on branch_1_2, the TSO needs to react

        TSO:
        The TSO needs to anticipate to prevent the market from violating the RSO constraint
        => needs to limit the flow on branch_1_2 to 35.

        option 1 : TSO cost=15.
            produce at least 35.MW on bus 2
            => P_prod_2_1 \in 35-200 => (unit was off) cost=35
            to avoid any impositions costs, P_prod_1_1 will stay at [20, 200] (imposition cost would be 180)
            Since P_prod_1_1 will produce 20MW but bus 1's production should not exceed 85MW
            TSO will also limit wind_1_1 to 65MW

            The market will need to cap 20MW on wind_1_1

        option 2 : cost=180.
            produce at most 85MW on bus 1
            => P_prod_1_1 is OFF => cost= (200-0) + (0-20) = 180
            to avoid any impositions costs, P_prod_1_1 will stay at [20, 200]

        TSO decisions : cost = 15
            P_wind_1_1 : <= 65 (due to EOD not to TSO restrictions)
            P_prod_1_1 : 20 - 200
            P_prod_2_1 : 35 - 200
        Market decisions : cost = 65 + 20*10 + 35*50 = 2015 (ignoring start cost)
            P_wind_1_1 = 65
            P_prod_1_1 = 20
            P_prod_2_1 = 35

        FIXME? we prefer ?
        TSO decisions : cost = currently 180 but should be 0 ?
            P_wind_1_1 : no limit
            P_prod_1_1 : 0 - 0 i.e. OFF
            P_prod_2_1 : 20 - 200
        Market decisions : cost = 85 + 0 + 35*50 = 1835 (ignoring start cost)
            P_wind_1_1 = 85
            P_prod_1_1 = 0.
            P_prod_2_1 = 35
        =#
        @testset "tsobilevel_balancemarket" begin

            context = create_instance(ECH, S, TS, Dates.Minute(1),
                                    load_1,
                                    load_2,
                                    wind_1,
                                    test2_preceding_decision_prod_1,
                                    test2_preceding_decision_prod_2,
                                    35.)

            @testset "tso_bilevel" begin
                PSCOPF.run_step!(context, PSCOPF.TSOBilevel(), ECH[1], ECH[2])

                @test 65. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1")
                @test 20. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
                @test 35. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_2_1", TS[1], "S1")

                #no limitation #Note : even if there was, they would only appear for units that reached the DP
                @test 65. ≈ PSCOPF.get_limitation(context.tso_actions, "wind_1_1", TS[1])
                #impositions only appear at DP
                @test !isempty(PSCOPF.get_impositions(context.tso_actions))
                @test 20. ≈ PSCOPF.get_imposition(context.tso_actions, "prod_1_1", TS[1], "S1")[1]
                @test 200. ≈ PSCOPF.get_imposition(context.tso_actions, "prod_1_1", TS[1], "S1")[2]
                @test 35. ≈ PSCOPF.get_imposition(context.tso_actions, "prod_2_1", TS[1], "S1")[1]
                @test 200. ≈ PSCOPF.get_imposition(context.tso_actions, "prod_2_1", TS[1], "S1")[2]

                @test 20. ≈ PSCOPF.get_capping(context.tso_schedule, "wind_1_1", TS[1], "S1")
                @test PSCOPF.get_loss_of_load(context.tso_schedule, "bus_1", TS[1], "S1") < 1e-09
                @test PSCOPF.get_loss_of_load(context.tso_schedule, "bus_2", TS[1], "S1") < 1e-09
            end

            @testset "balancemarket" begin
                PSCOPF.run_step!(context, PSCOPF.BalanceMarket(), ECH[1], ECH[2])

                @test 65. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
                @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1")
                @test 35. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_2_1", TS[1], "S1")

                @test 20. ≈ PSCOPF.get_capping(context.market_schedule, "wind_1_1", TS[1], "S1")
                @test PSCOPF.get_loss_of_load(context.market_schedule, "bus_1", TS[1], "S1") < 1e-09
                @test PSCOPF.get_loss_of_load(context.market_schedule, "bus_2", TS[1], "S1") < 1e-09
            end

        end

    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
DiffEqBase.@def iipnlcachefields begin
  nlcache = alg.nlsolve.cache
  @unpack κ,tol,max_iter,min_iter,new_W = nlcache
  z = similar(u)
  dz = similar(u); tmp = similar(u); b = similar(u)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  uToltype = real(uEltypeNoUnits)
  ηold = one(uToltype)

  if typeof(alg.nlsolve) <: NLNewton
    if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype !== nothing
      W = WOperator(f, zero(t))
      J = nothing # is J = W.J better?
    else
      J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
      W = similar(J)
    end
    du1 = zero(rate_prototype)
    uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
    linsolve = alg.linsolve(Val{:init},uf,u)
    z₊ = z
  elseif typeof(alg.nlsolve) <: NLFunctional
    J = nothing
    W = nothing
    du1 = rate_prototype
    uf = nothing
    jac_config = nothing
    linsolve = nothing
    z₊ = similar(z)
  end

  if κ !== nothing
    κ = uToltype(nlcache.κ)
  else
    κ = uToltype(1//100)
  end
  if tol !== nothing
    tol = uToltype(nlcache.tol)
  else
    reltol = 1e-1 # TODO: generalize
    tol = uToltype(min(0.03,first(reltol)^(0.5)))
  end
end

DiffEqBase.@def oopnlcachefields begin
  nlcache = alg.nlsolve.cache
  @unpack κ,tol,max_iter,min_iter,new_W = nlcache
  z = uprev
  uf = alg.nlsolve isa NLNewton ? DiffEqDiffTools.UDerivativeWrapper(f,t,p) : nothing
  uToltype = real(uEltypeNoUnits)
  ηold = one(uToltype)
  if DiffEqBase.has_jac(f) && alg.nlsolve isa NLNewton
    J = f.jac(uprev, p, t)
    if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
      J = DiffEqArrayOperator(J)
    end
    W = WOperator(f.mass_matrix, zero(t), J)
  else
    W = typeof(u) <: Number ? u : Matrix{uEltypeNoUnits}(undef, 0, 0) # uEltype?
  end

  if κ !== nothing
    κ = uToltype(κ)
  else
    κ = uToltype(1//100)
  end
  if tol === nothing
    reltol = 1e-1 # TODO: generalize
    tol = uToltype(min(0.03,first(reltol)^(0.5)))
  end
  z₊,dz,tmp,b,k = z,z,z,z,rate_prototype
end

DiffEqBase.@def iipnlsolve begin
  _nlsolve = typeof(alg.nlsolve).name.wrapper
  nlcache = NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,γ,c,ηold,z₊,dz,tmp,b,k)
  nlsolve = _nlsolve{true, typeof(nlcache)}(nlcache)
end
DiffEqBase.@def oopnlsolve begin
  _nlsolve = typeof(alg.nlsolve).name.wrapper
  nlcache = NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,γ,c,ηold,z₊,dz,tmp,b,k)
  nlsolve = _nlsolve{false, typeof(nlcache)}(nlcache)
end

@cache mutable struct ImplicitEMCache{uType,rateType,JType,WType,JC,UF,N,noiseRateType,F,dWType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  tmp::uType
  gtmp::noiseRateType
  gtmp2::rateType
  J::JType
  W::WType
  jac_config::JC
  linsolve::F
  uf::UF
  dW_cache::dWType
  nlsolve::N
end

function alg_cache(alg::ImplicitEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
  @iipnlcachefields

  gtmp = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    gtmp2 = gtmp
    dW_cache = nothing
  else
    gtmp2 = zero(rate_prototype)
    dW_cache = zero(ΔW)
  end

  γ, c = alg.theta,zero(t)
  @iipnlsolve
  ImplicitEMCache(u,uprev,du1,fsalfirst,k,z,dz,tmp,gtmp,gtmp2,J,W,jac_config,linsolve,uf,
                  dW_cache,nlsolve)
end

mutable struct ImplicitEMConstantCache{F,N} <: StochasticDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::ImplicitEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{false}})
  @oopnlcachefields
  γ, c = alg.theta,zero(t)
  @oopnlsolve
  ImplicitEMConstantCache(uf,nlsolve)
end

@cache mutable struct ImplicitEulerHeunCache{uType,rateType,JType,WType,JC,UF,N,noiseRateType,F,dWType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  tmp::uType
  gtmp::noiseRateType
  gtmp2::rateType
  gtmp3::noiseRateType
  J::JType
  W::WType
  jac_config::JC
  linsolve::F
  uf::UF
  nlsolve::N
  dW_cache::dWType
end

u_cache(c::ImplicitEulerHeunCache)    = (c.uprev2,c.z,c.dz)
du_cache(c::ImplicitEulerHeunCache)   = (c.k,c.fsalfirst)

function alg_cache(alg::ImplicitEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
  @iipnlcachefields
  gtmp = zero(noise_rate_prototype)
  gtmp2 = zero(rate_prototype)

  if is_diagonal_noise(prob)
      gtmp3 = gtmp2
      dW_cache = nothing
  else
      gtmp3 = zero(noise_rate_prototype)
      dW_cache = zero(ΔW)
  end

  γ, c = alg.theta,zero(t)
  @iipnlsolve
  ImplicitEulerHeunCache(u,uprev,du1,fsalfirst,k,z,dz,tmp,gtmp,gtmp2,gtmp3,
                         J,W,jac_config,linsolve,uf,nlsolve,dW_cache)
end

mutable struct ImplicitEulerHeunConstantCache{F,N} <: StochasticDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::ImplicitEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{false}})
  @oopnlcachefields
  γ, c = alg.theta,zero(t)
  @oopnlsolve
  ImplicitEulerHeunConstantCache(uf,nlsolve)
end

@cache mutable struct ImplicitRKMilCache{uType,rateType,JType,WType,JC,UF,N,noiseRateType,F} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  tmp::uType
  gtmp::noiseRateType
  gtmp2::noiseRateType
  gtmp3::noiseRateType
  J::JType
  W::WType
  jac_config::JC
  linsolve::F
  uf::UF
  nlsolve::N
end

function alg_cache(alg::ImplicitRKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
  @iipnlcachefields
  gtmp = zero(noise_rate_prototype)
  gtmp2 = zero(rate_prototype)
  gtmp3 = zero(rate_prototype)

  γ, c = alg.theta,zero(t)
  @iipnlsolve
  ImplicitRKMilCache(u,uprev,du1,fsalfirst,k,z,dz,tmp,gtmp,gtmp2,gtmp3,
                   J,W,jac_config,linsolve,uf,nlsolve)
end

mutable struct ImplicitRKMilConstantCache{F,N} <: StochasticDiffEqConstantCache
  uf::F
  nlsolve::N
end

function alg_cache(alg::ImplicitRKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{false}})
  @oopnlcachefields
  γ, c = alg.theta,zero(t)
  @oopnlsolve
  ImplicitRKMilConstantCache(uf,nlsolve)
end

ext: jl


text:
using lmp_str

file_out = "test.data"
file_pdb = "test.pdb"
file_psf = "test.psf"

data = converter_vmd(file_pdb, file_psf, potential_charmm36)
write_data(data, file_out)


text:
using ClimateMachine
using ClimateMachine.ConfigTypes
using ClimateMachine.Mesh.Topologies: BrickTopology
using ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid
using ClimateMachine.DGmethods: DGModel, init_ode_state, LocalGeometry
using ClimateMachine.DGmethods.NumericalFluxes:
    RusanovNumericalFlux,
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder
using ClimateMachine.ODESolvers
using ClimateMachine.GeneralizedMinimalResidualSolver:
    GeneralizedMinimalResidual
using ClimateMachine.VTK: writevtk, writepvtu
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.MPIStateArrays: euclidean_distance
using ClimateMachine.MoistThermodynamics:
    air_density, total_energy, internal_energy, soundspeed_air
using ClimateMachine.Atmos:
    AtmosModel,
    AtmosAcousticLinearModel,
    RemainderModel,
    NoOrientation,
    NoReferenceState,
    ReferenceState,
    DryModel,
    NoPrecipitation,
    NoRadiation,
    ConstantViscosityWithDivergence,
    vars_state_conservative
using ClimateMachine.VariableTemplates: @vars, Vars, flattenednames
import ClimateMachine.Atmos: atmos_init_aux!, vars_state_auxiliary

using CLIMAParameters
using CLIMAParameters.Planet: kappa_d
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates, Test

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const output_vtk = false

function main()
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    numlevels = integration_testing ? 4 : 1

    expected_error = Dict()
    expected_error[Float64, MRIGARKIRK21aSandu, 1] = 2.3236071337679274e+01
    expected_error[Float64, MRIGARKIRK21aSandu, 2] = 5.2652585224989430e+00
    expected_error[Float64, MRIGARKIRK21aSandu, 3] = 1.2100430848052603e-01
    expected_error[Float64, MRIGARKIRK21aSandu, 4] = 2.1974838909870273e-03

    expected_error[Float64, MRIGARKESDIRK34aSandu, 1] = 2.3235626679098608e+01
    expected_error[Float64, MRIGARKESDIRK34aSandu, 2] = 5.2672845223341218e+00
    expected_error[Float64, MRIGARKESDIRK34aSandu, 3] = 1.2097276468825705e-01
    expected_error[Float64, MRIGARKESDIRK34aSandu, 4] = 2.0920468129065205e-03

    @testset "$(@__FILE__)" begin
        for FT in (Float64,), dims in 2
            for mrigark_method in (MRIGARKIRK21aSandu, MRIGARKESDIRK34aSandu)
                @info @sprintf """Configuration
                                  ArrayType      = %s
                                  mrigark_method = %s
                                  FT             = %s
                                  dims           = %d
                                  """ ArrayType "$mrigark_method" "$FT" dims

                setup = IsentropicVortexSetup{FT}()
                errors = Vector{FT}(undef, numlevels)

                for level in 1:numlevels
                    numelems =
                        ntuple(dim -> dim == 3 ? 1 : 2^(level - 1) * 5, dims)
                    errors[level] = run(
                        mpicomm,
                        ArrayType,
                        polynomialorder,
                        numelems,
                        setup,
                        mrigark_method,
                        FT,
                        dims,
                        level,
                    )

                    @test errors[level] ≈
                          expected_error[FT, mrigark_method, level]
                end

                rates = @. log2(
                    first(errors[1:(numlevels - 1)]) /
                    first(errors[2:numlevels]),
                )
                numlevels > 1 && @info "Convergence rates\n" * join(
                    [
                        "rate for levels $l → $(l + 1) = $(rates[l])"
                        for l in 1:(numlevels - 1)
                    ],
                    "\n",
                )
            end
        end
    end
end

function run(
    mpicomm,
    ArrayType,
    polynomialorder,
    numelems,
    setup,
    mrigark_method,
    FT,
    dims,
    level,
)
    brickrange = ntuple(dims) do dim
        range(
            -setup.domain_halflength;
            length = numelems[dim] + 1,
            stop = setup.domain_halflength,
        )
    end

    topology = BrickTopology(
        mpicomm,
        brickrange;
        periodicity = ntuple(_ -> true, dims),
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorder,
    )

    model = AtmosModel{FT}(
        AtmosLESConfigType,
        param_set;
        orientation = NoOrientation(),
        ref_state = IsentropicVortexReferenceState{FT}(setup),
        turbulence = ConstantViscosityWithDivergence(FT(0)),
        moisture = DryModel(),
        source = nothing,
        boundarycondition = (),
        init_state_conservative = isentropicvortex_initialcondition!,
    )
    # This is a bad idea; this test is just testing how
    # implicit GARK composes with explicit methods
    # The linear model has the fast time scales but will be
    # treated implicitly (outer solver)
    slow_model = AtmosAcousticLinearModel(model)
    # The remainder will be treated explicitly in the inner loop
    fast_model = RemainderModel(model, (slow_model,))

    dg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )
    fast_dg = DGModel(
        fast_model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        state_auxiliary = dg.state_auxiliary,
    )
    slow_dg = DGModel(
        slow_model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        state_auxiliary = dg.state_auxiliary,
    )

    timeend = FT(2 * setup.domain_halflength / setup.translation_speed)

    # determine the time step
    elementsize = minimum(step.(brickrange))
    dt =
        elementsize / soundspeed_air(model.param_set, setup.T∞) /
        polynomialorder^2 / 5
    nsteps = ceil(Int, timeend / dt)
    dt = timeend / nsteps

    Q = init_ode_state(dg, FT(0))

    fastsolver = LSRK54CarpenterKennedy(fast_dg, Q; dt = dt)

    linearsolver = GeneralizedMinimalResidual(Q; M = 50, rtol = 1e-10)

    ode_solver = mrigark_method(
        slow_dg,
        LinearBackwardEulerSolver(linearsolver; isadjustable = true),
        fastsolver,
        Q;
        dt = dt,
        t0 = 0,
    )

    eng0 = norm(Q)
    dims == 2 && (numelems = (numelems..., 0))
    @info @sprintf """Starting refinement level %d
                      numelems  = (%d, %d, %d)
                      dt        = %.16e
                      norm(Q₀)  = %.16e
                      """ level numelems... dt eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            runtime = Dates.format(
                convert(DateTime, now() - starttime[]),
                dateformat"HH:MM:SS",
            )
            @info @sprintf """Update
                              simtime = %.16e
                              runtime = %s
                              norm(Q) = %.16e
                              """ gettime(ode_solver) runtime energy
        end
    end
    callbacks = (cbinfo,)

    if output_vtk
        # create vtk dir
        vtkdir =
            "vtk_isentropicvortex_mrigark" *
            "_poly$(polynomialorder)_dims$(dims)_$(ArrayType)_$(FT)" *
            "_$(FastMethod)_level$(level)"
        mkpath(vtkdir)

        vtkstep = 0
        # output initial step
        do_output(mpicomm, vtkdir, vtkstep, dg, Q, Q, model)

        # setup the output callback
        outputtime = timeend
        cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
            vtkstep += 1
            Qe = init_ode_state(dg, gettime(ode_solver))
            do_output(mpicomm, vtkdir, vtkstep, dg, Q, Qe, model)
        end
        callbacks = (callbacks..., cbvtk)
    end

    solve!(Q, ode_solver; timeend = timeend, callbacks = callbacks)

    # final statistics
    Qe = init_ode_state(dg, timeend)
    engf = norm(Q)
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished refinement level %d
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ level engf engf / eng0 engf - eng0 errf errf / engfe
    errf
end

Base.@kwdef struct IsentropicVortexSetup{FT}
    p∞::FT = 10^5
    T∞::FT = 300
    ρ∞::FT = air_density(param_set, FT(T∞), FT(p∞))
    translation_speed::FT = 150
    translation_angle::FT = pi / 4
    vortex_speed::FT = 50
    vortex_radius::FT = 1 // 200
    domain_halflength::FT = 1 // 20
end

struct IsentropicVortexReferenceState{FT} <: ReferenceState
    setup::IsentropicVortexSetup{FT}
end
vars_state_auxiliary(::IsentropicVortexReferenceState, FT) =
    @vars(ρ::FT, ρe::FT, p::FT, T::FT)
function atmos_init_aux!(
    m::IsentropicVortexReferenceState,
    atmos::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    setup = m.setup
    ρ∞ = setup.ρ∞
    p∞ = setup.p∞
    T∞ = setup.T∞

    aux.ref_state.ρ = ρ∞
    aux.ref_state.p = p∞
    aux.ref_state.T = T∞
    aux.ref_state.ρe = ρ∞ * internal_energy(atmos.param_set, T∞)
end

function isentropicvortex_initialcondition!(bl, state, aux, coords, t, args...)
    setup = bl.ref_state.setup
    FT = eltype(state)
    x = MVector(coords)

    ρ∞ = setup.ρ∞
    p∞ = setup.p∞
    T∞ = setup.T∞
    translation_speed = setup.translation_speed
    α = setup.translation_angle
    vortex_speed = setup.vortex_speed
    R = setup.vortex_radius
    L = setup.domain_halflength

    u∞ = SVector(translation_speed * cos(α), translation_speed * sin(α), 0)

    x .-= u∞ * t
    # make the function periodic
    x .-= floor.((x .+ L) / 2L) * 2L

    @inbounds begin
        r = sqrt(x[1]^2 + x[2]^2)
        δu_x = -vortex_speed * x[2] / R * exp(-(r / R)^2 / 2)
        δu_y = vortex_speed * x[1] / R * exp(-(r / R)^2 / 2)
    end
    u = u∞ .+ SVector(δu_x, δu_y, 0)

    _kappa_d::FT = kappa_d(param_set)
    T = T∞ * (1 - _kappa_d * vortex_speed^2 / 2 * ρ∞ / p∞ * exp(-(r / R)^2))
    # adiabatic/isentropic relation
    p = p∞ * (T / T∞)^(FT(1) / _kappa_d)
    ρ = air_density(bl.param_set, T, p)

    state.ρ = ρ
    state.ρu = ρ * u
    e_kin = u' * u / 2
    state.ρe = ρ * total_energy(bl.param_set, e_kin, FT(0), T)
end

function do_output(
    mpicomm,
    vtkdir,
    vtkstep,
    dg,
    Q,
    Qe,
    model,
    testname = "isentropicvortex_mrigark",
)
    ## name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state_conservative(model, eltype(Q)))
    exactnames = statenames .* "_exact"

    writevtk(filename, Q, dg, statenames, Qe, exactnames)

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(pvtuprefix, prefixes, (statenames..., exactnames...))

        @info "Done writing VTK: $pvtuprefix"
    end
end

main()


text:
using Test, StateSpaceModels, Statistics, CSV, LinearAlgebra, Random, Distributions

const SSM = StateSpaceModels

# Run tests
include("test_utils.jl")
include("test_userdefined.jl")
include("test_locallevel.jl")
include("test_lineartrend.jl")
include("test_structural.jl")
include("test_regression.jl")
include("test_simulate.jl")
include("test_opt_methods.jl")
include("test_results.jl")


text:
"""
    contact-implicit model-predictive control policy
"""

@with_kw mutable struct CIMPCOptions{T}
	altitude_update::Bool = false
	altitude_impact_threshold::T = 1.0
	altitude_verbose::Bool = false
    ip_max_time::T = 1e5     # maximum time allowed for an InteriorPoint solve
    live_plotting::Bool=false # Use the live plotting tool to debug
	gains::Bool=false
end

mutable struct CIMPC{T,NQ,NU,NW,NC,NB,NZ,Nθ,R,RZ,Rθ,Nν,W,FC,NQQ,NJ,NR,NI,OB,LS,NV} <: Policy{T}
	u::Vector{T}
	traj::ContactTraj{T,NQ,NU,NW,NC,NB,NZ,Nθ}
	traj_cache::ContactTraj{T,NQ,NU,NW,NC,NB,NZ,Nθ}
	ref_traj::ContactTraj{T,NQ,NU,NW,NC,NB,NZ,Nθ}
	K_traj::Vector{Matrix{T}}
	im_traj::ImplicitTrajectory{T,R,RZ,Rθ,NQQ}
	im_traj_cache::ImplicitTrajectory{T,R,RZ,Rθ,NQQ}
	H::Int
	stride::Vector{T}
	altitude::Vector{T}
	ϕ::Vector{T}
	κ::Vector{T}
	newton::Newton{T,NQ,NU,NW,NC,NB,NZ,Nθ,Nν,NJ,NR,NI,OB,LS,NV}
	newton_mode::Symbol
	s::Simulation{T,W,FC}
	q0::Vector{T}
	N_sample::Int
	cnt::Vector{Int}
	window::Vector{Int}
	opts::CIMPCOptions{T}
	buffer_time::T
	next_time_update::T
	times_reference::Vector{T}
	total_time_reference::T
end

function ci_mpc_policy(traj::ContactTraj, s::Simulation{T}, obj::Objective;
	H_mpc = traj.H,
	N_sample = 1,
	κ_mpc = traj.κ[1],
	mode = :configurationforce,
	newton_mode = :direct,
	n_opts = NewtonOptions(
		r_tol = 3e-4,
		max_iter = 5,
		verbose = false,
		live_plotting = false),
	mpc_opts = CIMPCOptions(),
	ip_opts = InteriorPointOptions(
				γ_reg = 0.1,
				undercut = 5.0,
				κ_tol = κ_mpc,
				r_tol = 1.0e-8,
				diff_sol = true,
				solver = :empty_solver,
				max_time = mpc_opts.ip_max_time,)
	) where T

	traj = deepcopy(traj)
	traj_cache = deepcopy(traj)
	ref_traj = deepcopy(traj)
	im_traj = ImplicitTrajectory(traj, s,
		κ = κ_mpc,
		max_time = mpc_opts.ip_max_time,
		opts=ip_opts,
		mode = mode)

	im_traj_cache = deepcopy(im_traj)

	stride = get_stride(s.model, traj)
	altitude = zeros(s.model.nc)
	ϕ = zeros(s.model.nc)
	if newton_mode == :direct
		newton = Newton(s, H_mpc, traj.h, traj, im_traj, obj = obj, opts = n_opts)
	elseif newton_mode == :structure
		newton = NewtonStructure(s, H_mpc, traj, obj, κ_mpc, opts = n_opts)
	else
		@error "invalid Newton solver specified"
	end

	window = zeros(Int, H_mpc + 2)

	K_traj = reference_gains(s, traj, obj)

	total_time_reference = traj.h * (traj.H - 1)
	times_reference = Array(range(0, stop=total_time_reference, length=traj.H))

	CIMPC(zeros(s.model.nu), traj, traj_cache, ref_traj, K_traj, im_traj, im_traj_cache,
		H_mpc, stride, altitude, ϕ, [κ_mpc], newton, newton_mode, s, copy(ref_traj.q[1]),
		N_sample, [N_sample], window, mpc_opts,
		0.0, 0.0,
		times_reference, total_time_reference)
end

function policy(p::CIMPC{T,NQ,NU,NW,NC}, traj::Trajectory{T}, t::Int) where {T,NQ,NU,NW,NC}
	# reset
	if t == 1
		p.cnt[1] = p.N_sample
		p.q0 .= p.ref_traj.q[1]
		p.altitude .= 0.0
		set_trajectory!(p.traj, p.ref_traj)
		set_implicit_trajectory!(p.im_traj, p.im_traj_cache)
		reset_window!(p.window)
	end

    if p.cnt[1] == p.N_sample
		# altitude_update
		(p.opts.altitude_update && t > 1) && (update_altitude!(p.altitude, p.ϕ, p.s,
									traj, t, NC, p.N_sample,
									threshold = p.opts.altitude_impact_threshold,
									verbose = p.opts.altitude_verbose))
		set_altitude!(p.im_traj, p.altitude)

		# # optimize
		q1 = traj.q[t+1]
		newton_solve!(p.newton, p.s, p.q0, q1,
			p.window, p.im_traj, p.traj, warm_start = t > 1)

		update!(p.im_traj, p.traj, p.s, p.altitude, p.κ[1], p.traj.H)

		# visualize
		p.opts.live_plotting && live_plotting(p.s.model, p.traj, traj, p.newton, p.q0, traj.q[t+1], t)

		# shift trajectory
		rot_n_stride!(p.traj, p.traj_cache, p.stride, p.window)

		# update
		update_window!(p.window, p.ref_traj.H)
		p.q0 .= q1

		# reset count
		p.cnt[1] = 0
    end

    p.cnt[1] += 1

	# scale control
	if p.newton_mode == :direct
		p.u .= p.newton.traj.u[1]
		p.u ./= p.N_sample
	elseif p.newton_mode == :structure
		p.u .= p.newton.u[1]
		p.u ./= p.N_sample
	else
		println("newton mode specified not available")
	end

	return p.u
end

function reset_window!(window)
	n = length(window)
	for i = 1:n
		window[i] = i
	end
	return
end

function update_window!(window, max_window)
	n = length(window)
	for i = 1:n
		window[i] += 1
		if window[i] > max_window
			window[i] = 1
		end
	end
	return window
end


text:
using Match
using DataStructures
include("CherokeeQuantumGates.jl")

# author: Younis Shah

stage_deque = StageOutputDeque(Deque{StageOutput}())

function map_inputs(input::String)
  mapreduce(state -> rstrip(lstrip(state)) == "zero" ? b₀ : b₁ , kron, split(input, " "))
end

function tensor(ψ::Array{Symbol})
  ψ = detect_gate(ψ)
  accumulator = ψ[1]
  for i = 2:length(ψ)
    accumulator = round_v1(kron(accumulator, ψ[i]))
  end
  accumulator
end

# Convert the string of quantum gates to a Symbol Table/Matrix
function create_symbolic_ckt_matrix(ckt_string::String, nqubits::Int64, stages::Int64)
  symbolic_ckt_mat = Array(Symbol, nqubits, stages)
  ckt = split(ckt_string,";")
  for column = 1:length(ckt)
    rows = split(ckt[column])
    for row = 1:length(rows)
      symbolic_ckt_mat[row, column] = symbol(rows[row])
    end
  end
  symbolic_ckt_mat
end

function detect_gate(gate_symbol::Symbol)
  @match gate_symbol begin
    :H => H
    :I => Ι
    :X => σₓ
    :Y => Y
    :Z => Z
    :S => S
     _ => "NA"
  end
end

@vectorize_1arg Symbol detect_gate

function generate_basis_vectors_v1(n::Int64)
  basis_vecs = [0 : 2^n-1]
  digits_ = Array(String, 2^n)
  for i in basis_vecs
    digits_[i+1] = bin(i,n)
  end
  digits_
end

function gate_CNOT_positions(column::Array{Symbol})
  control_list, target_pos = Int64[], 0
  for i = 1:length(column)
    if column[i] == :C
      push!(control_list, i)
    elseif column[i] == :N
      target_pos = i
    end
  end
  control_list,target_pos
end

function simulate_CNOT(basis_vectors::Array, control, target::Int64, nqubits::Int64)
  vectors = map(vec -> parseint(vec, 2), map(vec -> xor(vec, control, target), basis_vectors))
  CNOT = zeros(Int64, 2^nqubits, 2^nqubits) # as Sparse matrix
  j=1
  for i in vectors
    CNOT[i+1,j] = 1
    j +=1
  end
  CNOT
end

function xor(s::String, control_qubits, target_qubit::Int64)
  digits_ = split(s,"")
  and_of_c_qubits = reduce(&, create_control_qubit_array(digits_, control_qubits)) # "And" the qubits -> Ummr Parry
  digits_[target_qubit] = string(and_of_c_qubits $ int(digits_[target_qubit]))
  join(digits_)
end

function create_control_qubit_array(s::Array, control_qubits)
  c_qubits = Int64[]
  for c in control_qubits
    push!(c_qubits, int(s[c]))
  end
  c_qubits
end

function round_v1(accumulator::Array)
  if length(filter(x -> x != 0, imag(accumulator))) == 0 # means no imaginary part
      accumulator = round(accumulator, 3)
    else
      accumulator = complex(round(real(accumulator),3), round(imag(accumulator),3))
    end
  accumulator
end

# expects a column vector
function router(mat::Array{Symbol}, basis_vecs::Array, nqubits_stages::(Int64,Int64), stage_deque::StageOutputDeque)
  if :C in mat && :N in mat && !(:H in mat || :X in mat || :Y in mat || :Z in mat || :S in mat) # only CNOT
    control, target = gate_CNOT_positions(mat)
    stage = StageOutput(nqubits_stages[2], simulate_CNOT(basis_vecs, control, target, nqubits_stages[1]))
    push!(stage_deque.stages, stage)
  elseif !(:C in mat && :N in mat)  # No CNOT
    stage = StageOutput(nqubits_stages[2], tensor(mat))
    push!(stage_deque.stages, stage)
  else
    return "NA"
  end
end

function multiply_stages(stages::StageOutputDeque)
  accumulator = pop!(stages.stages).output;
  for i = 1:length(stages.stages)
    stage = pop!(stage_deque.stages).output
    accumulator = accumulator * stage
  end
 round_v1(accumulator)
end

function simulate(ckt::String, nqubits::Int64, stages::Int64)
  #print_with_color(:blue, "SIMULATING...")
  mat = create_symbolic_ckt_matrix(ckt, nqubits, stages)
  basis_vecs = generate_basis_vectors_v1(nqubits)
  for i = 1:stages
    if router(mat[:,i], basis_vecs, (nqubits,i), stage_deque) == "NA"
      return "NA"
    end
  end
 # print_with_color(:blue, "DONE\n")
  #saved_stages = deepcopy(stage_deque) ############################### for WEB
  multiply_stages(stage_deque) # OUTPUT! Yaaaaayyy!!!!
end

###############################################

#ckt = "H H H Y;Y X I I" # from web
##println(simulate(ckt, nqubits, stages))

################################################

max_stars_repo_name: younisshah/cherokee

max_issues_repo_name: younisshah/cherokee

max_forks_repo_name: younisshah/cherokee


text:
using ITensors

# TensorNetworkCircuit struct and public functions
export TensorNetworkCircuit, add_input!, add_output!, qubits
export contract_tn!, decompose_tensor!, create_test_tnc
export input_tensors, output_tensors, input_indices, output_indices, single_amplitude

const input_output_tensors = Dict{Char,Array{Float64, 1}}('0' => [1., 0.],
                                                          '1' => [0., 1.],
                                                          '+' => [1., 1.]./sqrt(2),
                                                          '-' => [1., -1.]./sqrt(2))

"""Tensor network circuit data-structure"""
mutable struct TensorNetworkCircuit
    qubits::Int64
    input_indices::Array{Index, 1}
    output_indices::Array{Index, 1}
    input_tensors::Array{Union{Symbol, Nothing}, 1}
    output_tensors::Array{Union{Symbol, Nothing}, 1}
    tn::TensorNetwork
end

function TensorNetworkCircuit(qubits::Int64)
    input_indices = [Index(2, tags="Qubit $i") for i in 1:qubits]
    output_indices = copy(input_indices)
    input_tensors = Array{Union{Symbol, Nothing}, 1}(nothing, qubits)
    output_tensors = Array{Union{Symbol, Nothing}, 1}(nothing, qubits)
    TensorNetworkCircuit(qubits, input_indices, output_indices, input_tensors, output_tensors, TensorNetwork())
end

function Base.copy(tnc::TensorNetworkCircuit)
    tn = TensorNetwork()
    f = x -> x === nothing ? nothing : push!(tn, tnc[x]; tid=x)
    my_input_tensors = convert(Array{Union{Symbol, Nothing}}, f.(input_tensors(tnc)))
    my_output_tensors = convert(Array{Union{Symbol, Nothing}}, f.(output_tensors(tnc)))
    for x in keys(tnc)
        if !(x in input_tensors(tnc)) && !(x in output_tensors(tnc))
            push!(tn, tnc[x]; tid=x)
        end
    end
    tn.next_id = tnc.tn.next_id

    TensorNetworkCircuit(qubits(tnc),
                         copy(input_indices(tnc)),
                         copy(output_indices(tnc)),
                         my_input_tensors,
                         my_output_tensors,
                         tn)
end

Base.eltype(::TensorNetworkCircuit) = ITensor
Base.getindex(tnc::TensorNetworkCircuit, i::Symbol) = tnc.tn[i]
Base.getindex(tnc::TensorNetworkCircuit, i::T) where T <: Index = tnc.tn[i]
Base.iterate(tnc::TensorNetworkCircuit) = iterate(values(tnc))
Base.iterate(tnc::TensorNetworkCircuit, i) = iterate(values(tnc), i)
Base.keys(tnc::TensorNetworkCircuit) = keys(tnc.tn)
Base.length(tnc::TensorNetworkCircuit) = length(tnc.tn)
Base.show(io::IO, ::MIME"text/plain", tnc::TensorNetworkCircuit) = print(io, "TensorNetworkCircuit(qubits => $(qubits(tnc)), gates => $(length(tnc)))")
Base.values(tnc::TensorNetworkCircuit) = values(tnc.tn)

input_tensors(tnc::TensorNetworkCircuit) = tnc.input_tensors
output_tensors(tnc::TensorNetworkCircuit) = tnc.output_tensors
input_indices(tnc::TensorNetworkCircuit) = tnc.input_indices
output_indices(tnc::TensorNetworkCircuit) = tnc.output_indices
bonds(tnc::TensorNetworkCircuit) = bonds(tnc.tn)
qubits(tnc::TensorNetworkCircuit) = tnc.qubits
tensor_data(tnc::TensorNetworkCircuit, i; kwargs...) = tensor_data(tnc.tn, i; kwargs...)
hyperindices(tnc::TensorNetworkCircuit, i; kwargs...) = hyperindices(tnc.tn, i; kwargs...)
get_hyperedges(tnc::TensorNetworkCircuit) = get_hyperedges(tnc.tn)
disable_hyperindices!(tnc::TensorNetworkCircuit) = disable_hyperindices!(tnc.tn)

"""
    push!(tnc::TensorNetworkCircuit,
          qubits::Vector{Int64},
          data::Array{T, 2}) where T

Function to add a gate to the tensor network circuit given
the qubits it acts on and an array of the matrix elements
"""
function Base.push!(tnc::TensorNetworkCircuit,
                    qubits::Vector{Int64},
                    data::Array{T, 2};
                    decompose::Bool=true) where T
    gate_input_indices = tnc.output_indices[qubits]
    tnc.output_indices[qubits] = gate_output_indices = [prime(x) for x in gate_input_indices]
    if length(qubits) == 1 || !decompose
        indices = [gate_output_indices..., gate_input_indices...]
        @assert prod(size(data)) == prod(dim.(indices)) "Data matrix dimension does not match indices"
        data = reshape(data, Tuple(dim.(indices)))
        push!(tnc, indices, data)
    elseif length(qubits) == 2 && decompose
        data = reshape(data, Tuple([dim.(gate_output_indices)..., dim.(gate_input_indices)...]))
        A, B = decompose_gate(data)
        if size(A)[3] == size(B)[1] == 1
            A = reshape(A, Tuple(size(A)[1:2]))
            B = reshape(A, Tuple(size(B)[2:3]))
            a_indices = [gate_output_indices[1], gate_input_indices[1]]
            b_indices = [gate_output_indices[2], gate_input_indices[2]]
        else
            virtual_index = Index(size(A)[3])
            a_indices = [gate_output_indices[1], gate_input_indices[1], virtual_index]
            b_indices = [virtual_index, gate_output_indices[2], gate_input_indices[2]]
        end
        push!(tnc, a_indices, A)
        push!(tnc, b_indices, B)
    else

        @error("Gates spanning more than two qubits not supported with decomposition yet.")
    end
end

Base.push!(tnc::TensorNetworkCircuit, indices::Vector{<:Index}, data::AbstractArray; kwargs...) = push!(tnc.tn, indices, data; kwargs...)

"""
    delete!(tnc::TensorNetworkCircuit, tensor_id::Symbol)

Function to remove a tensor from a tensor network circuit.
"""
function Base.delete!(tnc::TensorNetworkCircuit, tensor_id::Symbol)
    replace!(tnc.input_tensors, tensor_id => nothing)
    replace!(tnc.output_tensors, tensor_id => nothing)
    delete!(tnc.tn, tensor_id)
end

"""
    push_input!(tnc::TensorNetworkCircuit, tensor::Array{Elt, 1}, pos::Int64) where Elt

Function to add a single input tensor to the tensor network circuit at the given position
"""
function push_input!(tnc::TensorNetworkCircuit, tensor::Array{Elt, 1}, pos::Int64) where Elt
    index = tnc.input_indices[pos]
    tensor_sym = tnc.input_tensors[pos]
    if tensor_sym !== nothing
        delete!(tnc, tensor_sym)
    end
    tnc.input_tensors[pos] = push!(tnc, [index], tensor; tid=tensor_sym)
end

"""
    push_output!(tnc::TensorNetworkCircuit, tensor::Array{Elt, 1}, pos::Int64) where Elt

Function to add a single output tensor to the tensor network circuit at the given position
"""
function push_output!(tnc::TensorNetworkCircuit, tensor::Array{Elt, 1}, pos::Int64) where Elt
    index = tnc.output_indices[pos]
    tensor_sym = tnc.output_tensors[pos]
    if tensor_sym !== nothing
        delete!(tnc, tensor_sym)
    end
    tnc.output_tensors[pos] = push!(tnc, [index], tensor; tid=tensor_sym)
end

"""
    add_input!(tnc::TensorNetworkCircuit; input::Union{String, Nothing}=nothing)

Function to add input tensors to the circuit
"""
function add_input!(tnc::TensorNetworkCircuit, input::Union{String, Nothing}=nothing)
    if input === nothing input = "0"^qubits(tnc) end
    [push_input!(tnc, input_output_tensors[input[pos]], pos) for pos in 1:qubits(tnc)]
end

"""
    add_output!(tnc::TensorNetworkCircuit; output::Union{String, Nothing}=nothing)

Function to add output tensors to the circuit
"""
function add_output!(tnc::TensorNetworkCircuit, output::Union{String, Nothing}=nothing)
    if output === nothing output = "0"^qubits(tnc) end
    [push_output!(tnc, input_output_tensors[output[pos]], pos) for pos in 1:qubits(tnc)]
end

contract_tn!(tnc::TensorNetworkCircuit, plan) = contract_tn!(tnc.tn, plan)

contract_tn(tnc::TensorNetworkCircuit, plan) = contract_tn!(copy(tnc), plan)

function decompose_tensor!(tnc::TensorNetworkCircuit, args...; kwargs...)
    decompose_tensor!(tnc.tn, args...; kwargs...)
end

simple_contraction(tnc::TensorNetworkCircuit) = simple_contraction(tnc.tn)
simple_contraction!(tnc::TensorNetworkCircuit) = simple_contraction!(tnc.tn)

"""
    create_test_tnc(;input::Union{String, Nothing}=nothing,
                    output::Union{String, Nothing}=nothing,
                    no_input::Bool=false,
                    no_output::Bool=false,
                    kwargs...)

Create a tensor network circuit for a small example circuit, 3 qubit ghz preparation circuit in this case
"""
function create_test_tnc(;input::Union{String, Nothing}=nothing,
                         output::Union{String, Nothing}=nothing,
                         no_input::Bool=false,
                         no_output::Bool=false,
                         kwargs...)
    tnc = TensorNetworkCircuit(3)
    if !no_input add_input!(tnc, input) end
    push!(tnc, [1], Gates.h(); kwargs...)
    push!(tnc, [2, 1], Gates.cx(); kwargs...)
    push!(tnc, [3, 2], Gates.cx(); kwargs...)
    if !no_output add_output!(tnc, output) end
    tnc
end


text:
module TestPDA
    using Base.Test
    using DataArrays
    using Compat

    p = @pdata [9, 9, 8, NA, 1, 1]
    pcopy = copy(p)
    @assert levels(p) == [1, 8, 9]
    @assert levels(setlevels(p, ["a", "b", "c"])) == ["a", "b", "c"]
    @assert dropna(setlevels(p, (@data ["a", "b", NA]))) == ["b", "a", "a"]
    @assert dropna(setlevels(p, (@data ["a", "b", "a"]))) == ["a", "a", "b", "a", "a"]
    @assert levels(setlevels(p, (@data ["a", "b", "a"]))) == ["a", "b"]
    @assert levels(setlevels(p, Dict([(1, 111)]))) == [111, 8, 9]
    @assert levels(setlevels(p, Dict([(1, 111), (8, NA)]))) == [111, 9]
    @assert levels(PooledDataArray(p, [9, 8, 1])) == [9, 8, 1]
    @assert levels(PooledDataArray(p, [9, 8])) == [9, 8]
    @assert dropna(PooledDataArray(p, [9, 8])) == [9, 9, 8]
    @assert levels(PooledDataArray(p, levels(p)[[3,2,1]])) == [9,8,1]
    v = collect(1:6)
    @assert isequal(p, reorder(p))
    # @assert levels(reorder(p, v)) == [9,8,1]
    @assert isequal(p, pcopy)

    @assert levels(setlevels!(copy(p), [10,80,90])) == [10, 80, 90]
    @assert levels(setlevels!(copy(p), [1,8,1])) == [1, 8]
    @assert levels(setlevels!(copy(p), (@data [1, 8, NA]))) == [1, 8]
    @assert levels(setlevels!(copy(p), [1,8,9, 10])) == [1, 8, 9, 10]
    @assert levels(setlevels!(copy(p), Dict([(1, 111)]))) == [111, 8, 9]
    @assert levels(setlevels!(copy(p), Dict([(1, 111), (8, NA)]))) == [111, 9]
    # issue #201
    @assert levels(setlevels!(@pdata([1.0, 2.0]), [3,4])) == [3.0, 4.0]

    y = @pdata [1, NA, -2, 1, NA, 4, NA]
    @assert isequal(unique(y), @pdata [1, NA, -2, 4])
    @assert isequal(unique(reverse(y)), @data [NA, 4, 1, -2])
    @assert isequal(unique(dropna(y)), @data [1, -2, 4])
    @assert isequal(unique(reverse(dropna(y))), @data [4, 1, -2])

    z = @pdata ["frank", NA, "gertrude", "frank", NA, "herbert", NA]
    @assert isequal(unique(z), @pdata ["frank", NA, "gertrude", "herbert"])
    @assert isequal(unique(reverse(z)), @pdata [NA, "herbert", "frank", "gertrude"])
    @assert isequal(unique(dropna(z)), @pdata ["frank", "gertrude", "herbert"])
    @assert isequal(unique(reverse(dropna(z))), @pdata ["herbert", "frank", "gertrude"])

    # check case where only NA occurs in final position
    @assert isequal(unique(@pdata [1, 2, 1, NA]), @pdata [1, 2, NA])

    pp = PooledDataArray(Any[])
    @assert length(pp) == 0
    @assert length(levels(pp)) == 0

    # test construction with unordered types
    pim = @pdata [1 + im, 2 + im, 3 + im, 2 + im, 1 + im]
    @assert levels(pim) == [1 + im, 2 + im, 3 + im]

    # Test explicitly setting refs type
    testarray = [1, 1, 2, 2, 0, 0, 3, 3]
    testdata = @data [1, 1, 2, 2, 0, 0, 3, 3]
    for t in Any[testarray, testdata]
        for R in [UInt8, UInt16, UInt32, UInt64]
            @assert eltype(PooledDataArray(t, R).refs) == R
            @assert eltype(PooledDataArray(t, [1,2,3], R).refs) == R
            @assert eltype(PooledDataArray(t, [1,2,3], t .== 0, R).refs) == R
        end
    end

    pcopy = copy(p)
    @assert levels(append!(pcopy, @pdata [4, NA, 6, 5])) == [1, 8, 9, 4, 5, 6]

    x = PooledDataArray([9, 9, 8])
    y = PooledDataArray([1, 9, 3, 2, 2])
    @assert append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    x = PooledDataArray([9, 9, 8])
    y = [1, 9, 3, 2, 2]
    @assert append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    # convert methods
    for from in (@pdata(ones(5, 5)), @data(ones(5, 5)), ones(5, 5))
        for (to, totype) in ((PooledDataArray{Float32,UInt16,2}, PooledDataArray{Float32,UInt16,2}),
                             (PooledDataArray{Float32,UInt32,2}, PooledDataArray{Float32,UInt32,2}),
                             (PooledDataArray{Float64,UInt16,2}, PooledDataArray{Float64,UInt16,2}),
                             (PooledDataArray{Float64,UInt32,2}, PooledDataArray{Float64,UInt32,2}),
                             (PooledDataArray{Float32,UInt16}, PooledDataArray{Float32,UInt16,2}),
                             (PooledDataArray{Float32,UInt32}, PooledDataArray{Float32,UInt32,2}),
                             (PooledDataArray{Float64,UInt16}, PooledDataArray{Float64,UInt16,2}),
                             (PooledDataArray{Float64,UInt32}, PooledDataArray{Float64,UInt32,2}),
                             (PooledDataArray{Float32}, PooledDataArray{Float32,UInt32,2}),
                             (PooledDataArray{Float64}, PooledDataArray{Float64,UInt32,2}),
                             (PooledDataArray, PooledDataArray{Float64,UInt32,2}))
            rettype = typeof(convert(to, from))
            if rettype != totype
                error("convert($to, ::$(typeof(from))) returned $rettype (expected $totype)")
            end
        end
    end

    da = convert(DataArray, @pdata(ones(5, 5)))
    @assert isequal(da, @data(ones(5, 5)))
    @assert typeof(da) == DataArray{Float64,2}
    da = convert(DataArray{Float32}, @pdata(ones(5, 5)))
    @assert isequal(da, @data(ones(Float32, 5, 5)))
    @assert typeof(da) == DataArray{Float32,2}
    da = convert(DataArray{Float32,2}, @pdata(ones(5, 5)))
    @assert isequal(da, @data(ones(Float32, 5, 5)))
    @assert typeof(da) == DataArray{Float32,2}

    # permute
    pda = @pdata([NA, "A", "B", "C", "A", "B"])
    @test isequal(Base.permute!!(copy(pda), [2, 5, 3, 6, 4, 1]), @pdata(["A", "A", "B", "B", "C", NA]))
    @test isequal(Base.ipermute!!(copy(pda), [6, 1, 3, 5, 2, 4]), @pdata(["A", "A", "B", "B", "C", NA]))
end


text:
# This file defines struct type for keeping simulation data.
# It also defines functions for saving data in MAT (Matlab) file.

# module Results

using DrWatson
using MAT
using Statistics
using Suppressor

# using Attributes

#This type is to store simulation results
struct Result
    #parameters
    n::UInt32 #number of nodes
    g::UInt32 #number of attributes
    gammas::Array{Float64,1} #gammas
    maxtime::Float64 #simulation stopping time
    ode_fun_name::String #used ode_function

    attr_name::String #atribute type
    attr_threshold::Float64 #threshold (if any)
    attr_degeneracy::Int #number of possible values

    #true results
    HB::Array{Float64,1} #number of times HB and condition were achieved divided by number of repetitions
    HB_x::Array{Float64,1} #number of times HB was achieved (no matter the condition) divided by number of repetitions
    HB_attr::Array{Float64,1} #number of times HB in attribute similarity matrix was achieved divided by number of repetitions
    HB_only_weights::Array{Float64,1} #number of times HB was achieved but attribute similarity matrix is not in HB (divided by number of times HB was not achieved in ASM)
    BR::Array{Float64,1} #mean value of ratio of balanced triads
    paradise::Array{Float64,1} #number of times paradise was achieved divided by number of repetitions
    hell::Array{Float64,1} #number of times hell was achieved divided by number of repetitions
    weak_balance_in_complete_graph::Array{Float64,1} #number of times weak SB was achieved divided by number of repetitions
    local_polarization::Array{Float64,1} #local polarization measure
    global_polarization::Array{Float64,1} #number of times system was globally polarized

    initial_neg_links_count::Array{Float64,1} #mean number of initial negative links
    initial_neg_links_count_std::Array{Float64,1} #std of initial negative links
    links_destab_changed::Array{Float64,2} #array of mean numbers of positive links destabilized, negative links destab., positive links changes and neg links changed
    links_destab_changed_std::Array{Float64,2} #array of std of positive links destabilized, negative links destab., positive links changes and neg links changed
    Deltas::Array{Float64,2} # array with mean numbers of triad counts
    Deltas_std::Array{Float64,2} # std of triad counts

    sim::Array{Float64,1} #mean similarity of layers at the end of simulations
    x_attr_sim::Array{Float64,1} #mean similarity between weights and attributes at the end of simulations
    stab::Array{Float64,1} #mean stability at the end of simulations
    times::Array{Float64,1} #mean times to reach the end of simulation

    BR_std::Array{Float64,1} #std deviation of ratio of balanced triads
    local_polarization_std::Array{Float64,1} #std deviation of local polarization measure
    sim_std::Array{Float64,1} #std deviation of similarity of layers at the end of simulations
    x_attr_sim_std::Array{Float64,1} #std deviation of similarity of weights and attributes at the end of simulations
    times_std::Array{Float64,1} #std deviation of times to reach the end of simulation

    zmax::Array{UInt32,1} #number of repetitions

    data::Array{Float64,2} #all data listed in an array
    interpretation::Array{String,1} #interpretation of columns of above field
    # interpretation = ["N", "G", "gamma", "HB_in_attr",
    # "zmax", "HB", "HB_x", "paradise", "hell", "weak_balance_in_complete_graph",
    # "sim", "sim_std",
    # "x_attr_sim", "x_attr_sim_std",
    # "BR", "BR_std", "LP", "LP_std", "GP", "stab",
    # "times", "times_std",
    # "pos_links_destab", "pos_links_destab_std", "neg_links_destab", "neg_links_destab_std",
    # "pos_links_changed", "pos_links_changed_std", "neg_links_changed", "neg_links_changed_std",
    # "initial_neg_links_count", "initial_neg_links_count_std",
    # "Delta0", "Delta0_std", "Delta1", "Delta1_std", "Delta2",
    # "Delta2_std", "Delta3", "Delta3_std"]

    # Following is the different way of keeping and later saving data.
    # Previously (above) data was separated into AL-balanced and AL-unbalanced cases. Below it is not.
    data2::Array{Any,2} #all data listed in an array
    interpretation2::Array{String,1} #interpretation of columns of above field
    # interpretation = ["N", "G", "gamma", "maxtime", "ode_fun_name",
    # "attr", "attr_name", "attr_threshold",  "attr_degeneracy",
    # "zmax", "HB", "HB_x", "paradise", "hell", "weak_balance_in_complete_graph",
    # "sim", "sim_std",
    # "x_attr_sim", "x_attr_sim_std",
    # "BR", "BR_std", "LP", "LP_std", "GP", "stab",
    # "times", "times_std",
    # "pos_links_destab", "pos_links_destab_std", "neg_links_destab", "neg_links_destab_std",
    # "pos_links_changed", "pos_links_changed_std", "neg_links_changed", "neg_links_changed_std",
    # "initial_neg_links_count", "initial_neg_links_count_std",
    # "Delta0", "Delta0_std", "Delta1", "Delta1_std", "Delta2",
    # "Delta2_std", "Delta3", "Delta3_std"]
end
export Result

# Constructor
function Result(
    n::Int,
    attr::AbstractAttributes,
    gammas::Vector{Float64},
    maxtime::Float64,
    ode_fun_name::String,
)
    nb = length(gammas)

    HB = zeros(Float64, nb)
    HB_x = copy(HB)
    HB_attr = copy(HB)
    HB_only_weights = copy(HB)
    BR = copy(HB)
    paradise = copy(HB)
    hell = copy(HB)
    weak_balance_in_complete_graph = copy(HB)
    local_polarization = copy(HB)
    global_polarization = copy(HB)

    initial_neg_links_count = copy(HB)
    initial_neg_links_count_std = copy(HB)
    Deltas = zeros(Float64, 4, nb)
    links_destab_changed = copy(Deltas)

    sim = zeros(Float64, nb)
    x_attr_sim = copy(sim)
    stab = copy(sim)
    times = copy(sim)

    BR_std = copy(sim)
    local_polarization_std = copy(sim)
    links_destab_changed_std = copy(Deltas)
    Deltas_std = copy(Deltas)
    sim_std = copy(sim)
    x_attr_sim_std = copy(sim)
    times_std = copy(sim)

    zmax = zeros(UInt32, nb)

    interpretation = [
        "N",
        "G",
        "gamma",
        "HB_in_attr",
        "zmax",
        "HB",
        "HB_x",
        "paradise",
        "hell",
        "weak_balance_in_complete_graph",
        "sim",
        "sim_std",
        "x_attr_sim",
        "x_attr_sim_std",
        "BR",
        "BR_std",
        "LP",
        "LP_std",
        "GP",
        "stab",
        "times",
        "times_std",
        "pos_links_destab",
        "pos_links_destab_std",
        "neg_links_destab",
        "neg_links_destab_std",
        "pos_links_changed",
        "pos_links_changed_std",
        "neg_links_changed",
        "neg_links_changed_std",
        "initial_neg_links_count",
        "initial_neg_links_count_std",
        "Delta0",
        "Delta0_std",
        "Delta1",
        "Delta1_std",
        "Delta2",
        "Delta2_std",
        "Delta3",
        "Delta3_std",
    ]
    data = zeros(Float64, nb * 2, length(interpretation))

    interpretation2 = [
        "N",
        "G",
        "gamma",
        "maxtime",
        "ode_fun_name",
        "attr_name",
        "attr_threshold",
        "attr_degeneracy",
        "zmax",
        "HB",
        "HB_x",
        "paradise",
        "hell",
        "weak_balance_in_complete_graph",
        "sim",
        "sim_std",
        "x_attr_sim",
        "x_attr_sim_std",
        "BR",
        "BR_std",
        "LP",
        "LP_std",
        "GP",
        "stab",
        "times",
        "times_std",
        "pos_links_destab",
        "pos_links_destab_std",
        "neg_links_destab",
        "neg_links_destab_std",
        "pos_links_changed",
        "pos_links_changed_std",
        "neg_links_changed",
        "neg_links_changed_std",
        "initial_neg_links_count",
        "initial_neg_links_count_std",
        "Delta0",
        "Delta0_std",
        "Delta1",
        "Delta1_std",
        "Delta2",
        "Delta2_std",
        "Delta3",
        "Delta3_std",
    ]
    data2 = Array{Any,2}(undef, nb, length(interpretation2))
    data2[:, 1:4] .= 0
    data2[:, 5:6] .= ""
    data2[:, 7:end] .= 0

    Result(
        n,
        attr.g,
        gammas,
        maxtime,
        ode_fun_name,
        get_name(attr),
        get_threshold(attr),
        get_degeneracy(attr),
        HB,
        HB_x,
        HB_attr,
        HB_only_weights,
        BR,
        paradise,
        hell,
        weak_balance_in_complete_graph,
        local_polarization,
        global_polarization,
        initial_neg_links_count,
        initial_neg_links_count_std,
        links_destab_changed,
        links_destab_changed_std,
        Deltas,
        Deltas_std,
        sim,
        x_attr_sim,
        stab,
        times,
        BR_std,
        local_polarization_std,
        sim_std,
        x_attr_sim_std,
        times_std,
        zmax,
        data,
        interpretation,
        data2,
        interpretation2,
    )
end
export Result

function Result(
    n::Int,
    g::Int,
    gammas::Vector{Float64},
    maxtime::Float64,
    ode_fun_name::String,
)
    Result(n, BinaryAttributes(g), gammas, maxtime, ode_fun_name)
end
export Result

function Result(n::Int, g::Int, gammas::Vector{Float64}, maxtime::Float64)
    Result(n, g, gammas, maxtime, "Heider5!")
end
export Result

function Result(n::Int, g::Int, gammas::Vector{Float64}, ode_fun_name::String)
    Result(n, g, gammas, 1000.0, ode_fun_name)
end
export Result

function Result(n::Int, g::Int, gammas::Vector{Float64})
    Result(n, g, gammas, 1000.0)
end
export Result

#save results to MATLAB or jld2 file
function save_result(res, filename; ext = "mat")
    if ext == "mat"
        if !endswith(filename, "mat")
            filename = filename * ".mat"
        end
        file = matopen(filename, "w")
        write(file, "result", res)
        close(file)
    elseif ext == "jld2"
        if endswith(filename, r"\.[0-9a-z]+$")
            filename = replace(filename, (r"\.[0-9a-z]+$" => ""))
        end
        for (i, gamma) in enumerate(res.gammas)
            # This was not yet simulated. End of saving
            if res.zmax[i] == 0
                break
            end

            #creating filename
            local filename_g
            @suppress filename_g = savename(filename, (gamma = gamma,))
            filename_g2 = filename_g * "_1.jld2"

            # data_row = res.data[2*i-1, :]
            data_row = res.data2[i, :]
            @suppress @tagsave(
                filename_g2,
                # Dict(res.interpretation .=> data_row)
                Dict(res.interpretation2 .=> data_row)
            )
            # filename_g2 = filename_g*"_2.jld2"
            # data_row = res.data[2*i, :]
            # @suppress @tagsave(
            #     filename_g2,
            #     Dict(res.interpretation .=> data_row)
            # )
        end
    end

end
export save_result

function read_result(filename, varname)
    if !endswith(filename, ".mat")
        filename *= ".mat"
    end
    file = matopen(filename)
    out = read(file, varname) # note that this does NOT introduce a variable ``varname`` into scope
    close(file)
    out
end
export read_result

#function to update result fields based on elements of fields tuple
function update_result!(res::Result, fields)
    HB,
    HB_x,
    HB_attr,
    sim,
    x_attr_sim,
    BR,
    paradise,
    hell,
    initial_neg_links_count,
    links_destab_changed,
    Deltas,
    weak_balance_in_complete_graph,
    local_polarization,
    global_polarization,
    stab,
    times,
    i,
    rep,
    firstline = fields
    res.HB[i] = sum(HB) / rep
    res.HB_x[i] = sum(HB_x) / rep
    res.HB_attr[i] = sum(HB_attr) / rep
    res.HB_only_weights[i] = sum(HB[HB_attr.==0]) / ((1 - res.HB_attr[i]) * rep)
    res.BR[i] = sum(BR) / rep
    res.paradise[i] = sum(paradise) / rep
    res.hell[i] = sum(hell) / rep
    res.weak_balance_in_complete_graph[i] = sum(weak_balance_in_complete_graph) / rep
    res.local_polarization[i] = sum(local_polarization) / rep
    res.global_polarization[i] = sum(global_polarization) / rep

    res.initial_neg_links_count[i] = sum(initial_neg_links_count) / rep
    res.links_destab_changed[:, i] = sum(links_destab_changed, dims = 2) / rep
    res.Deltas[:, i] = sum(Deltas, dims = 2) / rep

    res.sim[i] = sum(sim) / rep
    res.x_attr_sim[i] = sum(x_attr_sim) / rep
    res.stab[i] = sum(stab) / rep
    res.times[i] = sum(times) / rep

    res.BR_std[i] = std(BR[1:rep])
    res.local_polarization_std[i] = std(local_polarization[1:rep])
    res.initial_neg_links_count_std[i] = std(initial_neg_links_count[1:rep])
    res.Deltas_std[:, i] = std(Deltas[:, 1:rep], dims = 2)
    res.links_destab_changed_std[:, i] = std(links_destab_changed[:, 1:rep], dims = 2)

    res.sim_std[i] = std(sim[1:rep])
    res.x_attr_sim_std[i] = std(x_attr_sim[1:rep])
    res.times_std[i] = std(times[1:rep])

    res.zmax[i] = rep

    #line when HB is in attribute matrix
    mask_true = (HB_attr .== 1)
    smt = sum(mask_true)
    mask_false = .!mask_true
    smf = rep - smt
    # try
    res.data[firstline, :] = [
        res.n,
        res.g,
        res.gammas[i],
        true,
        smt,
        sum(HB[mask_true]) / smt,
        sum(HB_x[mask_true]) / smt,
        sum(paradise[mask_true]) / smt,
        sum(hell[mask_true]) / smt,
        sum(weak_balance_in_complete_graph[mask_true]) / smt,
        sum(sim[mask_true]) / smt,
        std(sim[1:rep][mask_true[1:rep]]),
        sum(x_attr_sim[mask_true]) / smt,
        std(x_attr_sim[1:rep][mask_true[1:rep]]),
        sum(BR[mask_true]) / smt,
        std(BR[1:rep][mask_true[1:rep]]),
        sum(local_polarization[mask_true]) / smt,
        std(local_polarization[1:rep][mask_true[1:rep]]),
        sum(global_polarization[mask_true]) / smt,
        sum(stab[mask_true]) / smt,
        sum(times[mask_true]) / smt,
        std(times[1:rep][mask_true[1:rep]]),
        sum(links_destab_changed[1, mask_true]) / smt,
        std(links_destab_changed[1, 1:rep][mask_true[1:rep]]),
        sum(links_destab_changed[2, mask_true]) / smt,
        std(links_destab_changed[2, 1:rep][mask_true[1:rep]]),
        sum(links_destab_changed[3, mask_true]) / smt,
        std(links_destab_changed[3, 1:rep][mask_true[1:rep]]),
        sum(links_destab_changed[4, mask_true]) / smt,
        std(links_destab_changed[4, 1:rep][mask_true[1:rep]]),
        sum(initial_neg_links_count[mask_true]) / smt,
        std(initial_neg_links_count[1:rep][mask_true[1:rep]]),
        sum(Deltas[1, mask_true]) / smt,
        std(Deltas[1, 1:rep][mask_true[1:rep]]),
        sum(Deltas[2, mask_true]) / smt,
        std(Deltas[2, 1:rep][mask_true[1:rep]]),
        sum(Deltas[3, mask_true]) / smt,
        std(Deltas[3, 1:rep][mask_true[1:rep]]),
        sum(Deltas[4, mask_true]) / smt,
        std(Deltas[4, 1:rep][mask_true[1:rep]]),
    ]
    # catch y
    #     display([firstline, rep])
    #     display([length(mask_true), sum(mask_true), smt])
    #     display([length(HB_x), sum(HB_x[mask_true]), smt])
    #     display([length(sim), sum(sim[mask_true]), smt])
    #     display([length(HB_x), sum(HB_x[mask_true]), smt])
    #     display([length(sim), sum(sim[mask_true]), smt])
    #     throw(y)
    # end
    #line when HB is not in attribute matrix
    res.data[firstline+1, :] = [
        res.n,
        res.g,
        res.gammas[i],
        false,
        smf,
        sum(HB[mask_false]) / smf,
        sum(HB_x[mask_false]) / smf,
        sum(paradise[mask_false]) / smf,
        sum(hell[mask_false]) / smf,
        sum(weak_balance_in_complete_graph[mask_false]) / smf,
        sum(sim[mask_false]) / smf,
        std(sim[1:rep][mask_false[1:rep]]),
        sum(x_attr_sim[mask_false]) / smf,
        std(x_attr_sim[1:rep][mask_false[1:rep]]),
        sum(BR[mask_false]) / smf,
        std(BR[1:rep][mask_false[1:rep]]),
        sum(local_polarization[mask_false]) / smf,
        std(local_polarization[1:rep][mask_false[1:rep]]),
        sum(global_polarization[mask_false]) / smf,
        sum(stab[mask_false]) / smf,
        sum(times[mask_false]) / smf,
        std(times[1:rep][mask_false[1:rep]]),
        sum(links_destab_changed[1, mask_false]) / smf,
        std(links_destab_changed[1, 1:rep][mask_false[1:rep]]),
        sum(links_destab_changed[2, mask_false]) / smf,
        std(links_destab_changed[2, 1:rep][mask_false[1:rep]]),
        sum(links_destab_changed[3, mask_false]) / smf,
        std(links_destab_changed[3, 1:rep][mask_false[1:rep]]),
        sum(links_destab_changed[4, mask_false]) / smf,
        std(links_destab_changed[4, 1:rep][mask_false[1:rep]]),
        sum(initial_neg_links_count[mask_false]) / smf,
        std(initial_neg_links_count[1:rep][mask_false[1:rep]]),
        sum(Deltas[1, mask_false]) / smf,
        std(Deltas[1, 1:rep][mask_false[1:rep]]),
        sum(Deltas[2, mask_false]) / smf,
        std(Deltas[2, 1:rep][mask_false[1:rep]]),
        sum(Deltas[3, mask_false]) / smf,
        std(Deltas[3, 1:rep][mask_false[1:rep]]),
        sum(Deltas[4, mask_false]) / smf,
        std(Deltas[4, 1:rep][mask_false[1:rep]]),
    ]

    # "maxtime", "ode_fun_name",
    # "attr", "attr_name", "attr_threshold",  "attr_degeneracy",
    res.data2[i, :] = [
        res.n,
        res.g,
        res.gammas[i],
        res.maxtime,
        res.ode_fun_name,
        res.attr_name,
        res.attr_threshold,
        res.attr_degeneracy,
        rep,
        sum(HB) / rep,
        sum(HB_x) / rep,
        sum(paradise) / rep,
        sum(hell) / rep,
        sum(weak_balance_in_complete_graph) / rep,
        sum(sim) / rep,
        std(sim[1:rep]),
        sum(x_attr_sim) / rep,
        std(x_attr_sim[1:rep]),
        sum(BR) / rep,
        std(BR[1:rep]),
        sum(local_polarization) / rep,
        std(local_polarization[1:rep]),
        sum(global_polarization) / rep,
        sum(stab) / rep,
        sum(times) / rep,
        std(times[1:rep]),
        sum(links_destab_changed[1, mask_false]) / rep,
        std(links_destab_changed[1, 1:rep]),
        sum(links_destab_changed[2, mask_false]) / rep,
        std(links_destab_changed[2, 1:rep]),
        sum(links_destab_changed[3, mask_false]) / rep,
        std(links_destab_changed[3, 1:rep]),
        sum(links_destab_changed[4, mask_false]) / rep,
        std(links_destab_changed[4, 1:rep]),
        sum(initial_neg_links_count) / rep,
        std(initial_neg_links_count[1:rep]),
        sum(Deltas[1, mask_false]) / rep,
        std(Deltas[1, 1:rep]),
        sum(Deltas[2, mask_false]) / rep,
        std(Deltas[2, 1:rep]),
        sum(Deltas[3, mask_false]) / rep,
        std(Deltas[3, 1:rep]),
        sum(Deltas[4, mask_false]) / rep,
        std(Deltas[4, 1:rep]),
    ]
    # "zmax", "HB", "HB1", "HB2", "HB12", "sim", "sim_std", "stab",
    # "times", "times_std"]
end
export update_result!

DrWatson.allaccess(::Result) = (:attr_name, :n, :g, :attr_degeneracy)
# end

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
# Combines two 1-output circuits c1 and c2 with the same paramterers and number of gates into a 2-output circuit
# The output values of the result is the 2-output goal whose components are the output values of c1 and of c2.
function combine_circuits( c1::Chromosome, c2::Chromosome )
  # Two functions used only within this function
  inp1(i,numinputs) = i <= numinputs ? i : 2*(i-numinputs)+numinputs-1
  inp2(i,numinputs) = i <= numinputs ? i : 2*(i-numinputs)+numinputs
  @assert c1.params.numinputs == c2.params.numinputs
  @assert c1.params.numoutputs == 1
  @assert c2.params.numoutputs == 1
  @assert c1.params.numinteriors == c2.params.numinteriors
  @assert c1.params.numlevelsback == c2.params.numlevelsback
  if length(c1.interiors) < length(c2.interiors)
    println("  c1ints: ",length(c1.interiors),"  c2ints: ",length(c2.interiors))
    ctemp = c1
    c1 =c2
    c2 = ctemp
  end
  #println("c1ints: ",length(c1.interiors),"  c2ints: ",length(c2.interiors))
  p = Parameters( c1.params.numinputs, 2, c1.params.numinteriors+c2.params.numinteriors,
    c1.params.numlevelsback+c2.params.numlevelsback )
  inputs = c1.inputs
  interiors = InteriorNode[]
  new_numinteriors = 2*length(c2.interiors)
  #println("new_numints: ",new_numinteriors)
  #print_parameters(p)
  for i = 1:length(c2.interiors)
    int_node1 = deepcopy(c1.interiors[i])
    input1 = inp1(c1.interiors[i].inputs[1],p.numinputs)
    input2 = inp1(c1.interiors[i].inputs[2],p.numinputs)
    int_node1.inputs = [input1,input2]
    #println("i: ",i,"  int_node1: ",int_node1)
    push!(interiors,int_node1)
    int_node2 = deepcopy(c2.interiors[i])
    input1 = inp2(c2.interiors[i].inputs[1],p.numinputs)
    input2 = inp2(c2.interiors[i].inputs[2],p.numinputs)
    int_node2.inputs = [input1,input2]
    #println("i: ",i,"  int_node2: ",int_node2)
    push!(interiors,int_node2)
  end
  #println("interiors: ",interiors)
  outputs = [OutputNode(p.numinputs+new_numinteriors-1),OutputNode(p.numinputs+new_numinteriors)]
  Chromosome( p, c1.inputs, interiors, outputs, 0.0, 0.0 )
end

# Starting with circuit c (and its parameters), try to evolve a circuit whose output is the same as c
#   and whose number of active nodes is reduced.
# Version of 5/7/21
function evolve_reduce_numactive( c::Chromosome, maxsteps::Int64; num_mutations::Int64=1, funcs::Vector{Func}=default_funcs(c.params.numinputs)  )
  p = c.params
  goal = output_values(c)
  current_num_active = number_active_gates( c )
  for i = 1:maxsteps
    new_c = c
    for _ = 1:num_mutations
      new_c = deepcopy(c)
      mutate_chromosome!( new_c, funcs )
    end
    if output_values(new_c) == goal && number_active_gates( new_c ) <= current_num_active
      c = new_c
      if number_active_gates( c ) < current_num_active
        #println("number active reduced from ",current_num_active," to ",number_active_gates( new_c )," on step ",i,".")
        current_num_active = number_active_gates( c )
      end
    end
  end
  (c, current_num_active )
end

# Use neutral evolution to find a chromosome c that maximizes funct(c)
# funct( c::Chromosome ) returns a Float64.
function evolve_function( funct, p::Parameters, funcs::Vector{Func}, max_steps::Int64;
    goallist::Vector{Vector{MyInt}}=[MyInt[]], max_evolve_steps::Int64=10000 )
  println("goallist: ",goallist)
  use_goal = !(goallist==[MyInt[]])
  c = random_chromosome( p, funcs )
  orig_c = deepcopy(c)
  goal = goallist[1]
  println("goal: ",goal)
  if use_goal
    result = mut_evolve( c, goallist, funcs, max_evolve_steps, print_improvements = true )
    steps = result[2]
    if steps == max_evolve_steps
      error("evolution to goal failed in function evolve_function()")
    else
      println("evolution to goal suceeded in function evolve_function()")
    end
    c = result[1]
    goal = output_values(c)
    orig_c = deepcopy(c)
  end
  current_fitness = funct(c)
  println("starting fitness: ",current_fitness)
  for i = 1:max_steps
    prev_c = deepcopy(c)
    (c,active) = mutate_chromosome!(c,funcs)
    if use_goal && output_values(c) != goal
      c = prev_c
      continue
    end
    new_fitness = funct(c)
    if new_fitness < current_fitness
      c = prev_c
    elseif new_fitness > current_fitness
      println("i: ",i,"  new_fitness: ", new_fitness,"  funct(c): ",funct(c) )
      current_fitness = new_fitness
    end
    println("i: ",i,"  current_fitness: ", current_fitness, "  funct(c): ",funct(c)  )
  end
  (c,orig_c)
end

# to test circuit_evolve()
function test_evolve()
  p = Parameters(numinputs=2,numoutputs=1,numinteriors=3,numlevelsback=3);
  funcs=default_funcs(p.numinputs);c = random_chromosome(p,funcs); goal=output_values(c)
  res = mut_evolve(c,[goal],funcs,1000); c_dest=res[1]
  circuit_evolve(c,c_dest, 1000 )
end

# Attempts to reduce the number active of the chromosome c by mutational neutral evolution
# Older version written prior to 7/18/20.
# The chomosome c is assumed to be at optimal fitness, and this function will terminate with
#   an error if it is not optimal.
# Then mutation is done nreps times, and the mutation is accepted if number_active(c) is reduced.
function mut_reduce_numactive( c::Chromosome, goallist::GoalList, funcs::Vector{Func}, nreps::Integer;
      hamming_sel::Bool=true, num_mutations::Int64=1, print_improvements::Bool=false )
  output = output_values(c)   # Executes c if it has not already been executed
  goals_matched = hamming_sel ? goals_matched_hamming : goals_matched_exact
  ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs )
  if sort(output) != sort(goallist[matched_goals[1]])
    error("input chromosome must be have optimal fitness in function mut_reduce_numactive()")
  end
  best_num_active = number_active(c)
  #println("init num_active  num_active: ",best_num_active)
  orig_c = deepcopy(c)
  #println("initial fitness: ",fitness)
  step = 0
  worse = 0
  same = 0
  better = 0
  prev_c = deepcopy(c)  # previous generation c
  best_c = c
  @assert output_values(prev_c) == output_values(orig_c)
  (c, matched_goals, matched_goals_list ) = next_chromosome!(c, goallist, funcs,
      hamming_sel=hamming_sel, use_robustness=false, num_mutations=num_mutations )
  new_num_active = number_active(c)
  #println("before while new_num_active: ",new_num_active)
  #while step < nreps && new_fitness < c.params.numoutputs
  while step < nreps
    #println("step: ",step,"  output: ",output,"   ")
    if c.fitness > fitness  # This should never happen
      error("c.fitness > fitness")
      if print_improvements
        matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
        println("fitness improved from ",fitness," to ",c.fitness," at step: ",step, "  goals_matched: ", matched_goals, "  matched goal components: ", matched_goal_components)
      end
      fitness = c.fitness
      worse = 0
      same = 0
      better += 1
    elseif c.fitness == fitness && new_num_active < best_num_active
      #println("new chromosome with fitness: ",fitness, "  and num_active: ", new_num_active )
      if print_improvements
        matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
        println("num_active decreased from ",best_num_active," to ",new_num_active," at step: ",step,"  goals_matched: ", matched_goals, "  matched goal componets: ", matched_goal_components)
      end
      worse = 0
      same = 0
      better += 1
      best_c = deepcopy(c)
      @assert new_num_active == number_active(c)
      best_num_active = new_num_active
      #println("num_active_improved  new_num_active: ",new_num_active)
    elseif c.fitness == fitness && new_num_active == best_num_active
      #println("new chromosome with fitness: ",fitness, "  and num_active: ", new_num_active )
      same += 1
    elseif c.fitness == fitness && new_num_active > best_num_active
      #println("discarded chromosome with fitness: ",fitness, " and num_active: ", new_num_active )
      c = prev_c
      worse += 1
    elseif c.fitness < fitness
      #println("discarded chromosome with new_fitness: ",c.fitness)
      c = prev_c
      worse += 1
    end
    step += 1
    prev_c = deepcopy(c)
    if step < nreps
      (c, matched_goals, matched_goals_list ) = next_chromosome!(c, goallist, funcs,
            hamming_sel=hamming_sel, use_robustness=false, num_mutations=num_mutations )
      #output = output_values(c)   # Executes c if it has not already been executed
      new_num_active = number_active(c)
      #println("end while new_num_active: ",new_num_active)
    end
  end
  if step == nreps
    #println("mut_reduce_numactive finished at step limit ",nreps," with fitness: ", best_c.fitness, " and num_active: ", best_num_active )
  else
    error("step == nreps")
    println("mut_reduce_numactive finished in ",step," steps with fitness: ", c.fitness, " and num_active: ", new_num_active )
  end
  if orig_c.fitness > c.fitness
    error("orig_c.fitness > c.fitness")
    #println("(orig_c.fitness,c.fitness): ",(orig_c.fitness,c.fitness))
    output = output_values(orig_c)
    ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs )
    return(orig_c,step,0,0,0,output,matched_goals,matched_goals_list)
  end
  ##println("worse: ",worse,"  same: ",same,"  better: ",better)
  #println("matched_goals: ",matched_goals)
  if step < nreps
    if print_improvements
      matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
      println("num_active decreased from ",num_active," to ",new_num_active," at step: ",step, "  goals_matched: ", matched_goals, "  matched goal componets: ", matched_goal_components)
    end
    @assert sort(output) == sort(goallist[matched_goals[1]])
  end
  ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, best_c.params.numinputs )
  #println("assert  best_num_active: ",best_num_active,"  number_active(best_c): ",number_active(best_c))
  @assert best_num_active == number_active(best_c)
  (best_c,step,worse,same,better,output,matched_goals,matched_goals_list)
end

# Evaluate the performance of function evolve_reduce_numactive() which is in evolve_funct.jl.
# n_repeats is the maximum number of tries in mut_evolve_repeat()
# For each pair of chromosomes
#   Creates two chromosomes using parameter settings p.
#   Combines them into one chromosome with one more input
#   For each trial run evolve_reduce_numactive() and record average and minimum numactive.
function test_reduce_numactive( p::Parameters, numpairs::Int64, numtrials::Int64, n_repeats::Int64, maxsteps::Int64;
      num_mutations::Int64=1, csvfile::String="" )
  df = DataFrame()
  df.combined_goal = Goal[]
  df.numinputs = Int64[]
  df.numgates = Int64[]
  df.numlevelsback = Int64[]
  df.numpairs = Int64[]
  df.numtrials = Int64[]
  df.maxsteps = Int64[]
  df.nactive0 = Float64[]
  df.nactive1_min = Float64[]
  df.nactive1_mean = Float64[]
  df.nactive1_std = Float64[]
  if num_mutations >= 2
    df.nactive2_min = Float64[]
    df.nactive2_mean = Float64[]
    df.nactive2_std = Float64[]
  end
  if num_mutations == 3
    df.nactive3_min = Float64[]
    df.nactive3_mean = Float64[]
    df.nactive3_std = Float64[]
  end
  funcs = default_funcs(p.numinputs)
  for npair = 1:numpairs
    #c1 = random_chromosome(p); print(output_values(c1),"  "); print_circuit(c1)
    #c2 = random_chromosome(p); print(output_values(c2),"  "); print_circuit(c2)
    goal1 = randgoal(p)
    c = random_chromosome(p)
    (nc1,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve_repeat(n_repeats, p, [goal1], funcs, maxsteps )
    goal2 = randgoal(p)
    (nc2,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve_repeat(n_repeats, p, [goal2], funcs, maxsteps )
    c = combine_circuits(nc1,nc2)
    rowlist = [ output_values(c), p.numinputs, p.numinteriors, p.numlevelsback, numpairs, numtrials, maxsteps ]
    #@assert [output_values(nc1)[1],output_values(nc1)[1]] == output_values(c)
    println("output_values(c): ",output_values(c))
    println("goal: ",[output_values(nc1)[1],output_values(nc2)[1]])
    #nactive_lists = map(i->reduce_nactive_helper(c,maxsteps,num_mutations=num_mutations), collect(1:numtrials ) )
    nactive_lists = pmap(i->reduce_nactive_helper(c,maxsteps,num_mutations=num_mutations), collect(1:numtrials ) )
    nactive0 = vcat( [ na[1] for na in nactive_lists]... )
    nactive1 = vcat( [ na[2] for na in nactive_lists]... )
    if num_mutations == 1
      row_results = [mean(nactive0),minimum(nactive1),mean(nactive1),std(nactive1)]
    end
    if num_mutations >= 2
      nactive2 = vcat( [ na[3] for na in nactive_lists]... )
      row_results = [mean(nactive0),minimum(nactive2),mean(nactive1),std(nactive1),minimum(nactive2),mean(nactive2),std(nactive2)]
    end
    if num_mutations >= 3
      nactive3 = vcat( [ na[4] for na in nactive_lists]... )
      row_results = [mean(nactive0),minimum(nactive2),mean(nactive1),std(nactive1),minimum(nactive2),mean(nactive2),std(nactive2),minimum(nactive3),mean(nactive3),std(nactive3)]
    end
    row = deepcopy(rowlist)
    println("length(row): ",length(row),"  length(row_results): ",length(row_results),"  size(df): ",size(df))
    #println("npair: ",npair,"  mean(nactive1): ",mean(nactive1),"  mean(nactive2): ",mean(nactive2),"  mean(nactive3): ",mean(nactive3))
    push!(df,vcat(row,row_results))
  end
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# numpairs: ", numpairs)
      println(f,"# numtrials: ", numtrials)
      println(f,"# maxsteps: ", maxsteps)
      println(f,"# n_repeats: ", n_repeats)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end

# nactive0 is number_active_gates() after evolution to find circuit that maps to goal
#    but before running evolve_reduce_numactive()
function reduce_nactive_helper( c::Chromosome, maxsteps::Int64; num_mutations::Int64=1 )
  #println("num_mutations: ",num_mutations)
  goal = output_values( c )
  (ec,steps) = neutral_evolution( random_chromosome(c.params), goal, maxsteps )
  nactive0 = number_active_gates(ec)
  #println("nactive0: ",nactive0)
  nactive1 = Int64[]
  nc = deepcopy(c)
  (ec,cur_nactive) = evolve_reduce_numactive( nc, maxsteps, num_mutations=1 )
  @assert cur_nactive == number_active_gates(ec)
  push!(nactive1,number_active_gates(ec))
  if num_mutations == 1
    return(nactive0, nactive1)
  end
  nactive2 = Int64[]
  nc = deepcopy(c)
  (ec,cur_nactive) = evolve_reduce_numactive( nc, maxsteps, num_mutations=2 )
  push!(nactive2,number_active_gates(ec))
  if num_mutations == 2
    return(nactive0, nactive1, nactive2)
  end
  nactive3 = Int64[]
  nc = deepcopy(c)
  (ec,cur_nactive) = evolve_reduce_numactive( nc, maxsteps, num_mutations=3 )
  push!(nactive3,number_active_gates(ec))
  #println("trial: ",trial,"  mean(nactive1): ",mean(nactive1),"  mean(nactive2): ",mean(nactive2),"  mean(nactive2): ",mean(nactive2))
  (nactive0, nactive1, nactive2, nactive3)
end

# folder_path is relative to data folder
#function compare_dataframe( folder_path::String, suffix::String; csvfile::String="" )
function compare_dataframe( folder_path::String, suffix::String )
  data_prefix = "../data/"
  rdf = read_dataframe("$(data_prefix)$(folder_path)/reduce_numactive_$(folder_path)_$(suffix).csv")
  kdf = read_dataframe("$(data_prefix)$(folder_path)/k_complexity$(folder_path)$(suffix).csv")
  cdf = DataFrame()
  cdf.goal = kdf.goal
  cdf.kng = kdf.num_gates
  cdf.rng = rdf.nactive1_min
  cdf.k_minus_r_diff = cdf.kng .- cdf.rng
  csvfile = "$(data_prefix)$(folder_path)/compare_$(folder_path)_$(suffix).csv"
  lines = String[]
  open( "$(data_prefix)$(folder_path)/reduce_numactive_$(folder_path)_$(suffix).csv", "r" ) do infile
    line = readline(infile)
    while line[1] == '#'
      #println("line: ",line)
      push!(lines,line)
      line = readline(infile)
    end
  end
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      for line in lines
        println(f, line)
      end
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  cdf
end

# goal1 is the first uncombined goal
# goal2 is the first uncombined goal
# p is the parameters for these goals
function compare_two_vs_one_output( p::Parameters, numtrials::Int64, maxsteps::Int64 )
  @assert p.numoutputs == 1
  goal1 = randgoal(p)
  goal2 = randgoal(p)
  pc = Parameters( p.numinputs+1, 1, 2*p.numinteriors, 2*p.numlevelsback )
  pb = Parameters( p.numinputs, 2, 2*p.numinteriors, 2*p.numlevelsback )
  println("pc: ",pc)
  numbits = 2^p.numinputs
  goal = [(deepcopy(goal1)[1] <<= numbits) | goal2[1]]
  println("[goal]: ",[goal],"  [goal1,goal2]: ",[goal1[1],goal2[1]])
  csteps_sum = 0
  c_successes = 0
  bsteps_sum = 0
  b_successes = 0
  for i = 1:numtrials
    cr = random_chromosome(pc)
    #println("[goal]: ",[goal])
    (cnc,csteps) = neutral_evolution( cr, goal, maxsteps )
    if csteps < maxsteps
      csteps_sum += csteps
      c_successes += 1
    end
    cb = random_chromosome(pb)
    #println("[goal1,goal2]: ",[goal1[1],goal2[1]])
    (bnc,bsteps) = neutral_evolution( cb, [goal1[1],goal2[1]], maxsteps )
    if bsteps < maxsteps
      bsteps_sum += bsteps
      b_successes += 1
    end
    #println("(csteps,bsteps): ",(csteps,bsteps))
  end
  println("csteps_avg: ",c_successes > 0 ? csteps_sum/c_successes : maxsteps )
  println("bsteps_avg: ",b_successes > 0 ? bsteps_sum/b_successes : maxsteps )
  (goal, numtrials, c_successes,c_successes > 0 ? csteps_sum/c_successes : maxsteps, b_successes, b_successes > 0 ? bsteps_sum/b_successes : maxsteps )
end

function run_compare_two_vs_one_output( p::Parameters, numruns::Int64, numtrials::Int64, maxsteps::Int64; csvfile::String="" )
  df = DataFrame()
  df.goal = Goal[]
  df.numtrials = Int64[]
  df.c_successes = Int64[]
  df.c_avg = Float64[]
  df.b_successes = Int64[]
  df.b_avg = Float64[]
  rows = pmap( x->compare_two_vs_one_output( p, numtrials, maxsteps ), collect(1:numruns) )
  for r in rows
    push!( df, r )
  end
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# numruns: ",numruns)
      println(f,"# numtrials: ",numtrials)
      println(f,"# maxsteps: ",maxsteps)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end


text:
using SBMLToolkit
using Test

using SBML
using Catalyst
using OrdinaryDiffEq

@testset "SBMLToolkit.jl" begin
    @testset "Model to MTK conversions" begin include("reactionsystem.jl") end
    @testset "Rules" begin include("rules.jl") end
end


text:
using Revise

includet("src/Spring.jl")
IArray = InternalArray

# 1. Solving the equation using the euler method by approximation ------------------------

const timelen = 40

prob = SpringProblem(p0 = 0.8, v0 = 0, m = 1, k = 1.5, b = 0.3, A = 0.4, ω = 2.4)

x = eulermethod_by_approximation(prob, timelen, error=1e-4)

# 2. Plotting the result ------------------------------------------------------------------

timeline = [prob.t0+timelen*i/length(x) for i in 0:length(x)-1]
pos = plot(timeline, IArray(x, 1), label = "movement")
speed = plot(timeline, IArray(x, 2), label = "speed")
plot(pos, speed)

savefig("media/spring_by_approximation.png")

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations