text:
# Electrical Parameters of Machine 1, Milano's Book Page 526
################## Machine Data #####################
Basic = BaseMachine(   0.0, #R
                       0.2995, #Xd_p
                       1.05, #eq_p
                       615.0)  #MVABase

oneDoneQ = OneDOneQMachine(0.0, #R
                        0.8979, #Xd
                        0.646, #Xq
                        0.2995, #Xd_p
                        0.04, #Xq_p
                        7.4, #Td0_p
                        0.033, #Tq0_p
                        615.0)   #MVABase

AndersonFouad = AndersonFouadMachine(0.0, #R
                        0.8979, #Xd
                        0.646, #Xq
                        0.2995, #Xd_p
                        0.646, #Xq_p
                        0.23, #Xd_pp
                        0.4, #Xq_pp
                        7.4, #Td0_p
                        0.01, #Tq0_p #Data not available in Milano: Used 0.01
                        0.03, #Td0_pp
                        0.033, #Tq0_pp
                        615.0) #MVABase

KundurMachine = SimpleFullMachine(0.003, #R on Example 3.1 and 4.1 of Kundur
                                  0.0006, #R_f
                                  0.0284, #R_1d or RD in Machowski
                                  0.0062, #R_1q or RQ on Machowski
                                  1.81, #L_d
                                  1.76, #L_q
                                  1.66, #L_ad or k*M_f or k*M_D in Machowski
                                  1.61, #L_aq or k*M_Q in Machowski
                                  1.66, #L_f1d or L_fD in Machowski. Assumed to be equal to L_ad
                                  1.825, #L_ff
                                  0.1713, #L_1d or L_D in Machowski
                                  0.7525, #L_1q or L_Q in Machowski
                                  555.0) #MVABase

KundurFullMachine = FullMachine(0.003, #R on Example 3.1 and 4.1 of Kundur
                                0.0006, #R_f
                                #0.003, #R_f
                                0.0284, #R_1d or RD in Machowski
                                0.0062, #R_1q or RQ on Machowski
                                1.81, #L_d
                                1.76, #L_q
                                1.66, #L_ad or k*M_f or k*M_D in Machowski
                                1.61, #L_aq or k*M_Q in Machowski
                                1.66, #L_f1d or L_fD in Machowski. Assumed to be equal to L_ad
                                1.825, #L_ff
                                0.1713, #L_1d or L_D in Machowski
                                0.7525, #L_1q or L_Q in Machowski
                                555.0) #MVABase

Mach2_benchmark = OneDOneQMachine(0.0, #R
                        1.3125, #Xd
                        1.2578, #Xq
                        0.1813, #Xd_p
                        0.25, #Xq_p
                        5.89, #Td0_p
                        0.6, #Tq0_p
                        100.0)   #MVABase

Mach3_benchmark = OneDOneQMachine(0.0, #R
                        1.3125, #Xd
                        1.2578, #Xq
                        0.1813, #Xd_p
                        0.25, #Xq_p
                        5.89, #Td0_p
                        0.6, #Tq0_p
                        100.0)   #MVABase

################ Shaft Data #####################
BaseShaft = SingleMass(5.148, #H
                       2.0) #D

FiveShaft = FiveMassShaft(5.148,  #H
                          0.3348, #H_hp
                          0.7306, #H_ip
                          0.8154, #H_lp
                          0.0452, #H_ex,
                          2.0,    #D
                          0.5180, #D_hp
                          0.2240, #D_ip
                          0.2240, #D_lp
                          0.1450, #D_ex
                          0.0518, #D_12
                          0.0224, #D_23
                          0.0224, #D_34
                          0.0145, #D_45
                          33.07,  #K_hp
                          28.59,  #K_ip
                          44.68,  #K_lp
                          21.984) #K_ex

KundurShaft = SingleMass(3.525, #H given in example 3.5 in Kundur
                         2.0) #D: Damping is not given, so I assume the value 2.0

Shaft2_benchmark = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                       0.0) #D

Shaft3_benchmark = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                      0.0) #D

################# PSS Data #####################
no_pss = PSSFixed(0.0)

################ TG Data #####################
fixed_tg = TGFixed(1.0) #eff

typeI_tg = TGTypeI(0.02, #R
                   0.1, #Ts
                   0.45, #Tc
                   0.0, #T3
                   0.0, #T4
                   50.0, #T5
                   0.3, #P_min
                   1.2) #P_max

typeII_tg = TGTypeII(0.05, #R
                     0.3, #T1
                     0.1, #T2
                     1.0, #τ_max
                     0.1) #τ_min

################ AVR Data #####################
proportional_avr = AVRSimple(5000.0) #Kv

fixed_avr = AVRFixed(1.05) #Emf

typeI_avr = AVRTypeI(200.0, #Ka
                     1.0, #Ke
                     0.0012, #Kf
                     0.02, #Ta
                     0.19, #Te
                     1.0, #Tf
                     0.001, #Tr
                     9.9, #Vr_max
                     0.0, #Vr_min
                     0.0006, #Ae
                     0.9)

kundur_avr = AVRFixed(0.00043) #Emf
gen_avr_bench = AVRFixed(1.7613) #Vf

gen2_avr_benchmark_m = AVRTypeI(20.0, #Ka - Gain
                   0.01, #Ke
                   0.063, #Kf
                   0.2, #Ta
                   0.314, #Te
                   0.35, #Tf
                   0.001, #Tr
gen2_avr_benchmark = AVRTypeII(20.0, #K0 - Gain
                     0.2, #T1 - 1st pole
                     0.063, #T2 - 1st zero
                     0.35, #T3 - 2nd pole
                     0.01, #T4 - 2nd zero
                     0.314, #Te - Field current time constant
                     0.001, #Tr - Measurement time constant
                     5.0, #Vrmax
                     -5.0, #Vrmin
                     0.0039, #Ae - 1st ceiling coefficient
                     1.555) #Be - 2nd ceiling coefficient

gen3_avr_benchmark = AVRTypeII(20.0, #K0 - Gain
                  0.2, #T1 - 1st pole
                  0.063, #T2 - 1st zero
                  0.35, #T3 - 2nd pole
                  0.01, #T4 - 2nd zero
                  0.314, #Te - Field current time constant
                  0.001, #Tr - Measurement time constant
                  5.0, #Vrmax
                  -5.0, #Vrmin
                  0.0039, #Ae - 1st ceiling coefficient
                  1.555) #Be - 2nd ceiling coefficient

gen2_avr_benchmark_m = AVRTypeIIManual(20.0, #Ka - Gain
                   0.2, #Ta - 1st pole
                   0.063, #Kf - 1st zero
                   0.35, #Tf - 2nd pole
                   0.01, #Te - 2nd zero
                   0.314, #Tr - Field current time constant
                   5.0, #Vrmax
                   -5.0, #Vrmin
                   0.0039, #Ae - 1st ceiling coefficient
                   1.555) #Be - 2nd ceiling coefficient
gen3_avr_benchmark_m = AVRTypeI(20.0, #Ka - Gain
                                0.01, #Ke
                                0.063, #Kf
                                0.2, #Ta
                                0.314, #Te
                                0.35, #Tf
                                0.001, #Tr
                                5.0, #Vrmax
                                -5.0, #Vrmin
                                0.0039, #Ae - 1st ceiling coefficient
                                1.555) #Be - 2nd ceiling coefficient

gen2_avr_fixed = AVRFixed(2.1996)
gen3_avr_fixed = AVRFixed(1.7327)
gen3_avr_benchmark_m = AVRTypeIIManual(20.0, #Ka - Gain
                0.2, #Ta - 1st pole
                0.063, #Kf - 1st zero
                0.35, #Tf - 2nd pole
                0.01, #Te - 2nd zero
                0.314, #Tr - Field current time constant
                5.0, #Vrmax
                -5.0, #Vrmin
                0.0039, #Ae - 1st ceiling coefficient
                1.555) #Be - 2nd ceiling coefficient

######################### Generators ########################

Gen1AVR = DynGenerator(1, #Number
                 :TestGen,
                 nodes_OMIB[2],#bus
                 1.0, # ω_ref,
                 1.05,
                 0.4,
                 Basic,
                 BaseShaft,
                 proportional_avr, #avr
                 fixed_tg, #tg
                 no_pss)

Gen1AVRnoAVR = DynGenerator(1, #Number
                 :TestGen,
                 nodes_OMIB[2],#bus
                 1.0, # ω_ref,
                 1.05,
                 0.4,
                 Basic,
                 BaseShaft,
                 fixed_avr, #avr
                 fixed_tg, #tg
                 no_pss)

Gen2AVRnoAVR = DynGenerator(1, #Number
                 :TestGen,
                 nodes_OMIB[2],#bus
                 1.0, # ω_ref,
                 1.02,
                 0.4,
                 oneDoneQ,
                 BaseShaft,
                 fixed_avr, #avr
                 fixed_tg, #tg
                 no_pss)

Gen2AVR = DynGenerator(1, #Number
                 :TestGen,
                 nodes_OMIB[2],#bus
                 1.0, # ω_ref,
                 1.02,
                 0.4,
                 oneDoneQ,
                 BaseShaft,
                 proportional_avr, #avr
                 fixed_tg, #tg
                 no_pss)

Gen3AVR = DynGenerator(1, #number
                       :TestGen, #name
                       nodes_OMIB[2], #bus
                       1.0, # ω_ref
                       1.02,
                       0.4,
                       oneDoneQ, #machine
                       BaseShaft,
                       typeI_avr, #avr
                       fixed_tg, #tg
                       no_pss)

Gen4AVRTG = DynGenerator(1, #number
                       :TestGen, #name
                       nodes_OMIB[2], #bus
                       1.0, # ω_ref
                       1.02,
                       0.4,
                       oneDoneQ, #machine
                       BaseShaft,
                       typeI_avr, #avr
                       typeI_tg, #tg
                       no_pss)

Gen5AVR = DynGenerator(1, #number
                       :TestGen, #name
                       nodes_OMIB[2], #bus
                       1.0, # ω_ref
                       1.02,
                       0.4,
                       AndersonFouad, #machine
                       BaseShaft,
                       typeI_avr, #avr
                       fixed_tg, #tg
                       no_pss)

Gen6AVR = DynGenerator(1, #number
                       :TestGen, #name
                       nodes_OMIB[2], #bus
                       1.0, #ω_ref
                       1.02, #V_ref
                       0.4, #P_ref
                       oneDoneQ, #machine
                       FiveShaft, #shaft
                       typeI_avr, #avr
                       fixed_tg, #tg
                       no_pss)

GenKundur = DynGenerator(1, #number
                         :KundurGen, #name
                         nodes_kundur[1],
                         1.0, #ω_ref
                         1.02, #V_ref
                         0.2, #P_ref
                         KundurMachine, #machine
                         KundurShaft, #shaft
                         kundur_avr, #avr
                         fixed_tg, #tg
                         no_pss) #pss

GenFullKundur =  DynGenerator(1, #number
                              :KundurGen, #name
                              nodes_kundur[1],
                              1.0, #ω_ref
                              1.02, #V_ref
                              0.20, #P_ref
                              KundurFullMachine, #machine
                              KundurShaft, #shaft
                              kundur_avr, #avr
                              fixed_tg, #tg
                              no_pss) #pss

Gen2_benchmark = DynGenerator(1, #number
                              :Gen1, #name
                              nodes_benchmark[2],
                              1.0, #ω_ref
                              1.0142, #V_ref
                              1.0, #P_ref
                              Mach2_benchmark, #machine
                              Shaft2_benchmark, #shaft
                              gen2_avr_benchmark_m, #avr
                              fixed_tg, #tg
                              no_pss) #pss
Gen2_benchmark_tg = DynGenerator(1, #number
                              :Gen1, #name
                              nodes_benchmark[2],
                              1.0, #ω_ref
                              1.0142, #V_ref
                              1.0, #P_ref
                              Mach2_benchmark, #machine
                              Shaft2_benchmark, #shaft
                              gen2_avr_benchmark_m, #avr
                              typeII_tg, #tg
                              no_pss) #pss

Gen3_benchmark = DynGenerator(2, #number
                                :Gen2, #name
                                nodes_benchmark[3],
                                1.0, #ω_re
                                1.0059, #V_ref
                                1.0, #P_ref
                                Mach3_benchmark, #machine
                                Shaft3_benchmark, #shaft
                                gen3_avr_benchmark_m, #avr
                                fixed_tg, #tg
                                no_pss) #pss

Gen2_benchmark_fixedavr = DynGenerator(1, #number
                              :Gen1, #name
                              nodes_benchmark[2],
                              1.0, #ω_ref
                              1.0142, #V_ref
                              1.0, #P_ref
                              Mach2_benchmark, #machine
                              Shaft2_benchmark, #shaft
                              gen2_avr_fixed, #avr
                              fixed_tg, #tg
                              no_pss) #pss

Gen3_benchmark_fixedavr = DynGenerator(2, #number
                                :Gen2, #name
                                nodes_benchmark[3],
                                1.0, #ω_ref
                                1.00586, #V_ref
                                1.0, #P_ref
                                Mach3_benchmark, #machine
                                Shaft3_benchmark, #shaf
                                gen3_avr_fixed, #avr
                                fixed_tg, #tg
                                no_pss) #pss

Gen_bench = DynGenerator(1, #number
                              :Gen1, #name
                              nodes_bench[2],
                              1.0, #ω_ref
                              1.0, #V_ref
                              1.0, #P_ref
                              Mach2_benchmark, #machine
                              Shaft2_benchmark, #shaft
                              gen_avr_bench, #avr
                              fixed_tg, #tg
                              no_pss) #pss
                                gen3_avr_benchmark_m, #avr
                                fixed_tg, #tg
                                no_pss) #pss

################### Sources #####################

inf_gen = StaticSource(1, #number
                 :InfBus, #name
                 nodes_OMIB[1],#bus
                 1.05, #VR
                 0.0, #VI
                 0.000005) #Xth

inf_gen_kundur = StaticSource(1,
                        :InfBus,
                        nodes_kundur[2],
                        1.01,
                        0.0,
                        0.000005)

inf_gen_benchmark = StaticSource(1, #number
                 :InfBus, #name
                 nodes_benchmark[1],#bus
                 1.02, #VR
                 0.0, #VI
                 0.000001) #Xth

################### Sources #####################
inf_gen_kundur = Source(1,
                        :InfBus,
                        nodes_kundur[2],
                        1.01,
                        0.0,
                        0.000005)

inf_gen_benchmark = Source(1, #number
                 :InfBus, #name
                 nodes_benchmark[1],#bus
                 1.02, #VR
                 0.0, #VI
                 0.000001) #Xth

################### Sources #####################
OMIB =  DynamicSystem(nodes_OMIB, branch_OMIB, [Gen1AVR], [inf_gen], 100.0, 60.0)
OMIBA = DynamicSystem(nodes_OMIB, branch_OMIB, [Gen5AVR], [inf_gen], 100.0, 60.0)
OMIB_Kundur = DynamicSystem(nodes_kundur, branch_kundur, [GenKundur], [inf_gen_kundur], 100.0, 60.0)
OMIB_FullKundur = DynamicSystem(nodes_kundur, branch_kundur, [GenFullKundur], [inf_gen_kundur], 100.0, 60.0
TwoBus_Benchmark = DynamicSystem(nodes_bench, #buses
                                   branch_bench, #branches
                                   [Gen_bench], #dynamic injections
                                   [inf_gen_benchmark, loads_bench[1]], #static injections
                                   100.0, #Sbase
                                   60.0) #fbase

ThreeBus_Benchmark = DynamicSystem(nodes_benchmark, #buses
                                  branch_benchmark, #branches
                                  [Gen2_benchmark, Gen3_benchmark], #dynamic injections
                                  [inf_gen_benchmark, loads_benchmark[1], loads_benchmark[2], loads_benchmark[3]], #static injections
                                  100.0, #Sbase
                                  60.0) #fbase

ThreeBus_Benchmark = DynamicSystem(nodes_benchmark, #buses
                                  branch_benchmark, #branches
                                  [Gen2_benchmark, Gen3_benchmark], #dynamic injections
                                  [inf_gen_benchmark, loads_benchmark[1], loads_benchmark[2], loads_benchmark[3]], #static injections
                                  100.0, #Sbase
                                  60.0) #fbase

ThreeBus_Benchmark_tg = DynamicSystem(nodes_benchmark, #buses
                                  branch_benchmark, #branches
                                  [Gen2_benchmark_tg, Gen3_benchmark], #dynamic injections
                                  [inf_gen_benchmark, loads_benchmark[1], loads_benchmark[2], loads_benchmark[3]], #static injections
                                  100.0, #Sbase
                                  60.0) #fbase

ThreeBus_Benchmark_fixedavr = DynamicSystem(nodes_benchmark, #buses
                                branch_benchmark, #branches
                                [Gen2_benchmark_fixedavr, Gen3_benchmark_fixedavr], #dynamic injections
                                [inf_gen_benchmark,loads_benchmark[1], loads_benchmark[2], loads_benchmark[3]], #static injections
                                100.0, #Sbase
                                60.0) #fbase
ThreeBus_Benchmark = DynamicSystem(nodes_benchmark, #buses
                                   branch_benchmark, #branches
                                   [Gen2_benchmark, Gen3_benchmark], #dynamic injections
                                   [inf_gen_benchmark, loads_benchmark[1], loads_benchmark[2], loads_benchmark[3]], #static injections
                                   100.0, #Sbase
                                   60.0) #fbase


text:
abstract type YearPeriod end

struct YearQuarter <: YearPeriod
  y::UInt16
  q::UInt8

  function YearQuarter(y,q)
    q == 1 || q==2 || q==3 || q==4 || error("Invalid YearQuarter $y $q")
    return new(y,q)
  end
end

#creates a year quarter from the yyyy.q format or from yyyy (assuming eoy)
#unpredicatable results if the decimal is not 0 (year only case), 1,2,3, or 4
function YearQuarter(f::Real)::YearQuarter
  ix10::Int = Int(f*10)
  y::Int = ix10 ÷ 10
  q::Int = ix10 - y*10
  (q==0) && (q=4) #if a year is passed, assume it is the end of the year (q4)

  return YearQuarter(y,q)
end

function YearQuarter(d::Date)
  y::UInt16 = year(d)
  q::UInt8 = quarterofyear(d)

  return YearQuarter(y,q)
end

#fallback in case we need these
Float64(yq::YearQuarter) = yq.y + yq.q/10
Base.show(io::IO, yq::YearQuarter) = print(io, Float64(yq))

#quick way to split the year quarter
Tuple{Int,Int}(yq::YearQuarter)::Tuple{Int,Int} = (yq.y, yq.q)

#allows for adding quarters
function +(yq::YearQuarter, qadded::Int)::YearQuarter
  (y₀::Int, q₀::Int) = Tuple{Int,Int}(yq)

  q::Int = (q₀ + qadded) % 4
  (q≤0) && (q+=4) #only allow 1,2,3,4

  y::Int = (qadded - (q-q₀)) ÷ 4 + y₀

  @assert q-q₀ + 4*(y-y₀) == qadded #this should never fail

  return YearQuarter(y,q)
end

#helper methods to build on the above
+(qadded::Int, yq::YearQuarter)::YearQuarter = yq + qadded
-(yq::YearQuarter, qadded::Int)::YearQuarter = yq + (-1*qadded)

#gets the number of quarters difference
function -(yq1::YearQuarter, yq2::YearQuarter)::Int

  Δy::Int = yq1.y - yq2.y
  Δq::Int = yq1.q - yq2.q

  return Δy * 4 + Δq
end

#comparison operators (primary sorting functions, then operator symbols)
isless(yq1::YearQuarter, yq2::YearQuarter) = yq1.y*10+yq1.q < yq2.y*10 + yq2.q

boq(yq::YearQuarter)::Date = Date(yq.y, (yq.q-1)*3,1)
eoq(yq::YearQuarter)::Date = lastdayofquarter(boq(yq))

##################################Year-Month#####################
struct YearMonth <: YearPeriod
  y::UInt16
  m::UInt8

  function YearMonth(y,m)
    (m≥1 && m≤12) || error("Invalid YearMonth $y $m")
    return new(y,m)
  end
end

#creates a year month from yyyymm date style
function YearMonth(i::Integer)::YearMonth
  local y::UInt16
  local m::UInt8

  (i>1000_00 && i<10000_00) || error("Invalid YearMonth $i")
  (y,m) = (i ÷ 100, i % 100)

  return YearMonth(y,m)
end

#fallback in case we need these
Int(ym::YearMonth)::Int = ym.y*100 + ym.m
Base.show(io::IO, ym::YearMonth) = print(io, "$(ym.y)_$(ym.m)")

#quick way to split the year quarter
Tuple{Int,Int}(ym::YearMonth)::Tuple{Int,Int} = (ym.y, ym.m)

#allows for adding quarters
function +(ym::YearMonth, madded::Integer)::YearMonth
  (y₀::Int, m₀::Int) = Tuple{Int,Int}(ym)

  m::Int = (m₀ + madded) % 12
  (m≤0) && (m+=12) #only allow positive months

  y::Int = (madded - (m-m₀)) ÷ 12 + y₀

  @assert m-m₀ + 12*(y-y₀) == madded #this should never fail

  return YearMonth(y,m)
end

#helper methods to build on the above
+(qadded::Int, yq::YearMonth)::YearMonth = ym + qadded
-(ym::YearMonth, madded::Int)::YearMonth = ym + (-1*madded)

#gets the number of months difference
function -(ym1::YearMonth, ym2::YearMonth)::Int

  Δy::Int = ym1.y - ym2.y
  Δm::Int = ym1.m - ym2.m

  return Δy * 12 + Δm
end

#comparison operators (primary sorting functions, then operator symbols)
isless(ym1::YearMonth, ym2::YearMonth) = ym1.y*100+ym1.m < ym2.y*100 + ym2.m

#convert back to a date format
bom(ym::YearMonth)::Date = Date(ym.y, ym.m, 1)
eom(ym::YearMonth)::Date = lastdayofmonth(bom(ym))

Base.length(::YearPeriod) = 1
Base.iterate(yp::YearPeriod) = yp
Base.iterate(ym::YearPeriod, s::Integer) = nothing
Base.Broadcast.broadcastable(yp::YearPeriod) = Ref(yp)

ext: jl

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Flux
using PyCallChainRules.Torch: TorchModuleWrapper, torch

###
# WARNING: Doesn't work correctly yet.
###

input_dim = 4
output_dim = 2
hiddendim = 16

batchsize = 8

torch_module = torch.nn.Sequential(
                            torch.nn.Linear(input_dim, hiddendim), torch.nn.ReLU(),
                            torch.nn.Linear(hiddendim, hiddendim), torch.nn.ReLU(),
                            torch.nn.Linear(hiddendim, output_dim)
                        )
jlmodel = TorchModuleWrapper(torch_module)

input = randn(Float32, input_dim, batchsize)
target = randn(Float32, output_dim, batchsize)

loss(x, y) = Flux.Losses.mse(jlmodel(x), y)
@info "before" loss(input, target)

function train(model; nsteps=100)
    opt = Flux.ADAM(0.01)
    ps = Flux.Zygote.Params(Flux.params(model))
    for i in 1:nsteps
        gs, = Flux.Zygote.gradient(ps) do
            loss(input, target)
        end
        Flux.Optimise.update!(opt, ps, gs)
    end
    return model
end

jlmodel = train(jlmodel)

@info "after" loss(input, target)


text:
using ModelSanitizer
using Test

function _inner_runtests()::Nothing
    @testset "Unit tests" begin
        @debug("Running unit tests...")
        @testset "unit-tests/test-dataframes.jl" begin
            @debug("unit-tests/test-dataframes.jl")
            include("unit-tests/test-dataframes.jl")
        end
        @testset "unit-tests/test-elements.jl" begin
            @debug("unit-tests/test-elements.jl")
            include("unit-tests/test-elements.jl")
        end
        @testset "unit-tests/test-sanitize.jl" begin
            @debug("unit-tests/test-sanitize.jl")
            include("unit-tests/test-sanitize.jl")
        end
        @testset "unit-tests/test-utils.jl" begin
            @debug("unit-tests/test-utils.jl")
            include("unit-tests/test-utils.jl")
        end
        @testset "unit-tests/test-zero.jl" begin
            @debug("unit-tests/test-zero.jl")
            include("unit-tests/test-zero.jl")
        end
    end
    @testset "Integration tests" begin
        @debug("Running integration tests...")
        @testset "integration-tests/test-proof-of-concept-dataframes.jl" begin
            @debug("integration-tests/test-proof-of-concept-dataframes.jl")
            include("integration-tests/test-proof-of-concept-dataframes.jl")
        end
        @testset "integration-tests/test-proof-of-concept-linearmodel" begin
            @debug("integration-tests/test-proof-of-concept-linearmodel.jl")
            include("integration-tests/test-proof-of-concept-linearmodel.jl")
        end
        if Base.JLOptions().can_inline > 0
            @testset "integration-tests/test-proof-of-concept-mlj.jl" begin
                @debug("integration-tests/test-proof-of-concept-mlj.jl")
                include("integration-tests/test-proof-of-concept-mlj.jl")
            end
        end
    end
    return nothing
end

function runtests()::Nothing
    @sync begin
        _test_runner = @async begin
            _inner_runtests()
        end
        @async begin
            while !istaskdone(_test_runner)
                sleep(60)
                @debug("[[ModelSanitizer tests are still running...]]")
            end
        end
    end
    return nothing
end

@testset "ModelSanitizer.jl" begin
    runtests()
end


text:
# tests
using NearestNeighborDescent
using NearestNeighborDescent: NNTuple, _neighbors, brute_knn, _heappush!
using DataStructures
using Distances
using Test

include("utils_tests.jl")
include("nndescent_tests.jl")


text:
module test_julia_number

using Test

@test iszero(zero(1))
@test isone(one(0))

end # module test_julia_number

module test_julia_number_adjoint

using Test

@test im' == -im
@test im'' == im
@test conj(im) == (transpose ∘ conj)(im) == adjoint(im)

end # module test_julia_number_adjoint

module test_julia_number_isless

using Test

@test isless(-0.0, 0.0)
@test !(-0.0 < 0.0)

end # module test_julia_number_isless


text:
@kwdef struct BalancedAcceptance{I <: Integer} <: SpatialSampler
    numpoints::I = 50
    α = 1.0
end

function _generate!(piv::BalancedAcceptance, uncertainty::M) where {M<:AbstractMatrix}
    seed = Int32.([floor(10^7*rand()),floor(10^7*rand())])
    np, α = piv.numpoints, piv.α
    x,y = size(uncertainty)

    stduncert = StatsBase.transform(StatsBase.fit(ZScoreTransform, uncertainty, dims=2), uncertainty)
    reluncert = broadcast(x->exp(α*x)/(1+exp(α*x)), stduncert)
    coords = []

    ptct = 0
    while length(coords) < np
        i,j = haltonvalue(seed[1]+ptct,2), haltonvalue(seed[2]+ptct,3)
        candcoord = CartesianIndex(convert.(Int32, [ceil(x*i), ceil(y*j)])...)
        prob = reluncert[candcoord]
        rand() < prob && push!(coords, candcoord)
        ptct += 1
    end

    return coords
end


text:
## activate project environment
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

## compute relative performance
using DelimitedFiles
using Measurements

function compute_relative_performance(filename)
  header = "# Polydeg Primitive/conservative-variables-mean std"
  data = readdlm(filename, comments=true)
  cons = data[:, 2] .± data[:, 3]
  prim = data[:, 4] .± data[:, 5]
  ratio = prim ./ cons
  open(filename[1:end-4] * "_relative.dat", "w") do io
    println(io, header)
    writedlm(io, hcat(data[:, 1], # polydeg
                      Measurements.value.(ratio),
                      Measurements.uncertainty.(ratio)))
  end
end

compute_relative_performance(joinpath(@__DIR__, "primitive_2D_flux_shima_etal.dat"))
compute_relative_performance(joinpath(@__DIR__, "primitive_3D_flux_shima_etal.dat"))
compute_relative_performance(joinpath(@__DIR__, "primitive_2D_flux_ranocha.dat"))
compute_relative_performance(joinpath(@__DIR__, "primitive_3D_flux_ranocha.dat"))


text:
mutable struct Rotational{T,Nc} <: Joint{T,Nc}
    V12::SMatrix{2,3,T,6}
    V3::Adjoint{T,SVector{3,T}}
    qoff::Quaternion{T}
    cid::Int64

    function Rotational{T,Nc}(body1::AbstractBody{T}, body2::AbstractBody{T}; axis::AbstractVector{T} = zeros(3), offset::Quaternion{T} = Quaternion{T}()) where {T,Nc}
        axis = vrotate(SVector(axis...), inv(offset))
        if norm(axis) != 0
            axis = axis / norm(axis)
        end
        A = svd(skew(axis)).Vt # in frame of body1
        V12 = A[1:2,:]
        V3 = axis' # instead of A[3,:] for correct sign: abs(axis) = abs(A[3,:])
        cid = body2.id

        new{T,Nc}(V12, V3, offset, cid), body1.id, body2.id
    end
end

Rotational0 = Rotational{T,0} where T
Rotational1 = Rotational{T,1} where T
Rotational2 = Rotational{T,2} where T
Rotational3 = Rotational{T,3} where T


text:
using NeuralQuantum
using Test
using NeuralQuantum: LdagLSparseSuperopProblem, LRhoSparseOpProblem, LRhoKLocalOpProblem

Nsites = 4
T = Float64
lind = quantum_ising_lind(SquareLattice([Nsites],PBC=true), g=1.3, V=2.0, γ=1.0)

function test_ldagl_op(T, Nsites, net, lind)

    prob  = LRhoSparseOpProblem(T, lind);
    probL = LRhoKLocalOpProblem(T, lind);

    v    = state(prob, net)
    vL   = state(probL, net)

    L=liouvillian(lind)
    LdagL = L.data'*L.data
    rho   = dm(net, prob, false).data
    rhov  = vec(rho)
    Clocs_ex  = abs.((L.data*rhov)./rhov).^2

    ic     = NeuralQuantum.MCMCSRLEvaluationCache(net, prob);  zero!(ic)
    icL_fb = NeuralQuantum.MCMCSRLEvaluationCache(net, prob);  zero!(icL_fb)
    icL    = NeuralQuantum.MCMCSRLEvaluationCache(net, probL); zero!(icL)

    for i=1:spacedimension(v)
        set_index!(v, i)
        set_index!(vL, i)
        NeuralQuantum.sample_network!(ic,     prob,  net, v);
        NeuralQuantum.sample_network!(icL_fb, probL, net, v);
        NeuralQuantum.sample_network!(icL,    probL, net, vL);
    end

    SREval     = EvaluatedNetwork(SR(), net)
    SREvalL    = EvaluatedNetwork(SR(), net)
    SREvalL_fb = EvaluatedNetwork(SR(), net)

    evaluation_post_sampling!(SREval,     ic)
    evaluation_post_sampling!(SREvalL,    icL)
    evaluation_post_sampling!(SREvalL_fb, icL_fb)

    @test Clocs_ex ≈ ic.Evalues

    @testset "LookUp table evaluation" begin
        @test Clocs_ex ≈ icL.Evalues

        @test SREval.L ≈ SREvalL.L
        @test ic.Evalues ≈ icL.Evalues
        @test all([l≈r for (l,r)=zip(SREval.F, SREvalL.F)])
        @test all([l≈r for (l,r)=zip(SREval.S, SREvalL.S)])
    end
    @testset "Fallback (no LUT)" begin
        @test Clocs_ex ≈ icL_fb.Evalues

        @test SREval.L ≈ SREvalL_fb.L
        @test ic.Evalues ≈ icL_fb.Evalues
        @test all([l≈r for (l,r)=zip(SREval.F, SREvalL_fb.F)])
        @test all([l≈r for (l,r)=zip(SREval.S, SREvalL_fb.S)])
    end
end

@testset "LdagL (operator) problem" begin
    @testset "RBMSplit $netT" for netT=[T, Complex{T}]
        net  = cached(RBMSplit(netT, Nsites, 2))
        test_ldagl_op(T, Nsites, net, lind)
    end
    @testset "NDM $T" for T=[T]
        net  = cached(NDM(T, Nsites, 2, 1))
        test_ldagl_op(T, Nsites, net, lind)
    end
end


text:
module TestTasks

# using Revise
using Test
using MLJBase

XX = (Crim = [0.00632, 0.02731, 0.02729],
      Zn = [18.0, 0.0, 0.0],
      Indus = [2.31, 7.07, 7.07],
      Chas = [0, 0, 0],
      NOx = [0.538, 0.469, 0.469],
      Rm = [6.575, 6.421, 7.185],
      Age = [65.2, 78.9, 61.1],
      Dis = [4.09, 4.9671, 4.9671],
      Rad = [1.0, 2.0, 2.0],
      Tax = [296.0, 242.0, 242.0],
      PTRatio = [15.3, 17.8, 17.8],
      Black = [396.9, 396.9, 392.83],
      LStat = [4.98, 9.14, 4.03],
      MedV = [24.0, 21.6, 34.7],)

allnames = collect(MLJBase.schema(XX).names)

task = SupervisedTask(data=XX, target=[:Crim, :Zn], is_probabilistic=true, ignore=:Dis)

y = y_(task);
X = X_(task);
t = MLJBase.schema(X);
@test collect(t.names) == filter(allnames) do ftr
    !(ftr in [:Crim, :Zn, :Dis])
end
y1 = [(XX.Crim[i], XX.Zn[i]) for i in 1:length(y)];
@test y == y1

task = SupervisedTask(data=XX, target=:Crim, is_probabilistic=true, ignore=[:Dis, :Rm])
y = y_(task);
X = X_(task);
t = MLJBase.schema(X)
@test collect(t.names) == filter(allnames) do ftr
    !(ftr in [:Crim, :Dis, :Rm])
end

task = UnsupervisedTask(data=XX, ignore=[:Dis, :Rm])
X = task();
t = MLJBase.schema(X)
@test collect(t.names) == filter(allnames) do ftr
    !(ftr in [:Dis, :Rm])
end

task = UnsupervisedTask(data=XX, ignore=:Rm)
X = X_(task);
t = MLJBase.schema(X)
@test collect(t.names) == filter(allnames) do ftr
    !(ftr in [:Rm, ])
end

# single feature for input:
task = UnsupervisedTask(data=(Crim=XX.Crim,))
@test task.X == XX.Crim
task = SupervisedTask(data=MLJBase.selectcols(XX, [:Crim, :Zn]), target=:Zn, is_probabilistic=true)
@test task.X == XX.Crim
@test task.y == XX.Zn

end # module
true


text:
struct BipartiteMatching
    weight :: Int64
    match  :: Vector{Int64}
end

function bipartite_cardinality_matching(l_in::Vector{Int}, r_in::Vector{Int}, m, n; l_sorted=false)
    @assert length(l_in) == length(r_in)
    @assert m <= n
    l = l_in
    r = r_in
    if !l_sorted
        perm = sortperm(l_in)
        l = l_in[perm]
        r = r_in[perm]
    end

    len = length(l)
    matching_l = zeros(Int, m)
    matching_r = zeros(Int, n)

    # create initial matching
    match_len = 0
    for (li,ri) in zip(l,r)
        if matching_l[li] == 0 && matching_r[ri] == 0
            matching_l[li] = ri
            matching_r[ri] = li
            match_len += 1
        end
    end

    if match_len < m
        # creating indices be able to get edges a vertex is connected to
        # only works if l is sorted
        index_l = zeros(Int, m+1)
        last = l[1]
        c = 2
        for i = 2:len
            if l[i] != last
                index_l[last+1] = c
                last = l[i]
            end
            c += 1
        end
        index_l[l[end]+1] = c

        index_l[1] = 1

        process_nodes = zeros(Int, m+n)
        depths = zeros(Int, m+n)
        parents = zeros(Int, m+n)
        used_l = zeros(Bool, m)
        used_r = zeros(Bool, n)
        found = false

        # find augmenting path
        while match_len < m
            pend = 1
            pstart = 1
            for li in l
                # free vertex
                if matching_l[li] == 0
                    process_nodes[pstart] = li
                    depths[pstart] = 1
                    break
                end
            end

            begin
            while pstart <= pend
                node = process_nodes[pstart]
                depth = depths[pstart]

                # from left to right
                if depth % 2 == 1
                    used_l[node] = true
                    # only works if l is sorted
                    for ri=index_l[node]:index_l[node+1]-1
                        child_node = r[ri]
                        # don't use matching edge
                        if matching_r[child_node] != node && !used_r[child_node]
                            used_r[child_node] = true
                            pend += 1
                            depths[pend] = depth+1
                            process_nodes[pend] = child_node
                            parents[pend] = pstart
                        end
                    end
                else # right to left (only matching edge)
                    # if matching edge
                    match_to = matching_r[node]
                    if match_to != 0
                        if !used_l[match_to]
                            used_l[match_to] = true
                            pend += 1
                            depths[pend] = depth+1
                            process_nodes[pend] = match_to
                            parents[pend] = pstart
                        end
                    else
                        # found augmenting path
                        parent = pstart
                        last = 0
                        c = 0
                        while parent != 0
                            current = process_nodes[parent]
                            if last != 0
                                if c % 2 == 1
                                    matching_r[last] = current
                                    matching_l[current] = last
                                end
                            end
                            c += 1
                            last = current
                            parent = parents[parent]
                        end
                        # break because we found a path
                        found = true
                        break
                    end
                end
                pstart += 1
            end
            if found
                match_len += 1
                if match_len < m
                    used_l .= false
                    used_r .= false
                end
                found = false
            else
                break
            end
            end
        end
    end
    return BipartiteMatching(match_len, matching_l)
end


text:
using LinearAlgebraicRepresentation
using Plasm, LinearAlgebra
Lar = LinearAlgebraicRepresentation
using SparseArrays

include("./LARinput.jl")

# s = Lar.Struct(store)
# V,CV,FV,EV = Lar.struct2lar(s)
# V = Plasm.normalize3D(V) TODO:  solve MethodError bug
# Plasm.view(V,CV)

cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_EW = convert(Lar.ChainOp, cop_EV);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(
	W::Lar.Points, cop_EW::Lar.ChainOp, cop_FE::Lar.ChainOp)

cc = [copEV, copFE, copCF]
output = Lar.lar2obj(V::Lar.Points, cc::Lar.ChainComplex)

f = open("test/out3d.obj", "w")
print(f, output)
close(f)

V,EVs,FVs = Lar.obj2lar("test/out3d.obj")

# for k=1:length(FVs)
# 	Plasm.view(V, FVs[k])
# end

Plasm.view(V, FVs[1])


text:
struct CenterCropCat
    trms
end

function CenterCropCat(t, vol)
    t = vol ? (t, t, t) : (t, t)
    CenterCropCat((t..., 0, 0))
end

function (cc::CenterCropCat)(x1, x2) # x2 : input
    lo = cc.trms .+ 1
    up = size(x2) .- cc.trms
    chcat(x2[UnitRange.(lo, up)...], x1)
end

"""
    unet(; inchannels, nclasses = 1, volume = false, base = 64,
    batchnorm = false, padding = false, upsample = :convt, nlevels = 4)

Build a [U-Net](https://arxiv.org/abs/1505.04597v1) to process images with
`inchannels` channels (e.g. `inchannels` = 3 for RGB images). Default argument
values correspond to the original paper implementation with unpadded
convolutions.

* `nclasses`: number of pixel classes. If `nclasses` = 2 the final layer have
  one channel with a sigmoid activation function, otherwise the final layer have
  `nclasses` channels followed by a softmax function.
* `volume`: set to `true` to process tri-dimensional datas
* `basewidth`: base number of channels, the number of filters/channels is
  multiplied by two at each U-Net level. Notice that `basewidth` should be a
  power of two number, other values could lead to unfunctionnal U-Net.
* `batchnorm`: if `true` add a `BatchNorm` layer after each convolution
* `padding`: if `true` convolutions are padded
* `upsample`: either `:convt` (`ConvTranspose` layer), `:nearest` or `:bilinear`
  (`Upsample layer`)
* `nlevels`: number of level or depth of the U-Net
"""
function unet(; inchannels,
                nclasses = 2,
                volume = false,
                basewidth = 64,
                # expansion = 2,
                batchnorm = false,
                padding = false,
                upsample = :convt,
                nlevels = 4)
    if volume && upsample == :bilinear
        throw(ArgumentError(
            "Bilinear upsampling can not be used with tri-dimensionnal datas."))
    end

    # Pack parameters in a named tuple for convenience
    pars = (k₁ = volume ? (1, 1, 1) : (1, 1),
            k₂ = volume ? (2, 2, 2) : (2, 2),
            k₃ = volume ? (3, 3, 3) : (3, 3),
            nclasses = nclasses,
            inchannels = inchannels,
            basewidth = basewidth,
            batchnorm = batchnorm,
            padding = padding,
            upsample = upsample)

    # Build encoders and decoders for each level
    enc, dec, con = [], [], []
    for l ∈ 1:nlevels
        # e, d = level_enc_dec(l, pars)
        # push!(enc, e)
        # push!(dec, d)
        push!.((enc, dec), level_enc_dec(l, pars))
        if padding
            push!(con, chcat)
        else
            push!(con, CenterCropCat(trim(l, nlevels), volume))
        end
    end

    # Build bridge (i.e. bottom part)
    ec, dc = level_channels(nlevels + 1, pars)
    ic, mc, _ = ec
    _, _, oc = dc
    dw, up = samplers(oc, pars)
    bdg  = [dw, double_conv(ic, mc, oc, pars)..., up]

    uchain(encoders = enc, decoders = dec, bridge = bdg, connection = con)
end

# Build encoder and decoder blocks for a given level
function level_enc_dec(lvl, pars)
    # number of channels (input, middle, output)
    ec, dc = level_channels(lvl, pars)
    enc = []
    dec = []
    if lvl == 1
        enc = double_conv(ec, pars)
        # build classifier layers
        k₁ = pars[:k₁]
        nc = pars[:nclasses]
        if nc > 2
            d = length(k₁)
            cl = [Conv(k₁, dc[3] => nc), x -> softmax(x, dims = d)]
        else
            cl = [Conv(k₁, dc[3] => 1, sigmoid)]
        end
        dec = [double_conv(dc, pars)..., cl...]
    else
        dw, up = samplers(dc[3], pars)
        enc = [dw, double_conv(ec, pars)...]
        dec = [double_conv(dc, pars)..., up]
    end
    enc, dec
end

# Return two convolutional layers, with batchnom if needed
function double_conv(chs, pars)
    ic, mc, oc = chs
    padding, k₃, batchnorm = pars[:padding], pars[:k₃], pars[:batchnorm]
    p = padding ? 1 : 0
    cv1(chs) = [Conv(k₃, chs, pad = p, bias = false), BatchNorm(chs[2], relu)]
    cv2(chs) = [Conv(k₃, chs, relu, pad = p)]
    cnv(chs) =  batchnorm ? cv1(chs) : cv2(chs)
    [cnv(ic => mc)..., cnv(mc => oc)...]
end

double_conv(ic, mc, oc, pars) = double_conv((ic, mc, oc), pars)

# Compute encoder/decoder number of channels at a given level (lvl)
# return two tuples one for encoder and one for decoder
# formula :
#  ice = expansion^(lvl - 2) * basewidth
#  mce = expansion^(lvl - 1) * basewidth
function level_channels(lvl, pars)
    inchannels, basewidth = pars[:inchannels], pars[:basewidth]
    upsample = pars[:upsample]
    # encoder channels: input, middle, ouptput = (ice, mce, oce)
    ice = (lvl == 1) ? inchannels : 2^(lvl - 2) * basewidth
    mce = oce = 2^(lvl - 1) * basewidth
    # decoder channels: input, middle, ouptput = (icd, mcd, ocd)
    icd, mcd = 2 * oce, oce
    if upsample == :convt
        ocd = mcd
    else
        ocd = (lvl == 1) ? mcd : mcd ÷ 2
    end
    (ice, mce, oce), (icd, mcd, ocd)
end

# Downsampler, Upsampler
function samplers(chs, pars)
    k₂, upsample = pars[:k₂], pars[:upsample]
    dw = MaxPool(k₂)
    if upsample == :convt
        up =  ConvTranspose(k₂, chs => (chs ÷ 2), stride = 2)
    else
        up = Upsample(upsample, scale = k₂)
    end
    dw, up
end


text:
include("impl.jl")

for (front_name, backend) in (
        :conv          => :_nnpack,
        :∇conv_data    => :_nnpack,
        :∇conv_filter  => :_nnpack,
    )
    @eval begin
        function $(Symbol("$(front_name)$(backend)!"))(
                        out::Array{T1,4}, in1::Array{T2,4}, in2::Array{T3,4},
                        cdims::ConvDims; kwargs...) where {T1, T2, T3}
            @warn "Automatically converting input tensor to Float32. This will have performance implications" maxlog=1
            # Output must of the same type as in the function signature
            T1.($(Symbol("$(front_name)$(backend)!"))(Float32.(out), Float32.(in1),
                                                      Float32.(in2), cdims; kwargs...))
        end
    end
end

function conv_nnpack(x::Array{T1, 4}, w::Array{T2, 4}, cdims::ConvDims; kwargs...) where {T1, T2}
    y = similar(x, output_size(cdims)..., channels_out(cdims), size(x, 4))
    return conv_nnpack!(y, x, w, cdims; kwargs...)
end

function ∇conv_data(dy::Array{T1, 4}, w::Array{T2, 4}, cdims::ConvDims; kwargs...) where {T1, T2}
    dx = similar(dy, input_size(cdims)..., channels_in(cdims), size(dy, 4))
    return ∇conv_data!(dx, dy, w, cdims; kwargs...)
end

function ∇conv_filter(x::Array{T1, 4}, dy::Array{T2, 4}, cdims::ConvDims; kwargs...) where {T1, T2}
    dw = similar(x, kernel_size(cdims)..., channels_in(cdims), channels_out(cdims))
    return ∇conv_filter!(dw, x, dy, cdims; kwargs...)
end

function maxpool_nnpack!(y::Array{T1, 4}, x::Array{T2, 4}, pdims::PoolDims;
                         kwargs...) where {T1, T2}
    @warn "Automatically converting input tensor to Float32. This will have performance implications" maxlog=1
    # We want the output to be of the same type as desired
    T1.(maxpool_nnpack!(Float32.(y), Float32.(x), pdims; kwargs...))
end

function maxpool_nnpack(x::Array{T, 4}, pdims::PoolDims; kwargs...) where {T}
    y = similar(x, output_size(pdims)..., channels_out(pdims), size(x, 4))
    return maxpool_nnpack!(y, x, pdims; kwargs...)
end

"""
    check_supported_operation(x::Array, pdims::PoolDims)

Returns `true` if nnpack supports the pooling operation for the given input.
"""
function check_supported_operation(x::Array{T, 4}, pdims::PoolDims{2, K, S, P, (1, 1)}) where {T, K, S, P}
    val = size(x)[1:2] .+ (P[1] + P[2], P[3] + P[4]) .- K
    return val .% S == (0, 0) ? true : false
end


text:
using InstantiateFromURL

# activate the QuantEcon environment
activate_github("QuantEcon/QuantEconLectureAllPackages", tag = "v0.9.5");

max_stars_repo_name: fratrik/lecture-source-jl

max_issues_repo_name: fratrik/lecture-source-jl

max_forks_repo_name: fratrik/lecture-source-jl


text:
@safetestset "Properties and Basic Interface" begin
    @safetestset "FirstPoincareInvariant" begin
        using PoincareInvariants

        D = 6; N = 123; θ = canonical_one_form; P = FirstFinDiffPlan
        pinv = FirstPoincareInvariant{Float64, D}(θ, N, P)

        @test getdim(pinv) == D
        @test getplan(pinv) isa P
        @test N <= getpointnum(pinv) <= 2N
        @test getpointspec(pinv) == getpointspec(N, typeof(getplan(pinv)))
        @test getform(pinv) === θ

        pnts = getpoints(pinv) do x
            x, x, x, 0, 1, 2
        end
        @test pnts isa Matrix{Float64}
        @test size(pnts) == (getpointnum(pinv), D)

        @test size(getpoints(pinv)::Vector{Float64}) == (getpointnum(pinv),)
    end

    @safetestset "SecondPoincareInvariant" begin
        using PoincareInvariants

        D = 4; N = 321; ω = canonical_two_form; P = SecondFinDiffPlan
        pinv = SecondPoincareInvariant{Float64, D}(ω, N, P)

        @test getdim(pinv) == D
        @test getplan(pinv) isa P
        @test N <= getpointnum(pinv) <= 2N
        @test getpointspec(pinv) == getpointspec(N, typeof(getplan(pinv)))
        @test getform(pinv) === ω

        pnts = getpoints(pinv) do x, y
            x, x, y, y
        end
        @test pnts isa Matrix{Float64}
        @test size(pnts) == (getpointnum(pinv), D)

        @test size(getpoints(pinv)::Matrix{Float64}) == (getpointnum(pinv), 2)
    end
end

@safetestset "Basic Usage" begin
    @safetestset "FirstPoincareInvariant" begin
        using PoincareInvariants
        using DoubleFloats

        D = 2

        Ts = [
            Float64, Float64, Double64, Float64, Double64,
            Float64, Double64, Double64,
            Double64, Float64, Float64
        ]

        Ns = [
            500, 501, 518, 519, 520,
            541, 542, 543,
            566, 567, 568
        ]

        first_pinvs = [
            FirstPoincareInvariant{Ts[1], D}(canonical_one_form, Ns[1]),
            FirstPoincareInvariant{Ts[2], D}(canonical_one_form, Ns[2], FirstFinDiffPlan),
            FirstPoincareInvariant{Ts[3], D}(canonical_one_form, Ns[3], FirstFinDiffPlan),
            FirstPoincareInvariant{Ts[4], D}(canonical_one_form, Ns[4], FirstFourierPlan),
            FirstPoincareInvariant{Ts[5], D}(canonical_one_form, Ns[5], FirstFourierPlan),

            FirstPI{Ts[6], D}(canonical_one_form, Ns[6]),
            FirstPI{Ts[7], D}(canonical_one_form, Ns[7], FirstFinDiffPlan),
            FirstPI{Ts[8], D}(canonical_one_form, Ns[8], FirstFourierPlan),

            CanonicalFirstPI{Ts[9], D}(Ns[9]),
            CanonicalFirstPI{Ts[10], D}(Ns[10], FirstFinDiffPlan),
            CanonicalFirstPI{Ts[11], D}(Ns[11], FirstFourierPlan),
        ]

        to_circle_line(θ) = (cospi(-2θ), sinpi(-2θ))

        @testset "FirstPoincareInvariant $i" for (i, pinv) in enumerate(first_pinvs)
            T = Ts[i]
            N = Ns[i]

            pnts = getpoints(to_circle_line, pinv)
            @test pnts isa Matrix{T}
            @test size(pnts, 2) == D
            @test size(pnts, 1) == getpointnum(pinv)
            @test N ≤ getpointnum(pinv) ≤ 2N

            I = compute!(pinv, pnts, 0.1, nothing)
            @test I isa T
            @test π ≈ I atol=1e-4
        end
    end

    @safetestset "SecondPoincareInvariant" begin
        using PoincareInvariants
        using DoubleFloats

        D = 2

        Ts = [
            Float64, Float64, Double64, Float64, Double64,
            Float64, Float64, Double64,
            Float64, Double64, Float64,
            Double64, Float64, Float64
        ]

        ps = [
            100, 101, 118, (11, 12), 140,
            162, 173, 196,
            222, 230, 244,
            253, 289, (19, 16)
        ]

        second_pinvs = [
            SecondPoincareInvariant{Ts[1], D}(canonical_two_form, ps[1]),
            SecondPoincareInvariant{Ts[2], D}(canonical_two_form, ps[2], SecondChebyshevPlan),
            SecondPoincareInvariant{Ts[3], D}(canonical_two_form, ps[3], SecondChebyshevPlan),
            SecondPoincareInvariant{Ts[4], D}(canonical_two_form, ps[4], SecondFinDiffPlan),
            SecondPoincareInvariant{Ts[5], D}(canonical_two_form, ps[5], SecondFinDiffPlan),

            SecondPI{Ts[6], D}(canonical_two_form, ps[6]),
            SecondPI{Ts[7], D}(canonical_two_form, ps[7], SecondChebyshevPlan),
            SecondPI{Ts[8], D}(canonical_two_form, ps[8], SecondFinDiffPlan),

            SecondPI{Ts[9], D}(CanonicalSymplecticMatrix{Float64}(D), ps[9]),
            SecondPI{Ts[10], D}(CanonicalSymplecticMatrix{Double64}(D), ps[10], SecondChebyshevPlan),
            SecondPI{Ts[11], D}(CanonicalSymplecticMatrix{Double64}(D), ps[11], SecondFinDiffPlan),

            CanonicalSecondPI{Ts[12], D}(ps[12]),
            CanonicalSecondPI{Ts[13], D}(ps[13], SecondChebyshevPlan),
            CanonicalSecondPI{Ts[14], D}(ps[14], SecondFinDiffPlan)
        ]

        to_wavy_area(x, y) = (x, y + sinpi(2x))

        @testset "SecondPoincareInvariant $i" for (i, pinv) in enumerate(second_pinvs)
            T = Ts[i]
            psi = ps[i]
            N = getpointnum(psi)

            pnts = getpoints(to_wavy_area, pinv)
            @test pnts isa Matrix{T}
            @test size(pnts, 2) == D
            @test size(pnts, 1) == getpointnum(pinv)
            @test N ≤ getpointnum(pinv) ≤ 2N

            I = compute!(pinv, pnts, 0.1, nothing)
            @test I isa T
            @test 1 ≈ I atol=50eps()
        end
    end
end

@safetestset "Consistency Between Implementations" begin
    @safetestset "FirstPoincareInvariant" begin
        using PoincareInvariants

        D = 6
        θ(z, t, p) = [-z[4], -z[5], -p * z[6], t * z[1], p * z[2], z[3]]

        function f(ϕ)
            y, x = sincospi(2ϕ)
            return 0.5 * x, x, 3 * x, 0.1 * y, 0.85 * y, 2 * y
        end

        Idefault = let pinv = FirstPI{Float64, D}(θ, 1_000)
            compute!(pinv, getpoints(f, pinv), 0.3, 0.1)
        end

        Idiff = let pinv = FirstPI{Float64, D}(θ, 1_000, FirstFinDiffPlan)
            compute!(pinv, getpoints(f, pinv), 0.3, 0.1)
        end

        Ifreq = let pinv = FirstPI{Float64, D}(θ, 1_000, FirstFourierPlan)
            compute!(pinv, getpoints(f, pinv), 0.3, 0.1)
        end

        @test Ifreq ≈ Idefault atol=10eps()
        @test Ifreq ≈ Idiff atol=5e-4
    end

    @safetestset "SecondPoincareInvariant" begin
        using PoincareInvariants

        D = 6
        ω(z, t, p) = [
                0     0     0  z[1]  z[2]  z[3]
                0     0     0  z[4]  z[5]  z[6]
                0     0     0     t     p     1
            -z[1] -z[4]    -t     0     0     0
            -z[2] -z[5]    -p     0     0     0
            -z[3] -z[6]    -1     0     0     0
        ]

        f(x, y) = (
            exp(x) * y,
            cospi(x * y),
            x^2 - y^3 + x*y,
            5.0,
            exp(-(x^2 + y^2)),
            x + y
        )

        Idefault = let pinv = SecondPI{Float64, D}(ω, 1_000)
            compute!(pinv, getpoints(f, pinv), 5, 11)
        end

        Icheb = let pinv = SecondPI{Float64, D}(ω, 1_000, SecondChebyshevPlan)
            compute!(pinv, getpoints(f, pinv), 5, 11)
        end

        Ifindiff = let pinv = SecondPI{Float64, D}(ω, (31, 33), SecondFinDiffPlan)
            compute!(pinv, getpoints(f, pinv), 5, 11)
        end

        @test Icheb ≈ Idefault atol=10eps()
        @test Icheb ≈ Ifindiff atol=5e-3
    end
end

@safetestset "Iterable of Points" begin
    using PoincareInvariants

    D = 4
    N = 398
    Nsteps = 13
    times = [2.0 + n * 0.15 for n in 0:Nsteps-1]

    f1(z, t, p) = canonical_one_form(z, t, p) .* t
    pi1 = FirstPI{Float64, D}(f1, N)

    f2(z, t, p) = canonical_two_form(z, t, p) .* t
    pi2 = SecondPI{Float64, D}(f2, N)

    p1(θ) = (cospi(2θ), cospi(2θ), sinpi(2θ), sinpi(2θ))
    p2(x, y) = (x, x, y, y)

    testIs1 = [-2π * t for t in times]
    testIs2 = [2 * t for t in times]

    let data = [getpoints(p1, pi1) for t in times]
        Is = compute!(pi1, data, times, nothing)
        @test testIs1 ≈ Is atol=1e-13
    end

    let data = [getpoints(p2, pi2) for t in times]
        Is = compute!(pi2, data, times, nothing)
        @test testIs2 ≈ Is atol=1e-13
    end

    let data = (getpoints(p1, pi1) for t in times)
        Is = compute!(pi1, data, times, nothing)
        @test testIs1 ≈ Is atol=1e-13
    end

    let data = (getpoints(p2, pi2) for t in times)
        Is = compute!(pi2, data, times, nothing)
        @test testIs2 ≈ Is atol=1e-13
    end
end

@safetestset "Linear Symplectic Maps" begin
    @safetestset "In 2D" begin
        using PoincareInvariants

        diffpi1 = CanonicalFirstPI{Float64, 2}(1_000, FirstFinDiffPlan)
        freqpi1 = CanonicalFirstPI{Float64, 2}(1_000, FirstFourierPlan)

        chebpi2 = CanonicalSecondPI{Float64, 2}(1_000, SecondChebyshevPlan)
        diffpi2 = CanonicalSecondPI{Float64, 2}(1_000, SecondFinDiffPlan)

        A = [1 5;
             0 1]
        B = [2   0;
             0 0.5]

        # circle (negative orientation)
        f1(θ) = A * B * [cospi(2θ), sinpi(2θ)]

        # half circle with radius one
        f2(r, θ) = A * B * [r*cospi(θ), r*sinpi(θ)]

        @test -π ≈ compute!(diffpi1, getpoints(f1, diffpi1), 0, nothing) atol=5e-5
        @test -π ≈ compute!(freqpi1, getpoints(f1, freqpi1), 0, nothing) atol=10eps()

        @test π/2 ≈ compute!(chebpi2, getpoints(f2, chebpi2), 0, nothing) atol=10eps()
        @test π/2 ≈ compute!(diffpi2, getpoints(f2, diffpi2), 0, nothing) atol=5e-3
    end

    @safetestset "In 8D" begin
        using PoincareInvariants

        diffpi1 = CanonicalFirstPI{Float64, 8}(1_000, FirstFinDiffPlan)
        freqpi1 = CanonicalFirstPI{Float64, 8}(1_000, FirstFourierPlan)

        chebpi2 = CanonicalSecondPI{Float64, 8}(1_000, SecondChebyshevPlan)
        diffpi2 = CanonicalSecondPI{Float64, 8}(1_000, SecondFinDiffPlan)

        A = CanonicalSymplecticMatrix(8)

        B = [1 0 0 0 1 5 8 1;
             0 1 0 0 5 2 6 9;
             0 0 1 0 8 6 3 7;
             0 0 0 1 1 9 7 4;
             0 0 0 0 1 0 0 0;
             0 0 0 0 0 1 0 0;
             0 0 0 0 0 0 1 0;
             0 0 0 0 0 0 0 1]

        C = [1.0   0  2.0    0   0   0    0    0;
               0 2.0    0  1.0   0   0    0    0;
               0 1.0    0 -2.0   0   0    0    0;
             2.0   0 -1.0    0   0   0    0    0;
               0   0    0    0 0.2   0  0.4    0;
               0   0    0    0   0 0.4    0  0.2;
               0   0    0    0   0 0.2    0 -0.4;
               0   0    0    0 0.4   0 -0.2    0]

        # test if symplectic
        @test transpose(A * B * C) * (-A) * (A * B * C) ≈ -A

        # half circle with positive orientation
        f1(θ) = A * B * C * [0, cospi(-θ), 0, 0, 0, sinpi(-θ), 0, 0]

        # half circle with radius one and negative orientation
        f2(r, θ) = A * B * C * [0, r*cospi(θ), 0, 0, 0, -r*sinpi(θ), 0, 0]

        @test π/2 ≈ compute!(diffpi1, getpoints(f1, diffpi1), 0, nothing) atol=1e-2
        @test π/2 ≈ compute!(freqpi1, getpoints(f1, freqpi1), 0, nothing) atol=5e-3

        @test -π/2 ≈ compute!(chebpi2, getpoints(f2, chebpi2), 0, nothing) atol=10eps()
        @test -π/2 ≈ compute!(diffpi2, getpoints(f2, diffpi2), 0, nothing) atol=5e-3
    end
end

@safetestset "Henon Map" begin
    using PoincareInvariants

    init1(θ) = (cospi(2θ), cospi(2θ), 0, -sinpi(2θ))
    I1 = π

    init2(x, y) = 2 .* (x - 0.5, -x + 0.5, y - 0.5, -y + 0.5)
    I2 = 2^2 * 2

    # Area preserving Henon map
    function henon(q, p)
        a = 1.4; b = 0.3
        q, p = q, 1 + p - a * q^2 # fold
        q, p = b * q, p / b # press
        q, p = -p, q  # rotate
        return q, p
    end

    # henon plus symplectic intermixing of q1/p1 with q2/p2
    function f((q1, q2, p1, p2))
        q1, p1 = henon(q1, p1)
        q2, p2 = henon(q2, p2)
        return (q1, q2, p1 - q2, p2 - q1)
    end

    diffpi1 = CanonicalFirstPI{Float64, 4}(1_000, FirstFinDiffPlan)
    freqpi1 = CanonicalFirstPI{Float64, 4}(1_000, FirstFourierPlan)

    chebpi2 = CanonicalSecondPI{Float64, 4}(1_000, SecondChebyshevPlan)
    diffpi2 = CanonicalSecondPI{Float64, 4}(1_000, SecondFinDiffPlan)

    f1(x) = init1(x) |> f
    @test I1 ≈ compute!(diffpi1, getpoints(f1, diffpi1), 0, nothing) atol=5e-5
    @test I1 ≈ compute!(freqpi1, getpoints(f1, freqpi1), 0, nothing) atol=5e-15

    g1(x, y) = init2(x, y) |> f
    @test I2 ≈ compute!(chebpi2, getpoints(g1, chebpi2), 0, nothing) atol=1e-14
    @test I2 ≈ compute!(diffpi2, getpoints(g1, diffpi2), 0, nothing) atol=5e-13

    f2(x) = init1(x) |> f |> f
    @test I1 ≈ compute!(diffpi1, getpoints(f2, diffpi1), 0, nothing) atol=5e-5
    @test I1 ≈ compute!(freqpi1, getpoints(f2, freqpi1), 0, nothing) atol=1e-13

    g2(x, y) = init2(x, y) |> f |> f
    @test I2 ≈ compute!(chebpi2, getpoints(g2, chebpi2), 0, nothing) atol=1e-11
    @test I2 ≈ compute!(diffpi2, getpoints(g2, diffpi2), 0, nothing) atol=5e-12

    f3(x) = init1(x) |> f |> f |> f
    @test I1 ≈ compute!(freqpi1, getpoints(f3, freqpi1), 0, nothing) atol=1e-7

    g3(x, y) = init2(x, y) |> f |> f |> f
    @test I2 ≈ compute!(chebpi2, getpoints(g3, chebpi2), 0, nothing) atol=5e-6
end

@safetestset "Wavy Map" begin
    using PoincareInvariants

    function wavy(x, y)
        x, y = x, y + 2 * sinpi(x)
        x, y = -y, x
        return x, y
    end

    # make wavy and mix dims
    function f((x1, x2, x3, y1, y2, y3))
        x1, y1 = wavy(x1, y1)
        x2, y2 = wavy(x2, y2)
        x3, y3 = wavy(x3, y3)
        return (x1 + y2, x2 + y1 + y3, x3 + y2, y1, y2, y3)
    end

    init1(θ) = (cospi(2θ), sinpi(2θ), cospi(2θ), sinpi(2θ), 1, sinpi(2θ))
    I1 = -2π

    init2(x, y) = (2x, 2x-1, 2x-2, 2y-1, 2y-1, 2y-1)
    I2 = 3 * 2^2

    diffpi1 = CanonicalFirstPI{Float64, 6}(1_000, FirstFinDiffPlan)
    freqpi1 = CanonicalFirstPI{Float64, 6}(1_000, FirstFourierPlan)

    chebpi2 = CanonicalSecondPI{Float64, 6}(1_000, SecondChebyshevPlan)
    diffpi2 = CanonicalSecondPI{Float64, 6}(1_000, SecondFinDiffPlan)

    f1(x) = init1(x) |> f
    @test I1 ≈ compute!(diffpi1, getpoints(f1, diffpi1), 0, nothing) atol=5e-5
    @test I1 ≈ compute!(freqpi1, getpoints(f1, freqpi1), 0, nothing) atol=1e-15

    g1(x, y) = init2(x, y) |> f
    @test I2 ≈ compute!(chebpi2, getpoints(g1, chebpi2), 0, nothing) atol=1e-14
    @test I2 ≈ compute!(diffpi2, getpoints(g1, diffpi2), 0, nothing) atol=5e-15

    f2(x) = init1(x) |> f |> f
    @test I1 ≈ compute!(diffpi1, getpoints(f2, diffpi1), 0, nothing) atol=1e-3
    @test I1 ≈ compute!(freqpi1, getpoints(f2, freqpi1), 0, nothing) atol=5e-14

    g2(x, y) = init2(x, y) |> f |> f
    @test I2 ≈ compute!(chebpi2, getpoints(g2, chebpi2), 0, nothing) atol=5e-3
    @test I2 ≈ compute!(diffpi2, getpoints(g2, diffpi2), 0, nothing) atol=1e-3

    f3(x) = init1(x) |> f |> f |> f
    @test I1 ≈ compute!(diffpi1, getpoints(f3, diffpi1), 0, nothing) atol=5e-2
    @test I1 ≈ compute!(freqpi1, getpoints(f3, freqpi1), 0, nothing) atol=1e-13

    g3(x, y) = init2(x, y) |> f |> f |> f
    # neither of the two methods manage this
end

@safetestset "Pendulum" begin
    using PoincareInvariants

    # m, g, L = 1
    function update(points::AbstractMatrix, dt)
        out = copy(points)
        for pnt in eachrow(out)
            pnt[1] += dt * pnt[2]
            pnt[2] -= dt * sin(pnt[1])
        end
        out
    end

    integrate(initial::AbstractMatrix, dt, n) = accumulate(1:n, init=initial) do points, _
        update(points, dt)
    end

    init1(θ) = 2 .* sincospi(2θ)
    I1 = 4π

    init2(x, y) = (π * x - π/2, 4 * y - 2)
    I2 = 4π

    pi1 = CanonicalFirstPI{Float64, 2}(1_000)
    pi2 = CanonicalSecondPI{Float64, 2}(1_00)

    dt = 0.05
    tmax = 2
    times = 0:dt:tmax

    data1 = integrate(getpoints(init1, pi1), dt, length(times))
    data2 = integrate(getpoints(init2, pi2), dt, length(times))

    @test maximum(abs, I1 .- compute!(pi1, data1, times, nothing)) < 1e-14
    @test maximum(abs, I2 .- compute!(pi2, data2, times, nothing)) < 1e-3
end


text:
using Luxor

isinside(p, a, b) = (b.x - a.x) * (p.y - a.y) > (b.y - a.y) * (p.x - a.x)

function intersection(a, b, s, f)
    dc = [a.x - b.x, a.y - b.y]
    dp = [s.x - f.x, s.y - f.y]
    n1 = a.x * b.y - a.y * b.x
    n2 = s.x * f.y - s.y * f.x
    n3 = 1.0 / (dc[1] * dp[2] - dc[2] * dp[1])
    Point((n1 * dp[1] - n2 * dc[1]) * n3, (n1 * dp[2] - n2 * dc[2]) * n3)
end

function clipSH(spoly, cpoly)
    outarr = spoly
    q = cpoly[end]
    for p in cpoly
        inarr = outarr
        outarr = Point[]
        s = inarr[end]
        for vtx in inarr
            if isinside(vtx, q, p)
                if !isinside(s, q, p)
                    push!(outarr, intersection(q, p, s, vtx))
                end
                push!(outarr, vtx)
            elseif isinside(s, q, p)
                push!(outarr, intersection(q, p, s, vtx))
            end
            s = vtx
        end
        q = p
    end
    outarr
end

subjectp = [Point(50, 150), Point(200, 50), Point(350, 150), Point(350, 300),
    Point(250, 300), Point(200, 250), Point(150, 350), Point(100, 250), Point(100, 200)]

clipp = [Point(100, 100), Point(300, 100), Point(300, 300), Point(100, 300)]

Drawing(400, 400, "intersecting-polygons.png")
background("white")
sethue("red")
poly(subjectp, :stroke, close=true)
sethue("blue")
poly(clipp, :stroke, close=true)
clipped = clipSH(subjectp, clipp)
sethue("gold")
poly(clipped, :fill, close=true)
finish()
preview()
println(clipped)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import PowerModelsACDC;
const _PMACDC = PowerModelsACDC;
import PowerModels;
const _PM = PowerModels;
import InfrastructureModels;
const _IM = InfrastructureModels;
import JuMP
import Gurobi
using MAT
using XLSX
using JLD2
using Statistics, FileIO

include("basencont_nw.jl")
Total_sample = 10# samples per yearsum()
total_yr = 3# the years in horizon, data coming from excels
period = "multi" # single or multi
# Prot_system = "NS_CB"
# Prot_system_coll = ["NS_CB"]
Prot_system_coll = ["FS_HDCCB", "NS_CB"]
curtailed_gen = [1,2] #geneartor numbers # change also constraint max(), generating power,file name
syncarea = 2
max_curt = 1
cost = Any[]
for proti = 1:3
    Prot_system = Prot_system_coll[proti]
    file = "./test/data/4bus_TNEP.m"
    data_sp = _PM.parse_file(file)
    _PMACDC.process_additional_data!(data_sp)

    include("bkgrnd_cal.jl")
    gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "Presolve" => -1)

    no_nw = (length(data_sp["branchdc_ne"]) + 1) * Total_sample*total_yr
    data_mp = multi_network(file, no_nw)
    # year = ["NAT_2025_Generation", "NAT_2030_Generation", "NAT_2035_Generation", "NAT_2040_Generation", "NAT_2045_Generation", "NAT_2050_Generation"]
    # year_num = ["2025", "2030", "2035", "2040", "2045", "2050"]    # year = ["NAT_2045_Generation"]

    year = ["NAT_2025_Generation","NAT_2035_Generation", "NAT_2045_Generation"]
    year_num = ["2025", "2035", "2045"]
    # year_num = ["2035"]
    # year = ["NAT_2035_Generation"]
    yr = 1
    # @load "scenario.jld2"

    data_ps, cost = mp_datainputs_ne(deepcopy(data_mp), Total_sample, year, year_num,file)
    base_weight_ip = Array{Float64}(undef, 1, no_nw)
    @save "costt.jld2" cost
    @load "cluster.jld2"
    data_cont, Cont_list, base_list, base_weight = mp_contignecy_ne(deepcopy(data_ps), Total_sample*total_yr, base_weight_ip, weights)
    year_base = year_base_networks(Total_sample, total_yr, base_list)

    # year_base = year_base_networks(Total_sample, total_yr, base_list)
        if Prot_system == "FS_HDCCB" || Prot_system == "FS_MDCCB"
        s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => true, "NSprotection" => false,
        "Permanentloss" => false, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea, "weights" => base_weight, "year_base" => year_base, "total_yr" => total_yr)
        elseif Prot_system == "NS_CB" || Prot_system == "NS_FB"
        s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => false, "NSprotection" => true,
        "Permanentloss" => false, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea, "weights" => base_weight, "year_base" => year_base, "total_yr" => total_yr)
        else
        s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => false, "NSprotection" => false,
        "Permanentloss" => true, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea, "weights" => base_weight, "year_base" => year_base, "total_yr" => total_yr)
        end

    if period == "single"   push!(s, "multiperiod"=> false)
    elseif period == "multi"  push!(s, "multiperiod"=> true )
    end

    @assert ( s["FSprotection"] == true && (Prot_system == "FS_HDCCB" || Prot_system == "FS_MDCCB") ) || (s["NSprotection"] == true && (Prot_system == "NS_CB" || Prot_system == "NS_FB"))|| (Prot_system == "Permanentloss" && s["NSprotection"] == false && s["FSprotection"] == false)

    if occursin("4bus", file)
        filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\TNEP\\4bus\\",Prot_system,"_poscl.xlsx")
    elseif occursin("6bus", file)
        filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\TNEP\\6bus\\",Prot_system,"_poscl.xlsx")
    end

    for (n,nw) in data_cont["nw"]
                for (r,reserves) in nw["reserves"]
                     if Prot_system == "FS_HDCCB"; reserves["Tcl"] = 0.148
                     elseif Prot_system == "FS_MDCCB"; reserves["Tcl"] = 0.293
                     elseif Prot_system == "NS_FB"; reserves["Tcl"] = 0.100
                     elseif Prot_system == "NS_CB"; reserves["Tcl"] = 0.150
                     elseif Prot_system == "Permanentloss"; reserves["Tcl"] = 100
                     end
                 end
    end

    data_cont1 = Protectionsystemcost_4bus(deepcopy(data_cont), Prot_system, no_nw)
    resultDC1 = _PMACDC.run_mp_tnepscopf(data_cont1, _PM.DCPPowerModel, gurobi, multinetwork = true;  setting = s)
    display_keyindictionary(resultDC1, "isbuilt", "Pgg", base_list)
    built_cv, built_br = _PMACDC.display_results_tnep_mp(resultDC1)
    # curtail, maxFFR = curtailment(data_cont1, base_list, resultDC1, curtailed_gen)
     XLSX.openxlsx(filepath, mode="w") do xf
        sheet = xf[1]
        sheet["$(string("A",1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
        sheet["$(string("A",2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
        sheet["$(string("A",3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
        sheet["$(string("A",4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
        sheet["$(string("A",5))"] = mean(resultDC1["solution"]["nw"]["1"]["Curt"])
        # sheet["$(string("A",5))"] = sum(resultDC1["solution"]["nw"]["1"]["Cont"])
        sheet["$(string("A",6))"] = resultDC1["objective"]
        sheet["$(string("A",7))"] = resultDC1["solution"]["nw"]["1"]["Cont"]
        sheet["$(string("A",8))"] = 0
        sheet["$(string("A",9))"] = resultDC1["solution"]["nw"]["1"]["Inv_cost"]
        sheet["$(string("A",10))"] = built_cv
        sheet["$(string("A",11))"] = built_br
        sheet["$(string("A",12))"] = resultDC1["objective_lb"]
    end
end
######################start inertia############################
# inertia = 1:0.5:4
# column = ["A" "B" "C" "D" "E" "F" "G" "H"]
# filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\inertia\\",Prot_system,".xlsx")
# XLSX.openxlsx(filepath, mode="w") do xf
#     for i = 1:length(inertia)
#         for (n,nw) in data_cont["nw"]
#             for (r,reserves) in nw["reserves"]
#                  display(reserves)
#                  if Prot_system == "FS_HDCCB"; reserves["Tcl"] = 0.148
#                  elseif Prot_system == "FS_MDCCB"; reserves["Tcl"] = 0.293
#                  elseif Prot_system == "NS_FB"; reserves["Tcl"] = 0.100
#                  elseif Prot_system == "NS_CB"; reserves["Tcl"] = 0.150
#                  elseif Prot_system == "Permanentloss"; reserves["Tcl"] = 100
#                  end
#                  reserves["H"] = inertia[i]
#             end
#         end
#          display(data_cont["nw"]["1"]["gen"]["3"]["pmax"])
#          data_cont["nw"]["1"]["gen"]["1"]["pmin"] = data_cont["nw"]["1"]["gen"]["1"]["pmax"]
#          data_cont["nw"]["1"]["gen"]["2"]["pmin"] = data_cont["nw"]["1"]["gen"]["2"]["pmax"]
#          resultDC1 = _PMACDC.run_acdcscopf(data_cont, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
#          curtail, maxFFR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
#          sheet = xf[1]
#          cell_no = string(column[i],i)
#          sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
#          sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
#          sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
#          sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
#          sheet["$(string(column[i],5))"] = sum(curtail)
#          sheet["$(string(column[i],6))"] = resultDC1["objective"]
#          sheet["$(string(column[i],7))"] = maxFFR
#         end
# end
######################end inertia############################

#####################start fuel cost############################
# fuel_cost = 8:1:12
# column = ["A" "B" "C" "D" "E" "F" "G" "H"]
# filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\6bus\\fuel_cost\\",Prot_system,".xlsx", )
# XLSX.openxlsx(filepath, mode = "w") do xf
#         for i = 1:length(fuel_cost)
#             for (n, nw) in data_cont["nw"]
#                 nw["gen"]["3"]["cost"][3] = fuel_cost[i] * 100
#             end
#             resultDC1 = _PMACDC.run_acdcscopf(data_cont, _PM.DCPPowerModel, gurobi, multinetwork = true; setting = s)
#             curtail, maxFFR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
#             cell_no = string(column[i],i)
#             sheet = xf[1]
#             sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
#             sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
#             sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
#             sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
#             sheet["$(string(column[i],5))"] = sum(curtail)
#             sheet["$(string(column[i],6))"] = resultDC1["objective"]
#             sheet["$(string(column[i],7))"] = maxFFR
#             sheet["$(string(column[i],8))"] = resultDC1["objective_lb"]
#         end
#     end
# end
######################end fuel cost############################

######################start FFR cost############################
# FFR_cost = 64:8:96
# column = ["A" "B" "C" "D" "E" "F" "G" "H"]
# filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\6bus\\FFR_cost\\",Prot_system,".xlsx")
# XLSX.openxlsx(filepath, mode="w") do xf
#     for i = 1:length(FFR_cost)
#         for (n,nw) in data_cont["nw"]
#             for (r,reserves) in nw["reserves"]
#            reserves["Cf"] = FFR_cost[i]
#             end
#         end
#             resultDC1 = _PMACDC.run_acdcscopf(data_cont, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
#             curtail, maxFFR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
#             cell_no = string(column[i],i)
#             sheet = xf[1]
#             sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
#             sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
#             sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
#             sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
#             sheet["$(string(column[i],5))"] = sum(curtail)
#             sheet["$(string(column[i],6))"] = resultDC1["objective"]
#             sheet["$(string(column[i],7))"] = maxFFR
#             sheet["$(string(column[i],8))"] = resultDC1["objective_lb"]
#         end
# end
# end
######################end FFR cost############################

######################start FCR cost############################
# FCR_cost = [4.24 4.77 5.3 5.83 6.36]
# column = ["A" "B" "C" "D" "E" "F" "G" "H"]
# filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\6bus\\FCR_cost\\",Prot_system,".xlsx")
# XLSX.openxlsx(filepath, mode="w") do xf
#     for i = 1:length(FCR_cost)
#         for (n,nw) in data_cont["nw"]
#             for (r,reserves) in nw["reserves"]
#                  reserves["Cg"] = FCR_cost[i]
#             end
#         end
#               resultDC1 = _PMACDC.run_acdcscopf(data_cont, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
#               curtail, maxFFR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
#               cell_no = string(column[i],i)
#               sheet = xf[1]
#               sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
#               sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
#               sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
#               sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
#               sheet["$(string(column[i],5))"] = sum(curtail)
#               sheet["$(string(column[i],6))"] = resultDC1["objective"]
#               sheet["$(string(column[i],7))"] = maxFFR
#               sheet["$(string(column[i],8))"] = resultDC1["objective_lb"]
#     end
# end
# end
######################end FCR cost############################

######################start FFR time############################
# FFR_time = [0.1 0.5 1 1.5]
# column = ["A" "B" "C" "D" "E" "F" "G" "H"]
# filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\6bus\\FFR_time\\",Prot_system,".xlsx")
# XLSX.openxlsx(filepath, mode="w") do xf
#         for i = 1:length(FFR_time)
#             for (n,nw) in data_cont["nw"]
#                 for (r,reserves) in nw["reserves"]
#                     reserves["Tf"] = FFR_time[i]
#                 end
#             end
#              resultDC1 = _PMACDC.run_acdcscopf(data_cont, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
#              curtail, maxFFR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
#              sheet = xf[1]
#              cell_no = string(column[i],i)
#              sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
#              sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
#              sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
#              sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
#              sheet["$(string(column[i],5))"] = sum(curtail)
#              sheet["$(string(column[i],6))"] = resultDC1["objective"]
#              sheet["$(string(column[i],7))"] = maxFFR
#              sheet["$(string(column[i],8))"] = resultDC1["objective_lb"]
#         end
# end
# end
######################end FFR time############################

######################start FCR time############################
# FCR_time = [5  10]
# column = ["A" "B" "C" "D" "E" "F" "G" "H"]
# filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\6bus\\FCR_time\\",Prot_system,".xlsx")
# XLSX.openxlsx(filepath, mode="w") do xf
#         for i = 1:length(FCR_time)
#             for (n,nw) in data_cont["nw"]
#                 for (r,reserves) in nw["reserves"]
#                      reserves["Tg"] = FCR_time[i]
#                 end
#             end
#          resultDC1 = _PMACDC.run_acdcscopf(data_cont, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
#          curtail, maxFFR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
#          sheet = xf[1]
#                  cell_no = string(column[i],i)
#                  sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
#                  sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
#                  sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
#                  sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
#                  sheet["$(string(column[i],5))"] = sum(curtail)
#                  sheet["$(string(column[i],6))"] = resultDC1["objective"]
#                  sheet["$(string(column[i],7))"] = maxFFR
#                  sheet["$(string(column[i],8))"] = resultDC1["objective_lb"]
#         end
# end
# end
######################end FCR time############################

######################start Curtailemnt percentage############################
#     curtl = [0 0.1 0.2 0.3]
#     column = ["A" "B" "C" "D" "E" "F" "G" "H"]
#     filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\6bus\\Curtailment\\",Prot_system,".xlsx")
#     XLSX.openxlsx(filepath, mode="w") do xf
#              for i = 1:length(curtl)
#                 s["max_curt"] = curtl[i]
#                 resultDC1 = _PMACDC.run_acdcscopf(data_cont, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
#                 curtail, maxFFR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
#                 sheet = xf[1]
#                 cell_no = string(column[i],i)
#                 sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
#                 sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
#                 sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
#                 sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
#                 sheet["$(string(column[i],5))"] = sum(curtail)
#                 sheet["$(string(column[i],6))"] = resultDC1["objective"]
#                 sheet["$(string(column[i],7))"] = maxFFR
#                 sheet["$(string(column[i],8))"] = resultDC1["objective_lb"]
#             end
#    end
# end

# year = ["NAT_2025_Generation", "NAT_2030_Generation", "NAT_2035_Generation", "NAT_2040_Generation", "NAT_2045_Generation", "NAT_2050_Generation"]
# yr = 1
# nw_no = 0
# fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[1],"\\output.mat")
# vars = matread(fname)
# fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[1],"\\output_MC.mat")
# vars1 = matread(fname1)
# fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
# vars2 = matread(fname2)

# function mp_datainputs(data_ip,Total_sample)
#     kk= 1
#     for i = 1:Total_sample
#         sample = Int(round(8760*(1-rand())))
#         for tt = 1:(length(data_sp["branchdc_ne"])+1)
#             display("kk:$kk")
#             display("sample:$sample")
#             data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = vars["Heq"][sample]
#             data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2] = vars1["MC1"][sample]*100
#             rating = deepcopy(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
#             data_ip["nw"]["$kk"]["gen"]["1"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
#             data_ip["nw"]["$kk"]["gen"]["2"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
#             what about load?
#             kk +=1
#         end
#     end
#     return data_ip
# end

################multiperiod function###################check with Hakan
# extradata = mp_datainputs(data_sp,100)
# function mp_datainputs(data_sp,Total_sample)
#     extradata = Dict{String,Any}()
#     extradata["dim"] = 100
#     extradata["gen"] = Dict{String,Any}()
#     extradata["reserves"] = Dict{String,Any}()
#
#     extradata["gen"]["1"] = Dict{String,Any}()
#     extradata["gen"]["2"] = Dict{String,Any}()
#     extradata["gen"]["3"] = Dict{String,Any}()
#     extradata["reserves"]["2"] = Dict{String,Any}()
#
#     dim = Total_sample
#     extradata["gen"]["1"]["pmax"] = Array{Float64,2}(undef, 1, dim)
#     extradata["gen"]["2"]["pmax"] = Array{Float64,2}(undef, 1, dim)
#     extradata["gen"]["3"]["cost"] = Array{Float64,2}(undef, 2, dim)
#     extradata["reserves"]["2"]["H"] = Array{Float64,2}(undef, 1, dim)
#
#         f_index = 1
#         while f_index <= dim
#             sample = Int(round(8760*(1-rand())))
#             rating = deepcopy(data_sp["gen"]["1"]["pmax"])
#             extradata["gen"]["1"]["pmax"][1, f_index] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
#             extradata["gen"]["2"]["pmax"][1, f_index] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
#             extradata["gen"]["3"]["cost"][2, f_index] =  vars1["MC1"][sample]*100
#             extradata["reserves"]["2"]["H"][1, f_index] = vars["Heq"][sample]
#             f_index  = f_index  + 1
#         end
#         return extradata
# end
# mn_data = PowerModelsACDC.multinetwork_data(data_sp, extradata, Set{String}(["source_type", "name", "source_version", "per_unit"]))

# fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[1],"\\output.mat")

############curtailment##############################

# curtail = curtailment(data_cont, base_list)


text:
using Dates

struct SeriesInfo
    id::String
    category::String
    freq::UInt16
    period::String
    horizon::UInt16
    start_time::DateTime
end

const SeriesValue = Float32

const Series = Vector{SeriesValue}

const SeriesDict = Dict{String, Series}

const SeriesInfoDict = Dict{String, SeriesInfo}

function read_series(file_path)
    data = SeriesDict()
    is_header = true

    for line in eachline(file_path)
        if is_header
            is_header = false
        else
            parts = split(line, ",")
            parts = map(p->p[2:end - 1], parts)
            id = parts[1]
            values = [parse(Float64, p) for p in parts[2:end] if !isempty(p)]
            data[id] = values
        end
    end

    data
end

function parse_time(value, formats)
    for format in formats
        try
            return DateTime(value, format)
        catch
            continue
        end
    end
end

function read_info(file_path)
    info_dict = SeriesInfoDict()
    time_formats = [DateFormat("d-m-y H:M"), DateFormat("y-m-d H:M:S")]
    is_header = true

    for line in eachline(file_path)
        if is_header
            is_header = false
        else
            parts = split(line, ",")
            id = parts[1]
            category = parts[2]
            freq = parse(Int8, parts[3])
            period = parts[5]
            horizon = parse(Int8, parts[4])
            start_time = parse_time(parts[6], time_formats)
            info_dict[id] = SeriesInfo(id, category, freq, period, horizon, start_time)
        end
    end

    info_dict
end