text:
function read_files_basecase(PAROTS_flag, raw_path, rop_path, inl_path, con_path)
    if (PAROTS_flag) println("Reading input DATA: "); end
    contingency=zeros(Int32,0);
    (network, bus_shift)=process_raw(PAROTS_flag, raw_path);
    network=process_inl(PAROTS_flag, inl_path, network);
    network=process_rop(PAROTS_flag, rop_path, network, bus_shift);
    #contingency=process_con(PAROTS_flag, con_path, network, bus_shift);
    #network=add_LODFs(network, contingency)
    return(network, contingency);
end

function read_files_contingency(PAROTS_flag, raw_path, rop_path, inl_path, con_path)
    if (PAROTS_flag) println("Reading input DATA: "); end
    (network, bus_shift)=process_raw(PAROTS_flag, raw_path);
    network=process_inl(PAROTS_flag, inl_path, network);
    #network=process_rop(rop_path, network, bus_shift);
    contingency=process_con(PAROTS_flag, con_path, network, bus_shift);
    return(network, contingency);
end

"--------------------------------------------------
      Readiang 'raw' file (Network Data File)
---------------------------------------------------"
function add_buses(net_bus, RawString, start, size)
    bus_index=Dict{String,Any}();
    for i=1:size
        net_bus=merge(net_bus, Dict(string(i)=>Dict{String,Any}()));
        bus=split(RawString[i+start-1], ",");
        I=parse(Int64, split(bus[1]," "; limit=0, keepempty=false)[1]);
        VM=parse(Float64, split(bus[8]," "; limit=0, keepempty=false)[1]);
        VA=parse(Float64, split(bus[9]," "; limit=0, keepempty=false)[1]);
        NVHI=parse(Float64, split(bus[10]," "; limit=0, keepempty=false)[1]);
        NVLO=parse(Float64, split(bus[11]," "; limit=0, keepempty=false)[1]);
        EVHI=parse(Float64, split(bus[12]," "; limit=0, keepempty=false)[1]);
        EVLO=parse(Float64, split(bus[13]," "; limit=0, keepempty=false)[1]);
        bus_index=merge(bus_index, Dict(string(I)=>i));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("index"=>I));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("gen"=>false));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("vm"=>VM));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("va"=>deg2rad(VA)));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("vmax"=>NVHI));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("vmin"=>NVLO));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("evmax"=>EVHI));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("evmin"=>EVLO));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("bus_type"=>1));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("Pd"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("Qd"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("GL"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("BL"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("sh_max"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("sh_min"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("bcs"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("gfM"=>0.0));
        net_bus[string(i)]=merge(net_bus[string(i)], Dict("bfM"=>0.0));
    end
    return(net_bus, bus_index)
end
function add_loads(net_bus, bus_index, BaseMVA, RawString, start, size)
    for i=1:size
        load=split(RawString[i+start-1], ",");
        I_index=parse(Int64, split(load[1]," "; limit=0, keepempty=false)[1]);
        I=bus_index[string(I_index)];
        ST=parse(Int64, split(load[3]," "; limit=0, keepempty=false)[1]);
        PL=ST*parse(Float64, split(load[6]," "; limit=0, keepempty=false)[1])/BaseMVA+net_bus[string(I)]["Pd"];
        QL=ST*parse(Float64, split(load[7]," "; limit=0, keepempty=false)[1])/BaseMVA+net_bus[string(I)]["Qd"];
        net_bus[string(I)]=merge(net_bus[string(I)], Dict("Pd"=>PL));
        net_bus[string(I)]=merge(net_bus[string(I)], Dict("Qd"=>QL));
    end
    return(net_bus)
end

function add_shunts(net_bus, bus_index, BaseMVA, RawString, start, size)
    for i=1:size
        shunt=split(RawString[i+start-1], ",");
        I_index=parse(Int64, split(shunt[1]," "; limit=0, keepempty=false)[1]);
        I=bus_index[string(I_index)];
        ST=parse(Int64, split(shunt[3]," "; limit=0, keepempty=false)[1]);
        GL=ST*parse(Float64, split(shunt[4]," "; limit=0, keepempty=false)[1])/BaseMVA+net_bus[string(I)]["GL"];
        BL=ST*parse(Float64, split(shunt[5]," "; limit=0, keepempty=false)[1])/BaseMVA+net_bus[string(I)]["BL"];
        net_bus[string(I)]=merge(net_bus[string(I)], Dict("GL"=>GL));
        net_bus[string(I)]=merge(net_bus[string(I)], Dict("BL"=>BL));
    end
    return(net_bus)
end
function add_gens(net_bus, net_gen, bus_index, BaseMVA, RawString, start, size)
    temp_pamx_id=1; gen_buses=zeros(Int32,0);
    for i=1:size
        net_gen=merge(net_gen, Dict(string(i)=>Dict{String,Any}()));
        gen=split(RawString[i+start-1], ",");
        I_index=parse(Int64, split(gen[1]," "; limit=0, keepempty=false)[1]);
        I=bus_index[string(I_index)];
        ID=split(gen[2],r"\W"; limit=0, keepempty=false);
        if (length(ID)>0) ID=split(ID[1]," "; limit=0, keepempty=false)[1]; end
        St=parse(Int64, split(gen[15]," "; limit=0, keepempty=false)[1]);
        PG=parse(Float64, split(gen[3]," "; limit=0, keepempty=false)[1])/BaseMVA;
        QG=parse(Float64, split(gen[4]," "; limit=0, keepempty=false)[1])/BaseMVA;
        QT=parse(Float64, split(gen[5]," "; limit=0, keepempty=false)[1])/BaseMVA;
        QB=parse(Float64, split(gen[6]," "; limit=0, keepempty=false)[1])/BaseMVA;
        PT=parse(Float64, split(gen[17]," "; limit=0, keepempty=false)[1])/BaseMVA;
        PB=parse(Float64, split(gen[18]," "; limit=0, keepempty=false)[1])/BaseMVA;

        net_gen[string(i)]=merge(net_gen[string(i)], Dict("gen_bus_index"=>I));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("gen_bus"=>I_index));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("id"=>ID));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("gen_status"=>St));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("pg"=>PG));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("qg"=>QG));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("pmax"=>PT));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("qmax"=>QT));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("pmin"=>PB));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("qmin"=>QB));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("index"=>i));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("gov_droop"=>0.0));
        net_gen[string(i)]=merge(net_gen[string(i)], Dict("total_apf"=>0.0));

        if (St==1) net_bus[string(I)]["bus_type"]=2; net_bus[string(I)]["gen"]=true; append!(gen_buses,I); end
        if (St*PT>net_gen[string(temp_pamx_id)]["pmax"]) temp_pamx_id=i; end
    end
    net_bus[string(net_gen[string(temp_pamx_id)]["gen_bus_index"])]["bus_type"]=3;
    return(net_bus, net_gen, net_gen[string(temp_pamx_id)]["gen_bus_index"], gen_buses);
end

function add_Lbranches(net_br, bus_index, BaseMVA, RawString, start, size)
    i=0;
    for ind=1:size
        Lbr=split(RawString[ind+start-1], ",");
        St=parse(Int64, split(Lbr[14]," "; limit=0, keepempty=false)[1]);
        if (St==1) i+=1;
            net_br=merge(net_br, Dict(string(i)=>Dict{String,Any}()));
            # println("string_i: ", string(i))
            I_index=parse(Int64, split(Lbr[1]," "; limit=0, keepempty=false)[1]);
            I=bus_index[string(I_index)];
            J_index=parse(Int64, split(Lbr[2]," "; limit=0, keepempty=false)[1]);
            J=bus_index[string(J_index)];
            CKT=split(Lbr[3],r"\W"; limit=0, keepempty=false);
            if (length(CKT)>0) CKT=split(CKT[1]," "; limit=0, keepempty=false)[1]; end
            R=parse(Float64, split(Lbr[4]," "; limit=0, keepempty=false)[1]);
            X=parse(Float64, split(Lbr[5]," "; limit=0, keepempty=false)[1]);
            B=parse(Float64, split(Lbr[6]," "; limit=0, keepempty=false)[1]);
            RATEA=parse(Float64, split(Lbr[7]," "; limit=0, keepempty=false)[1])/BaseMVA;
            RATEC=parse(Float64, split(Lbr[9]," "; limit=0, keepempty=false)[1])/BaseMVA;

            net_br[string(i)]=merge(net_br[string(i)], Dict("br_r"=>R));
            net_br[string(i)]=merge(net_br[string(i)], Dict("ckt"=>CKT));
            net_br[string(i)]=merge(net_br[string(i)], Dict("br_status"=>St));
            net_br[string(i)]=merge(net_br[string(i)], Dict("br_x"=>X));
            net_br[string(i)]=merge(net_br[string(i)], Dict("b_fr"=>B/2));
            net_br[string(i)]=merge(net_br[string(i)], Dict("b_to"=>B/2));
            net_br[string(i)]=merge(net_br[string(i)], Dict("rate_a"=>RATEA));
            net_br[string(i)]=merge(net_br[string(i)], Dict("rate_c"=>RATEC));
            net_br[string(i)]=merge(net_br[string(i)], Dict("f_bus"=>I_index));
            net_br[string(i)]=merge(net_br[string(i)], Dict("t_bus"=>J_index));
            net_br[string(i)]=merge(net_br[string(i)], Dict("f_bus_index"=>I));
            net_br[string(i)]=merge(net_br[string(i)], Dict("t_bus_index"=>J));
            net_br[string(i)]=merge(net_br[string(i)], Dict("transformer"=>false));
            net_br[string(i)]=merge(net_br[string(i)], Dict("tap"=>1.0));
            net_br[string(i)]=merge(net_br[string(i)], Dict("shift"=>0.0));
            net_br[string(i)]=merge(net_br[string(i)], Dict("index"=>i));
            net_br[string(i)]=merge(net_br[string(i)], Dict("gfM"=>0.0));
            net_br[string(i)]=merge(net_br[string(i)], Dict("bfM"=>0.0));
            net_br[string(i)]=merge(net_br[string(i)], Dict("g_fr"=>0.0));
            net_br[string(i)]=merge(net_br[string(i)], Dict("g_to"=>0.0));
        end
    end
    return(net_br)
end
function add_Tbranches(net_bus, net_br, bus_index, BaseMVA, RawString, pre, start, size)
    # println("pre: ", pre);
    ind=0
    for i in 1:4:size
        Tran1=split(RawString[i+start-1], ",");
        Tran2=split(RawString[i+start], ",");
        Tran3=split(RawString[i+start+1], ",");
        Tran4=split(RawString[i+start+2], ",");
        St=parse(Int64, split(Tran1[12]," "; limit=0, keepempty=false)[1]);
        if (St==1)  ind+=1;
            net_br=merge(net_br, Dict(string(ind+pre)=>Dict{String,Any}()));
            I_index=parse(Int64, split(Tran1[1]," "; limit=0, keepempty=false)[1]);
            I=bus_index[string(I_index)];
            J_index=parse(Int64, split(Tran1[2]," "; limit=0, keepempty=false)[1]);
            J=bus_index[string(J_index)];
            CKT=split(Tran1[4],r"\W"; limit=0, keepempty=false);
            if (length(CKT)>0) CKT=split(CKT[1]," "; limit=0, keepempty=false)[1]; end
            MAG1=parse(Float64, split(Tran1[8]," "; limit=0, keepempty=false)[1]);
            MAG2=parse(Float64, split(Tran1[9]," "; limit=0, keepempty=false)[1]);

            R12=parse(Float64, split(Tran2[1]," "; limit=0, keepempty=false)[1]);
            X12=parse(Float64, split(Tran2[2]," "; limit=0, keepempty=false)[1]);
            WINDV1=parse(Float64, split(Tran3[1]," "; limit=0, keepempty=false)[1]);
            ANG1=parse(Float64, split(Tran3[3]," "; limit=0, keepempty=false)[1]);
            RATEA1=parse(Float64, split(Tran3[4]," "; limit=0, keepempty=false)[1])/BaseMVA;
            RATEC1=parse(Float64, split(Tran3[6]," "; limit=0, keepempty=false)[1])/BaseMVA;
            WINDV2=parse(Float64, split(Tran4[1]," "; limit=0, keepempty=false)[1]);

            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("br_r"=>R12));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("ckt"=>CKT));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("br_status"=>St));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("br_x"=>X12));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("b_fr"=>0.0));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("b_to"=>0.0));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("rate_a"=>RATEA1));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("rate_c"=>RATEC1));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("f_bus"=>I_index));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("t_bus"=>J_index));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("f_bus_index"=>I));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("t_bus_index"=>J));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("transformer"=>true));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("tap"=>WINDV1/WINDV2));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("shift"=>ANG1*pi/180));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("index"=>ind+pre)); #index=i
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("g_fr"=>0.0));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("g_to"=>0.0));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("gfM"=>MAG1));
            net_br[string(ind+pre)]=merge(net_br[string(ind+pre)], Dict("bfM"=>MAG2));

            gfM_total=net_bus[string(I)]["gfM"]+MAG1*St;
            bfM_total=net_bus[string(I)]["bfM"]+MAG2*St;
            net_bus[string(I)]=merge(net_bus[string(I)], Dict("gfM"=>gfM_total));
            net_bus[string(I)]=merge(net_bus[string(I)], Dict("bfM"=>bfM_total));
        end
    end
    return(net_bus,net_br)
end

function add_switch_shunts(net_bus, net_sshunt, bus_index, BaseMVA, RawString, start, size)
    net_sshunt=merge(net_sshunt, Dict("bus_index"=>zeros(Int32,0)));
    for i=1:size
        net_sshunt=merge(net_sshunt, Dict(string(i)=>Dict{String,Any}()));
        shunt=split(RawString[i+start-1], ",");
        I_index=parse(Int64, split(shunt[1]," "; limit=0, keepempty=false)[1]);
        I=bus_index[string(I_index)];
        STATUS=parse(Int64, split(shunt[4]," "; limit=0, keepempty=false)[1]);

        N=zeros(8); B=zeros(8); BL=zeros(8); Shmin=0; Shmax=0;
        for i=1:8
            N[i]=parse(Int64, split(shunt[2*(4+i)+1]," "; limit=0, keepempty=false)[1]);
            B[i]=parse(Float64, split(shunt[2*(5+i)]," "; limit=0, keepempty=false)[1]);
            BL[i]=N[i]*B[i]/BaseMVA; if (BL[i]>0) Shmax+=BL[i]*STATUS; elseif (BL[i]<0) Shmin+=BL[i]*STATUS; end
        end

        net_sshunt[string(i)]=merge(net_sshunt[string(i)], Dict("status"=>STATUS), Dict("shunt_bus"=>I_index), Dict("shunt_bus_index"=>I));
        net_sshunt[string(i)]=merge(net_sshunt[string(i)], Dict("shmax"=>Shmax), Dict("shmin"=>Shmin));

        Shmax+=net_bus[string(I)]["sh_max"]*STATUS;
        Shmin+=net_bus[string(I)]["sh_min"]*STATUS;
        net_bus[string(I)]=merge(net_bus[string(I)], Dict("sh_max"=>Shmax), Dict("sh_min"=>Shmin));
        push!(net_sshunt["bus_index"],I);
    end
    return(net_bus, net_sshunt);
end

function process_raw(PAROTS_flag, path_raw)
    network1=Dict{String,Any}();
    raw_string = readlines(open(path_raw));
    temp=split(raw_string[1], ","); baseMVA=parse(Float64, split(temp[2]," "; limit=0, keepempty=false)[1]);
    network1=merge(network1, Dict("baseMVA"=>baseMVA), Dict("bus"=>Dict{String,Any}()), Dict("sshunt"=>Dict{String,Any}()), Dict("gen"=>Dict{String,Any}()));
    network1=merge(network1, Dict("branch"=>Dict{String,Any}()), Dict("solution"=>Dict()));

    bus_start=4; bus_end=0; bus_check=true; bus_col=13;
    load_start=4; load_end=0; load_check=true; load_col=14;
    shunt_start=4; shunt_end=0; shunt_check=true; shunt_col=5;
    gen_start=4; gen_end=0; gen_check=true; gen_col=28;
    branch_start=4; branch_end=0; branch_check=true; branch_col=24;
    Sshunt_start=4; Sshunt_end=0; Sshunt_check=true; Sshunt_col=26;
    trans_start=4; trans_end=0; trans_check=true; Tcol1=21; Tcol2=3; Tcol3=17; Tcol4=2;
    for i=4:length(raw_string)
        temp=split(raw_string[i],","); temp_pre=split(raw_string[i-1],","); temp1=0; temp2=0; temp3=0; temp_pre_sh=0; temp_sh=0;
        if (length(temp)>=2) temp_sh=split(temp[2]," "; limit=0, keepempty=false)[1]; temp_sh=split(temp[2],""; limit=0, keepempty=false)[1]; end
        if (length(temp_pre)>=2) temp_pre_sh=split(temp_pre[2]," "; limit=0, keepempty=false)[1]; temp_pre_sh=split(temp_pre[2],""; limit=0, keepempty=false)[1]; end
        if (i<length(raw_string)-3) temp1=split(raw_string[i+1],","); temp2=split(raw_string[i+2],","); temp3=split(raw_string[i+3],","); end

        if (length(temp)==bus_col && bus_check) bus_start=i; bus_check=false; end
        if (length(temp_pre)==bus_col && length(temp)!=bus_col) bus_end=i-1; end

        if (length(temp)==load_col && load_check) load_start=i; load_check=false; end
        if (length(temp_pre)==load_col && length(temp)!=load_col) load_end=i-1; end

        if (length(temp)==shunt_col && temp_sh=="'" && shunt_check) shunt_start=i; shunt_check=false; end
        if (length(temp_pre)==shunt_col && (length(temp)!=shunt_col || temp_sh!="'") && temp_pre_sh=="'") shunt_end=i-1; end

        if (length(temp)==gen_col && gen_check) gen_start=i; gen_check=false; end
        if (length(temp_pre)==gen_col && length(temp)!=gen_col) gen_end=i-1; end

        if (length(temp)==branch_col && branch_check) branch_start=i; branch_check=false; end
        if (length(temp_pre)==branch_col && length(temp)!=branch_col) branch_end=i-1; end

        if (length(temp)==Tcol1 && length(temp1)==Tcol2 && length(temp2)==Tcol3 && length(temp3)==Tcol4 && trans_check) trans_start=i; trans_check=false; end
        if (length(temp)==Tcol4 && length(temp_pre)==Tcol3 && length(temp1)!=Tcol1) trans_end=i; end

        if (length(temp)==Sshunt_col && Sshunt_check) Sshunt_start=i; Sshunt_check=false; end
        if (length(temp_pre)==Sshunt_col && length(temp)!=Sshunt_col) Sshunt_end=i-1; end
    end

    (network1["bus"], bus_shift)=add_buses(network1["bus"], raw_string, bus_start, bus_end-bus_start+1);
    network1["bus"]=add_loads(network1["bus"], bus_shift, network1["baseMVA"], raw_string, load_start, load_end-load_start+1);
    network1["bus"]=add_shunts(network1["bus"], bus_shift, network1["baseMVA"], raw_string, shunt_start, shunt_end-shunt_start+1);
    (network1["bus"], network1["gen"], slack, gen_buses) =add_gens(network1["bus"], network1["gen"], bus_shift, network1["baseMVA"], raw_string, gen_start, gen_end-gen_start+1);
    network1=merge(network1, Dict("gen_buses"=>gen_buses, "slack"=>slack));
    network1["branch"]=add_Lbranches(network1["branch"], bus_shift, network1["baseMVA"], raw_string, branch_start, branch_end-branch_start+1);
    (network1["bus"], network1["branch"])=add_Tbranches(network1["bus"], network1["branch"], bus_shift, network1["baseMVA"], raw_string, length(network1["branch"]), trans_start, trans_end-trans_start+1);
    (network1["bus"], network1["sshunt"])=add_switch_shunts(network1["bus"], network1["sshunt"], bus_shift, network1["baseMVA"], raw_string, Sshunt_start, Sshunt_end-Sshunt_start+1);
    network1=merge(network1, Dict("nbal"=>Dict{String,Any}(), "source_type"=>"arpae"));
    for i=1:length(network1["bus"])
        network1["nbal"]=merge(network1["nbal"], Dict(string(i)=>Dict()));
        network1["nbal"][string(i)]=merge(network1["nbal"][string(i)], Dict("gen"=>zeros(Int32,0)), Dict("Sij"=>zeros(Int32,0)), Dict("Sji"=>zeros(Int32,0)), Dict("shunt"=>zeros(Int32,0)));
    end

    for g=1:length(network1["gen"]) gen=network1["gen"][string(g)];
        if (gen["gen_bus_index"]*gen["gen_status"]>0) push!(network1["nbal"][string(gen["gen_bus_index"])]["gen"],g); end
    end
    # println("nbr: ", length(network1["branch"]));
    for b=1:length(network1["branch"]) br=network1["branch"][string(b)];
        if (br["f_bus_index"]*br["br_status"]>0) push!(network1["nbal"][string(br["f_bus_index"])]["Sij"], b); end
        if (br["t_bus_index"]*br["br_status"]>0) push!(network1["nbal"][string(br["t_bus_index"])]["Sji"], b); end
    end
    for sh=1:length(network1["sshunt"])-1 shunt=network1["sshunt"][string(sh)];
        if (shunt["shunt_bus_index"]*shunt["status"]>0) push!(network1["nbal"][string(shunt["shunt_bus_index"])]["shunt"],sh); end
    end
    # if (PAROTS_flag>0) println("The RAW file was read successfully ...") end
    return(network1, bus_shift)
end

"--------------------------------------------------
   Readiang 'rop' file (Generator Cost Data File)
---------------------------------------------------"
function Generator_Dispatch_data(rop_string, gen_index, bus_index, GDD_in_start, size, generator)
    for i=1:size
        temp = split(rop_string[GDD_in_start+i-1],",");
        temp[2]=split(temp[2],r"\W"; limit=0, keepempty=false)[1];
        ID=split(temp[2]," "; limit=0, keepempty=false)[1];
        bus_no=parse(Int64,temp[1]);
        bus_no_index=bus_index[string(bus_no)];
        DSPTBL=parse(Int64,temp[4]);
        generator[string(gen_index[string(DSPTBL)])]=merge(generator[string(gen_index[string(DSPTBL)])], Dict("GenID"=>ID), Dict("bus"=>bus_no), Dict("bus_index"=>bus_no_index));
    end
    return(generator)
end

function process_rop(PAROTS_flag, path_rop, network, bus_index)
    rop_string = readlines(open(path_rop));
    generator=Dict{String,Any}();
    "Generator Dispatch Data"
    GDU_start=0; GDU_check=true; GDU_col=4; GDU_end=0;
    APDT_start=0; APDT_check=true; APDT_col=7; #APDT_end=0;
    PLCT_start=0; PLCT_check=true; PLCT_col=3; #PLCT_end=0;

    for i=1:length(rop_string)
        temp=split(rop_string[i],",");
        temp_pre= (i>1)  ? split(rop_string[i-1],",") : ["0 "," 0"] ;
	       temp_1=split(temp[1]," "; limit=0, keepempty=false)[1];
	          temp_1_pre=split(temp_pre[1]," "; limit=0, keepempty=false)[1];
        if (length(temp)==GDU_col && GDU_check && temp_1!="0") GDU_start=i; GDU_check=false; end
        if (length(temp_pre)==GDU_col && length(temp)!=GDU_col && temp_1_pre!="0") GDU_end=i-1; end
        if (length(temp)==APDT_col && APDT_check && temp_1!="0") APDT_start=i; APDT_check=false; end
        if (length(temp)==PLCT_col && PLCT_check && temp_1!="0") PLCT_start=i; PLCT_check=false; end
    end
    size=GDU_end-GDU_start+1;
    NPAIRSit=PLCT_start;
    NPAIRS=0;
    gen_index=Dict{String,Any}();
    for i=1:size
        temp = split(rop_string[NPAIRSit],","); NPAIRS = parse(Int64, temp[3]); nTable=parse(Int64, temp[1]);
        temp=split(rop_string[NPAIRSit+1],","); X=parse(Float64, temp[1]); Y=parse(Float64, temp[2]);
        generator=merge(generator, Dict(string(nTable)=>Dict("LinPoints"=>NPAIRS, "LinCost"=>Dict())));
        baseMVA=network["baseMVA"];
        generator[string(nTable)]["LinCost"]=merge(generator[string(nTable)]["LinCost"], Dict("1"=>[X/baseMVA,Y/baseMVA]));
        for c=2:NPAIRS
            temp=split(rop_string[NPAIRSit+c],","); X=parse(Float64, temp[1]); Y=parse(Float64, temp[2]);
            generator[string(nTable)]["LinCost"]=merge(generator[string(nTable)]["LinCost"], Dict(string(c)=>[X/baseMVA,Y/baseMVA]));
        end
        NPAIRSit+=NPAIRS+1;
    end
    for i=1:size
        temp1 = split(rop_string[APDT_start+i-1],",");
        index = parse(Int64, temp1[7]); TBL = parse(Int64, temp1[1]); gen_index=merge(gen_index, Dict(string(TBL)=>index));
        generator[string(index)]=merge(generator[string(index)], Dict("gov_droop"=>0.0));
    end
    for i=1:size
        temp2 = split(rop_string[GDU_start+i-1],",");
        temp2[2]=split(temp2[2],r"\W"; limit=0, keepempty=false)[1];
        ID=split(temp2[2]," "; limit=0, keepempty=false)[1];
        bus_no=parse(Int64,temp2[1]);
        bus_no_index=bus_index[string(bus_no)];
        DSPTBL=parse(Int64,temp2[4]);
        generator[string(gen_index[string(DSPTBL)])]=merge(generator[string(gen_index[string(DSPTBL)])], Dict("GenID"=>ID), Dict("bus"=>bus_no), Dict("bus_index"=>bus_no_index));
    end
    for i=1:length(generator) if (!haskey(generator[string(i)], "LinPoints") && (PAROTS_flag)) println("generator missing info at ", i) end end
    if (size==length(network["gen"])) Gen_Dict=deepcopy(network["gen"]);
        for i=1:size
            for j in Gen_Dict
                bg=j[2]["gen_bus"]; Ug=j[2]["id"];
                if (generator[string(i)]["bus"]==bg && generator[string(i)]["GenID"]==Ug)
                    LC=generator[string(i)]["LinCost"]; LP=generator[string(i)]["LinPoints"];
                    network["gen"][string(j[1])]=merge(network["gen"][string(j[1])], Dict("LinCost"=>LC), Dict("LinPoints"=>LP));
                    delete!(Gen_Dict, j[1]);
                end
            end
        end
        # if (PAROTS_flag>0) println("The ROP file was read successfully ...") end
    end
    # elseif (PAROTS_flag) println("ERROR | The number of generators in RAW and ROP files don't match") end
    return(network)
end

"--------------------------------------------------
Readiang 'inl' file (Governor Response Data File)
---------------------------------------------------"
function process_inl(PAROTS_flag, path_inl, network)
    inl_string = readlines(open(path_inl));
    Gen_Dict=deepcopy(network["gen"]);
    for i=1:length(inl_string)
        temp=split(inl_string[i],",");
        if (length(temp)==7)
            bus=parse(Int64, split(temp[1]," "; limit=0, keepempty=false)[1]);
            ID_temp=split(temp[2],r"\W"; limit=0, keepempty=false)[1]; genID=split(ID_temp," "; limit=0, keepempty=false)[1];
            for j in Gen_Dict
                bg=j[2]["gen_bus"]; Ug=j[2]["id"];
                if (bus==bg && genID==Ug)
                    gov_droop=parse(Float64, temp[6]);
                    network["gen"][string(j[1])]=merge(network["gen"][string(j[1])], Dict("gov_droop"=>gov_droop));
                    delete!(Gen_Dict, j[1]);
                end
            end
        end
    end
    sum_apf=0
    for i=1:length(network["gen"]) gen=network["gen"][string(i)]; st=gen["gen_status"]
        sum_apf+=st*gen["gov_droop"];
    end
    for i=1:length(network["gen"]) gen=network["gen"][string(i)]; st=gen["gen_status"]
        total_apf=sum_apf-st*gen["gov_droop"];
        network["gen"][string(i)]=merge(network["gen"][string(i)], Dict("total_apf"=>total_apf));
    end
    # if (PAROTS_flag>0) println("The INL file was read successfully ..."); end
    return(network);
end

"-------------------------------------------------------
Readiang 'con' file (Contingency Description Data File)
--------------------------------------------------------"
function adding_con_line_info(contingency, counts, Branch_Dict)
    fc=contingency[string(counts)]["f_bus"];
    tc=contingency[string(counts)]["t_bus"];
    CTc=contingency[string(counts)]["CT"];
    for j in Branch_Dict
        fb=j[2]["f_bus"];
        tb=j[2]["t_bus"];
        CTb=j[2]["ckt"];
        if (fb==fc && tb==tc && CTb==CTc)
            index=j[2]["index"];
            contingency[string(counts)]=merge(contingency[string(counts)], Dict("branch"=>index));
        end
    end
    return(contingency)
end

function adding_con_gen_info(contingency, counts, contingency_G, Gcounts, Gen_Dict)
    bc=contingency[string(counts)]["bus"]; Uc=contingency[string(counts)]["unit"];
    for j in Gen_Dict
        bg=j[2]["gen_bus"]; Ug=j[2]["id"];
        if (bc==bg && Uc==Ug)
            contingency[string(counts)]=merge(contingency[string(counts)], Dict("GenNum"=>j[2]["index"]));
            contingency_G[string(Gcounts)]=merge(contingency_G[string(Gcounts)], Dict("GenNum"=>j[2]["index"]));
        end
    end
    return(contingency, contingency_G);
end

function process_con(PAROTS_flag, path_con, network, bus_index)
    con_string = readlines(open(path_con));
    contingency=Dict{String,Any}(); contingency_G=Dict{String,Any}();
    counts=0; Gcounts=0
    for i=1:length(con_string)-1
        temp1 = split(con_string[i]," "; limit=0, keepempty=false); temp2 = split(con_string[i+1]," "; limit=0, keepempty=false);

        if (length(temp1)==2 && length(temp2)<6)
            label=temp1[2]; counts+=1;
            contingency[string(counts)]=Dict("status"=>false);
            contingency[string(counts)]=merge(contingency[string(counts)], Dict("label"=>label), Dict("type"=>NaN));
        elseif (length(temp1)==2 && length(temp2)==10)
            label=temp1[2]; type = split(temp1[2],""; limit=0, keepempty=false)[1]; counts+=1; f_bus_index=parse(Int64,temp2[5]); f_bus=bus_index[string(f_bus_index)];
            contingency[string(counts)]=Dict("type"=>type, "status"=>true, "label"=>label, "f_bus"=>f_bus_index, "f_bus_index"=>f_bus);
            t_bus_index=parse(Int64,temp2[8]); t_bus=bus_index[string(t_bus_index)]; CT_bus=temp2[10];
            contingency[string(counts)]=merge(contingency[string(counts)], Dict("t_bus"=>t_bus_index, "t_bus_index"=>t_bus, "CT"=>CT_bus));
            contingency=adding_con_line_info(contingency, counts, network["branch"]);
        elseif (length(temp1)==2 && length(temp2)==6)
            label=temp1[2];
            counts+=1;
            Gcounts+=1;
            unit=temp2[3];
            f_bus_index=parse(Int64,temp2[6]);
            f_bus=bus_index[string(f_bus_index)];
            contingency[string(counts)]=Dict("type"=>"G", "status"=>true, "label"=>label, "unit"=>unit, "bus"=>f_bus_index, "bus_index"=>f_bus);
            contingency_G[string(Gcounts)]=Dict("type"=>"G", "label"=>label, "unit"=>unit, "bus"=>f_bus_index, "bus_index"=>f_bus);
            (contingency, contingency_G)=adding_con_gen_info(contingency, counts, contingency_G, Gcounts, network["gen"]);
        end
    end

    # if (PAROTS_flag>0) println("The CON file was read successfully ..."); end
    return(contingency);
end

"--------------------------------------------------
   Readiang 'raw' file for LODF calculation
---------------------------------------------------"
function process_net_LODF(PAROTS_flag, path_raw)
    network1=Dict{String,Any}();
    raw_string = readlines(open(path_raw));
    temp=split(raw_string[1], ","); BaseMVA=parse(Float64, split(temp[2]," "; limit=0, keepempty=false)[1]);
    network1=merge(network1, Dict("baseMVA"=>BaseMVA), Dict("bus"=>Dict{String,Any}()), Dict("sshunt"=>Dict{String,Any}()), Dict("gen"=>Dict{String,Any}()));
    network1=merge(network1, Dict("branch"=>Dict{String,Any}()), Dict("solution"=>Dict()));

    bus_start=4; bus_end=0; bus_check=true; bus_col=13;
    gen_start=4; gen_end=0; gen_check=true; gen_col=28;
    branch_start=4; branch_end=0; branch_check=true; branch_col=24;
    trans_start=4; trans_end=0; trans_check=true; Tcol1=21; Tcol2=3; Tcol3=17; Tcol4=2;
    for i=4:length(raw_string)
        temp=split(raw_string[i],","); temp_pre=split(raw_string[i-1],","); temp1=0; temp2=0; temp3=0; temp_pre_sh=0; temp_sh=0;
        if (length(temp)>=2) temp_sh=split(temp[2]," "; limit=0, keepempty=false)[1]; temp_sh=split(temp[2],""; limit=0, keepempty=false)[1]; end
        if (length(temp_pre)>=2) temp_pre_sh=split(temp_pre[2]," "; limit=0, keepempty=false)[1]; temp_pre_sh=split(temp_pre[2],""; limit=0, keepempty=false)[1]; end
        if (i<length(raw_string)-3) temp1=split(raw_string[i+1],","); temp2=split(raw_string[i+2],","); temp3=split(raw_string[i+3],","); end

        if (length(temp)==bus_col && bus_check) bus_start=i; bus_check=false; end
        if (length(temp_pre)==bus_col && length(temp)!=bus_col) bus_end=i-1; end

        if (length(temp)==gen_col && gen_check) gen_start=i; gen_check=false; end
        if (length(temp_pre)==gen_col && length(temp)!=gen_col) gen_end=i-1; end

        if (length(temp)==branch_col && branch_check) branch_start=i; branch_check=false; end
        if (length(temp_pre)==branch_col && length(temp)!=branch_col) branch_end=i-1; end

        if (length(temp)==Tcol1 && length(temp1)==Tcol2 && length(temp2)==Tcol3 && length(temp3)==Tcol4 && trans_check) trans_start=i; trans_check=false; end
        if (length(temp)==Tcol4 && length(temp_pre)==Tcol3 && length(temp1)!=Tcol1) trans_end=i; end
    end
    bus_size=bus_end-bus_start+1; bus_shift=Dict{String,Any}();
    for i=1:bus_size
        bus=split(raw_string[i+bus_start-1], ",");
        I=parse(Int64, split(bus[1]," "; limit=0, keepempty=false)[1]);
        bus_shift=merge(bus_shift, Dict(string(I)=>i));
    end
    gen_size=gen_end-gen_start+1;
    temp_pamx_id=1; max_gen=0.0;
    for i=1:gen_size
        gen=split(raw_string[i+gen_start-1], ",");
        I_index=parse(Int64, split(gen[1]," "; limit=0, keepempty=false)[1]);
        I=bus_shift[string(I_index)];
        St=parse(Int64, split(gen[15]," "; limit=0, keepempty=false)[1]);
        PT=parse(Float64, split(gen[17]," "; limit=0, keepempty=false)[1])/BaseMVA;
        if (max_gen<St*PT) max_gen=St*PT; temp_pamx_id=i; end
    end
    network1=merge(network1, Dict("slack"=>temp_pamx_id));
    network1=merge(network1, Dict("nbus"=>bus_size));
    branch_size=branch_end-branch_start+1; i=0;
    for ind=1:branch_size
        Lbr=split(raw_string[ind+branch_start-1], ",");
        St=parse(Int64, split(Lbr[14]," "; limit=0, keepempty=false)[1]);
        if (St==1) i+=1;
            network1["branch"]=merge(network1["branch"], Dict(string(i)=>Dict{String,Any}()));
            I_index=parse(Int64, split(Lbr[1]," "; limit=0, keepempty=false)[1]);
            I=bus_shift[string(I_index)];
            J_index=parse(Int64, split(Lbr[2]," "; limit=0, keepempty=false)[1]);
            J=bus_shift[string(J_index)];
            CKT=split(Lbr[3],r"\W"; limit=0, keepempty=false);
            if (length(CKT)>0) CKT=split(CKT[1]," "; limit=0, keepempty=false)[1]; end

            R=parse(Float64, split(Lbr[4]," "; limit=0, keepempty=false)[1]);
            X=parse(Float64, split(Lbr[5]," "; limit=0, keepempty=false)[1]);
            B=parse(Float64, split(Lbr[6]," "; limit=0, keepempty=false)[1]);
            RATEC=parse(Float64, split(Lbr[9]," "; limit=0, keepempty=false)[1])/BaseMVA;

            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("br_r"=>R));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("br_status"=>St, "index"=>i));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("br_x"=>X,"ratec"=>RATEC));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("b_fr"=>B/2));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("b_to"=>B/2));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("f_bus"=>I_index, "t_bus"=>J_index, "ckt"=>CKT));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("f_bus_index"=>I));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("t_bus_index"=>J));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("tap"=>1.0));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("shift"=>0.0));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("gfM"=>0.0));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("bfM"=>0.0));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("g_fr"=>0.0));
            network1["branch"][string(i)]=merge(network1["branch"][string(i)], Dict("g_to"=>0.0));
        end
    end
    branch_size=length(network1["branch"]);
    tbranch_size=trans_end-trans_start+1; ind=0;
    for i in 1:4:tbranch_size
        Tran1=split(raw_string[i+trans_start-1], ",");
        Tran2=split(raw_string[i+trans_start], ",");
        Tran3=split(raw_string[i+trans_start+1], ",");
        Tran4=split(raw_string[i+trans_start+2], ",");
        St=parse(Int64, split(Tran1[12]," "; limit=0, keepempty=false)[1]);
        if (St==1) ind+=1;
            network1["branch"]=merge(network1["branch"], Dict(string(ind+branch_size)=>Dict{String,Any}()));
            I_index=parse(Int64, split(Tran1[1]," "; limit=0, keepempty=false)[1]);
            I=bus_shift[string(I_index)];
            J_index=parse(Int64, split(Tran1[2]," "; limit=0, keepempty=false)[1]);
            J=bus_shift[string(J_index)];
            CKT=split(Tran1[4],r"\W"; limit=0, keepempty=false);
            if (length(CKT)>0) CKT=split(CKT[1]," "; limit=0, keepempty=false)[1]; end
            MAG1=parse(Float64, split(Tran1[8]," "; limit=0, keepempty=false)[1]);
            MAG2=parse(Float64, split(Tran1[9]," "; limit=0, keepempty=false)[1]);
            R12=parse(Float64, split(Tran2[1]," "; limit=0, keepempty=false)[1]);
            X12=parse(Float64, split(Tran2[2]," "; limit=0, keepempty=false)[1]);
            WINDV1=parse(Float64, split(Tran3[1]," "; limit=0, keepempty=false)[1]);
            ANG1=parse(Float64, split(Tran3[3]," "; limit=0, keepempty=false)[1]);
            WINDV2=parse(Float64, split(Tran4[1]," "; limit=0, keepempty=false)[1]);
            RATEC1=parse(Float64, split(Tran3[6]," "; limit=0, keepempty=false)[1])/BaseMVA;

            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("br_r"=>R12));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("br_status"=>St, "index"=>ind+branch_size));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("br_x"=>X12));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("b_fr"=>0.0,"ratec"=>RATEC1));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("b_to"=>0.0));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("f_bus"=>I_index, "t_bus"=>J_index, "ckt"=>CKT));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("f_bus_index"=>I));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("t_bus_index"=>J));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("transformer"=>true));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("tap"=>WINDV1/WINDV2));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("shift"=>deg2rad(ANG1)));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("g_fr"=>0.0));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("g_to"=>0.0));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("gfM"=>MAG1));
            network1["branch"][string(ind+branch_size)]=merge(network1["branch"][string(ind+branch_size)], Dict("bfM"=>MAG2));
        end
    end
    return(network1, bus_shift)
end
"--------------------------------------------------
        find slack bus from 'raw' file
---------------------------------------------------"
function find_slack(path_raw)
    raw_string = readlines(open(path_raw));
    temp=split(raw_string[1], ","); BaseMVA=parse(Float64, split(temp[2]," "; limit=0, keepempty=false)[1]);
    bus_start=4; bus_end=0; bus_check=true; bus_col=13;
    gen_start=4; gen_end=0; gen_check=true; gen_col=28;
    for i=4:length(raw_string)
        temp=split(raw_string[i],","); temp_pre=split(raw_string[i-1],",");
        if (length(temp)==bus_col && bus_check) bus_start=i; bus_check=false; end
        if (length(temp_pre)==bus_col && length(temp)!=bus_col) bus_end=i-1; end
        if (length(temp)==gen_col && gen_check) gen_start=i; gen_check=false; end
        if (length(temp_pre)==gen_col && length(temp)!=gen_col) gen_end=i-1; end
    end
    bus_size=bus_end-bus_start+1; bus_shift=Dict{String,Any}();
    for i=1:bus_size
        bus=split(raw_string[i+bus_start-1], ",");
        I=parse(Int64, split(bus[1]," "; limit=0, keepempty=false)[1]);
        bus_shift=merge(bus_shift, Dict(string(I)=>i));
    end
    gen_size=gen_end-gen_start+1; index=0;
    for i=1:gen_size
        gen=split(raw_string[i+gen_start-1], ",");
        I_index=parse(Int64, split(gen[1]," "; limit=0, keepempty=false)[1]);
        I=bus_shift[string(I_index)];
        St=parse(Int64, split(gen[15]," "; limit=0, keepempty=false)[1]);
        if (St==1) index=I; break; end
    end
    return(index)
end

"--------------------------------------------------
 Readiang 'raw' file for node4 (Network Data File)
---------------------------------------------------"
function process_raw_node4(PAROTS_flag, path_raw)
    network1=Dict{String,Any}();
    raw_string = readlines(open(path_raw));
    temp=split(raw_string[1], ","); baseMVA=parse(Float64, split(temp[2]," "; limit=0, keepempty=false)[1]);
    network1=merge(network1, Dict("baseMVA"=>baseMVA), Dict("bus"=>Dict{String,Any}()), Dict("sshunt"=>Dict{String,Any}()), Dict("gen"=>Dict{String,Any}()));
    network1=merge(network1, Dict("branch"=>Dict{String,Any}()), Dict("solution"=>Dict()));

    bus_start=4; bus_end=0; bus_check=true; bus_col=13;
    load_start=4; load_end=0; load_check=true; load_col=14;
    shunt_start=4; shunt_end=0; shunt_check=true; shunt_col=5;
    gen_start=4; gen_end=0; gen_check=true; gen_col=28;
    branch_start=4; branch_end=0; branch_check=true; branch_col=24;
    Sshunt_start=4; Sshunt_end=0; Sshunt_check=true; Sshunt_col=26;
    trans_start=4; trans_end=0; trans_check=true; Tcol1=21; Tcol2=3; Tcol3=17; Tcol4=2;
    for i=4:length(raw_string)
        temp=split(raw_string[i],","); temp_pre=split(raw_string[i-1],","); temp1=0; temp2=0; temp3=0; temp_pre_sh=0; temp_sh=0;
        if (length(temp)>=2) temp_sh=split(temp[2]," "; limit=0, keepempty=false)[1]; temp_sh=split(temp[2],""; limit=0, keepempty=false)[1]; end
        if (length(temp_pre)>=2) temp_pre_sh=split(temp_pre[2]," "; limit=0, keepempty=false)[1]; temp_pre_sh=split(temp_pre[2],""; limit=0, keepempty=false)[1]; end
        if (i<length(raw_string)-3) temp1=split(raw_string[i+1],","); temp2=split(raw_string[i+2],","); temp3=split(raw_string[i+3],","); end

        if (length(temp)==bus_col && bus_check) bus_start=i; bus_check=false; end
        if (length(temp_pre)==bus_col && length(temp)!=bus_col) bus_end=i-1; end

        if (length(temp)==load_col && load_check) load_start=i; load_check=false; end
        if (length(temp_pre)==load_col && length(temp)!=load_col) load_end=i-1; end

        if (length(temp)==shunt_col && temp_sh=="'" && shunt_check) shunt_start=i; shunt_check=false; end
        if (length(temp_pre)==shunt_col && (length(temp)!=shunt_col || temp_sh!="'") && temp_pre_sh=="'") shunt_end=i-1; end

        if (length(temp)==gen_col && gen_check) gen_start=i; gen_check=false; end
        if (length(temp_pre)==gen_col && length(temp)!=gen_col) gen_end=i-1; end

        if (length(temp)==branch_col && branch_check) branch_start=i; branch_check=false; end
        if (length(temp_pre)==branch_col && length(temp)!=branch_col) branch_end=i-1; end

        if (length(temp)==Tcol1 && length(temp1)==Tcol2 && length(temp2)==Tcol3 && length(temp3)==Tcol4 && trans_check) trans_start=i; trans_check=false; end
        if (length(temp)==Tcol4 && length(temp_pre)==Tcol3 && length(temp1)!=Tcol1) trans_end=i; end

        if (length(temp)==Sshunt_col && Sshunt_check) Sshunt_start=i; Sshunt_check=false; end
        if (length(temp_pre)==Sshunt_col && length(temp)!=Sshunt_col) Sshunt_end=i-1; end
    end

    (network1["bus"], bus_shift)=add_buses(network1["bus"], raw_string, bus_start, bus_end-bus_start+1);
    (network1["bus"], network1["gen"], slack) =add_gens(network1["bus"], network1["gen"], bus_shift, network1["baseMVA"], raw_string, gen_start, gen_end-gen_start+1);
    # @save "slack.jld2" slack
    network1["bus"]=add_loads(network1["bus"], bus_shift, network1["baseMVA"], raw_string, load_start, load_end-load_start+1);
    network1["bus"]=add_shunts(network1["bus"], bus_shift, network1["baseMVA"], raw_string, shunt_start, shunt_end-shunt_start+1);
    network1=merge(network1, Dict("slack"=>slack));
    network1["branch"]=add_Lbranches(network1["branch"], bus_shift, network1["baseMVA"], raw_string, branch_start, branch_end-branch_start+1);
    (network1["bus"], network1["branch"])=add_Tbranches(network1["bus"], network1["branch"], bus_shift, network1["baseMVA"], raw_string, branch_end-branch_start+1, trans_start, trans_end-trans_start+1);
    (network1["bus"], network1["sshunt"])=add_switch_shunts(network1["bus"], network1["sshunt"], bus_shift, network1["baseMVA"], raw_string, Sshunt_start, Sshunt_end-Sshunt_start+1);
    network1=merge(network1, Dict("nbal"=>Dict{String,Any}()));
    for i=1:length(network1["bus"])
        network1["nbal"]=merge(network1["nbal"], Dict(string(i)=>Dict()));
        network1["nbal"][string(i)]=merge(network1["nbal"][string(i)], Dict("gen"=>zeros(Int32,0)), Dict("Sij"=>zeros(Int32,0)), Dict("Sji"=>zeros(Int32,0)), Dict("shunt"=>zeros(Int32,0)));
    end

    for g=1:length(network1["gen"]) gen=network1["gen"][string(g)];
        if (gen["gen_bus_index"]*gen["gen_status"]>0) push!(network1["nbal"][string(gen["gen_bus_index"])]["gen"],g); end
    end
    for b=1:length(network1["branch"]) br=network1["branch"][string(b)];
        if (br["f_bus_index"]*br["br_status"]>0) push!(network1["nbal"][string(br["f_bus_index"])]["Sij"], b); end
        if (br["t_bus_index"]*br["br_status"]>0) push!(network1["nbal"][string(br["t_bus_index"])]["Sji"], b); end
    end
    for sh=1:length(network1["sshunt"])-1 shunt=network1["sshunt"][string(sh)];
        if (shunt["shunt_bus_index"]*shunt["status"]>0) push!(network1["nbal"][string(shunt["shunt_bus_index"])]["shunt"],sh); end
    end
    return(network1, bus_shift)
end

"--------------------------------------------------
 Readiang 'con' file for node5 (Network Data File)
---------------------------------------------------"
function adding_con_line_info_node5(contingency, Branch_Dict, LconList)
    fc=contingency["f_bus"];
    tc=contingency["t_bus"];
    CTc=contingency["CT"];
    for j in Branch_Dict
        fb=j[2]["f_bus"];
        tb=j[2]["t_bus"];
        CTb=j[2]["ckt"];
        if (fb==fc && tb==tc && CTb==CTc)
            index=j[2]["index"];
            append!(LconList,index);
        end
    end
    return(LconList)
end
function process_con_node5(PAROTS_flag, path_con, branch, bus_index)
    con_string = readlines(open(path_con));
    LconList=zeros(Int32,0);
    for i=1:length(con_string)-1
        temp1 = split(con_string[i]," "; limit=0, keepempty=false); temp2 = split(con_string[i+1]," "; limit=0, keepempty=false);
        if (length(temp1)==2 && length(temp2)==10)
            f_bus_index=parse(Int64,temp2[5]); f_bus=bus_index[string(f_bus_index)];
            t_bus_index=parse(Int64,temp2[8]); t_bus=bus_index[string(t_bus_index)]; CT_bus=temp2[10];
            contingency=Dict("f_bus"=>f_bus_index, "t_bus"=>t_bus_index, "CT"=>CT_bus);
            LconList=adding_con_line_info_node5(contingency, branch, LconList);
        end
    end
    # data=Dict("data"=>LconList);
    return(LconList)
    # @save "LconList.bson" data
    # @save "LconList.bson" Dict("data"=>LconList)
end


text:
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

function reachable(adjacency::Matrix, src::Int, dest::Int)
    n = size(adjacency, 1)
    reachable = Vector(undef, n)
    reachable .= false
    reachable[src] = true
    for i=1:n
        for node=1:n
            reachable[node] |= reduce(|, reachable .& adjacency[:, node])
        end
    end
    reachable[dest]
end

n = 5

# run on sampled deterministic graph
adjacency_sampled = rand(Bool, n, n)
r = reachable(adjacency_sampled, 1, n)
println("Sampled graph reachability: ", r)

# run on random graph
code = @dice begin
    adjacency_rand = [flip(0.5) for i=1:n, j=1:n]
    reachable(adjacency_rand, 1, n)
end

# BDD analysis
bdd = compile(code)
num_flips(bdd), num_nodes(bdd)
@assert num_flips(bdd) == n*n-3n+3
infer(code, :bdd)

# IR analysis
to_dice_ir(code)

# # TODO fix loopy let identifiers
# has_dice_binary() && rundice(code)
# has_dice_binary() && infer(code, :ocaml)

ext: jl


text:
__precompile__()

module ForwardDiff

using DiffRules, DiffResults
using DiffResults: DiffResult, MutableDiffResult, ImmutableDiffResult
using StaticArrays
using Compat

import NaNMath
import SpecialFunctions
import CommonSubexpressions

include("prelude.jl")
include("partials.jl")
include("dual.jl")
include("config.jl")
include("apiutils.jl")
include("derivative.jl")
include("gradient.jl")
include("jacobian.jl")
include("hessian.jl")
include("deprecated.jl")

# This is a deprecation binding and should be removed in the next minor release.
const DiffBase = DiffResults

export DiffBase, DiffResults

end # module


text:
module AtomicGraphNets

export AGNConv, AGNPool#, AGNConvDEQ
include("layers.jl")
using .layers: AGNConv, AGNPool

include("models.jl")
export Xie_model

end


text:
using Test
using DataConvenience

x = Vector{Union{Missing, Int64}}(undef, 1_000_000)

@testset "count_missing" begin
    @test count(ismissing, x) == count_missing(x)
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
@testset "GEMM" begin
    # using LoopVectorization, LinearAlgebra, Test; T = Float64
    Unum, Tnum = LoopVectorization.register_count() == 16 ? (2, 6) : (3, 9)
    Unumt, Tnumt = LoopVectorization.register_count() == 16 ? (2, 6) : (5, 5)
    if LoopVectorization.register_count() != 8
        @test @inferred(LoopVectorization.matmul_params()) == (Unum, Tnum)
    end
    AmulBtq1 = :(for m ∈ axes(A,1), n ∈ axes(B,2)
                 C[m,n] = zeroB
                 for k ∈ axes(A,2)
                 C[m,n] += A[m,k] * B[n,k]
                 end
                 end);
    lsAmulBt1 = LoopVectorization.loopset(AmulBtq1);
    if LoopVectorization.register_count() != 8
        @test LoopVectorization.choose_order(lsAmulBt1) == (Symbol[:n,:m,:k], :m, :n, :m, Unum, Tnum)
    end
    AmulBq1 = :(for n ∈ axes(B,2), m ∈ axes(A,1)
                C[m,n] = 0.0
                for k ∈ axes(A,2)
                C[m,n] += A[m,k] * B[k,n]
                end
                end);
    lsAmulB1 = LoopVectorization.loopset(AmulBq1);
    if LoopVectorization.register_count() != 8
        @test LoopVectorization.choose_order(lsAmulB1) == (Symbol[:n,:m,:k], :m, :n, :m, Unum, Tnum)
    end
    AmulBq2 = :(for m ∈ 1:M, n ∈ 1:N
                C[m,n] = zero(eltype(B))
                for k ∈ 1:K
                C[m,n] += A[m,k] * B[k,n]
                end
                end)
    lsAmulB2 = LoopVectorization.loopset(AmulBq2);
    if LoopVectorization.register_count() != 8
        @test LoopVectorization.choose_order(lsAmulB2) == (Symbol[:n,:m,:k], :m, :n, :m, Unum, Tnum)
    end
    AmulBq3 = :(for m ∈ axes(A,1), n ∈ axes(B,2)
                ΔCₘₙ = zero(eltype(C))
                for k ∈ axes(A,2)
                ΔCₘₙ += A[m,k] * B[k,n]
                end
                C[m,n] += ΔCₘₙ
                end)
    lsAmulB3 = LoopVectorization.loopset(AmulBq3);
    if LoopVectorization.register_count() != 8
        @test LoopVectorization.choose_order(lsAmulB3) == (Symbol[:n,:m,:k], :m, :n, :m, Unum, Tnum)
    end
    if LoopVectorization.register_count() != 8
        for (fA,fB,v,Un,Tn) ∈ [(identity,identity,:m,Unum,Tnum),(adjoint,identity,:k,Unumt,Tnumt),(identity,adjoint,:m,Unum,Tnum),(adjoint,adjoint,:n,Unum,Tnum)]
            A = fA(rand(2,2))
            B = fB(rand(2,2))
            C = similar(A)
            ls = LoopVectorization.@avx_debug for m ∈ axes(A,1), n ∈ axes(B,2)
                ΔCₘₙ = zero(eltype(C))
                for k ∈ axes(A,2)
                    ΔCₘₙ += A[m,k] * B[k,n]
                end
                C[m,n] += ΔCₘₙ
            end
            (m, n) = v === :n ? (:n, :m) : (:m, :n)
            @test LoopVectorization.choose_order(ls) == (Symbol[:n,:m,:k], m, n, v, Un, Tn)
        end
    end
    function AmulB!(C, A, B)
        C .= 0
        for k ∈ axes(A,2), j ∈ axes(B,2)
            @simd ivdep for i ∈ axes(A,1)
                @inbounds C[i,j] += A[i,k] * B[k,j]
            end
        end
    end
    function AmulBavx1!(C, A, B)
        @avx for m ∈ 1:size(A,1), n ∈ axes(B,2)
            Cₘₙ = zero(eltype(C))
            for k ∈ 1:size(A,2)
                Cₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end
    function AmulBavx2!(C, A, B)
        z = zero(eltype(C))
        @avx unroll=(2,1) for m ∈ axes(A,1), n ∈ axes(B,2)
            C[m,n] = z
            for k ∈ axes(A,2)
                C[m,n] += A[m,k] * B[k,n]
            end
        end
    end
    function AmulBavx3!(C, A, B)
        @avx unroll=(2,2) for m ∈ axes(A,1), n ∈ axes(B,2)
            C[m,n] = zero(eltype(C))
            for k ∈ axes(A,2)
                C[m,n] += A[m,k] * B[k,n]
            end
        end
    end
    myzero(A) = zero(eltype(A))
    # function AmulBavx4!(C, A, B)
    #     @avx for m ∈ axes(A,1), n ∈ axes(B,2)
    #         C[m,n] = myzero(C)
    #         for k ∈ axes(A,2)
    #             C[m,n] += A[m,k] * B[k,n]
    #         end
    #     end
    # end
    # C = Cs; A = Ats'; B = Bs; factor = 1;
    # ls = LoopVectorization.@avx_debug for m ∈ axes(A,1), n ∈ axes(B,2)
    #         ΔCₘₙ = zero(eltype(C))
    #         for k ∈ axes(A,2)
    #             ΔCₘₙ += A[m,k] * B[k,n]
    #         end
    #         C[m,n] += ΔCₘₙ * factor
    #     end;
    function AmuladdBavx!(C, A, B, α = one(eltype(C)))
        @avx unroll=(2,2) for m ∈ indices((A,C),1), n ∈ indices((B,C),2)
            ΔCₘₙ = zero(eltype(C))
            for k ∈ indices((A,B),(2,1))
                ΔCₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] += α * ΔCₘₙ
        end
    end
    function AmuladdBavx!(C, A, B, α, β)# = zero(eltype(C)))
        @avx unroll=(1,1) for m ∈ axes(A,1), n ∈ axes(B,2)
            ΔCₘₙ = zero(eltype(C))
            for k ∈ axes(A,2)
                ΔCₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = α * ΔCₘₙ + β * C[m,n]
        end
    end
    Amuladdq = :(for m ∈ axes(A,1), n ∈ axes(B,2)
                 ΔCₘₙ = zero(eltype(C))
                 for k ∈ axes(A,2)
                 ΔCₘₙ += A[m,k] * B[k,n]
                 end
                 C[m,n] = α * ΔCₘₙ + β * C[m,n]
                 end);
    lsAmuladd = LoopVectorization.loopset(Amuladdq);
    if LoopVectorization.register_count() != 8
        @test LoopVectorization.choose_order(lsAmuladd) == (Symbol[:n,:m,:k], :m, :n, :m, Unum, Tnum)
    end
    Atmuladdq = :(for m ∈ axes(A,2), n ∈ axes(B,2)
                 ΔCₘₙ = zero(eltype(C))
                 for k ∈ axes(A,1)
                 ΔCₘₙ += A[k,m] * B[k,n]
                  end
                  C[m,n] += α * ΔCₘₙ
                 end);
    lsAtmuladd = LoopVectorization.loopset(Atmuladdq);
    # LoopVectorization.lower(lsAtmuladd, 2, 2)
    # lsAmuladd.operations
    # LoopVectorization.loopdependencies.(lsAmuladd.operations)
    # LoopVectorization.reduceddependencies.(lsAmuladd.operations)
    if LoopVectorization.register_count() != 8
        @test LoopVectorization.choose_order(lsAtmuladd) == (Symbol[:n,:m,:k], :m, :n, :k, Unumt, Tnumt)
    end

    function AmulB_avx1!(C, A, B)
        @_avx unroll=(2,2) for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
            Cₘₙ = zero(eltype(C))
            for k ∈ axes(A,2)
                Cₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end
    fq = :(for m ∈ axes(A,1), n ∈ axes(B,2)
           Cₘₙ = zero(eltype(C))
           for k ∈ axes(A,2)
           Cₘₙ += A[m,k] * B[k,n]
           end
           C[m,n] = Cₘₙ
           end);
    # exit()
    #         using LoopVectorization, Test
    #         T = Float64
    #         M = 77
    #         A = rand(M, M); B = rand(M, M); C = similar(A);
    function AmulB_avx2!(C, A, B)
        z = zero(eltype(C))
        @_avx unroll=(2,2) for m ∈ axes(A,1), n ∈ axes(B,2)
            C[m,n] = z
            for k ∈ axes(A,2)
                C[m,n] += A[m,k] * B[k,n]
            end
        end
    end
    # AmulB_avx2!(C, A, B)
    # gq = :(for m ∈ axes(A,1), n ∈ axes(B,2)
    # C[m,n] = z
    # for k ∈ axes(A,2)
    # C[m,n] += A[m,k] * B[k,n]
    # end
    # end);
    # ls = LoopVectorization.loopset(gq);
    # # ls.preamble_symsym
    # ls.operations[1]
    function AmulB_avx3!(C, A, B)
        Kmin = first(axes(A,2)); Kmax = last(axes(A,2))
        @_avx unroll=(2,2) for m ∈ axes(A,1), n ∈ axes(B,2)
            C[m,n] = zero(eltype(C))
            for k ∈ Kmin:Kmax
                C[m,n] += A[m,k] * B[k,n]
            end
        end
    end
    # function AmulB_avx4!(C, A, B)
    #     @_avx for m ∈ axes(A,1), n ∈ axes(B,2)
    #         C[m,n] = myzero(C)
    #         for k ∈ axes(A,2)
    #             C[m,n] += A[m,k] * B[k,n]
    #         end
    #     end
    # end
    # q = :(for m ∈ axes(A,1), n ∈ axes(B,2)
    #       C[m,n] = myzero(C)
    #       for k ∈ axes(A,2)
    #       C[m,n] += A[m,k] * B[k,n]
    #       end
    #       end)
    # ls = LoopVectorization.loopset(q);
    function AmuladdB_avx!(C, A, B, factor = 1)
        @_avx unroll=(2,2) for m ∈ axes(A,1), n ∈ axes(B,2)
            ΔCₘₙ = zero(eltype(C))
            for k ∈ axes(A,2)
                ΔCₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] += ΔCₘₙ * factor
        end
    end

    function AmulB2x2avx!(C, A, B)
        @avx unroll=(2,2) for m ∈ axes(A,1), n ∈ axes(B,2)
            ΔCₘₙ = zero(eltype(C))
            for k ∈ axes(A,2)
                ΔCₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = ΔCₘₙ
        end
    end
    function AmulB2x2_avx!(C, A, B)
        @_avx unroll=(2,2) for m ∈ axes(A,1), n ∈ axes(B,2)
            ΔCₘₙ = zero(eltype(C))
            for k ∈ axes(A,2)
                ΔCₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = ΔCₘₙ
        end
    end

    # function AtmulB!(C, A, B)
    #     for j ∈ axes(C,2), i ∈ axes(C,1)
    #         Cᵢⱼ = zero(eltype(C))
    #         @simd ivdep for k ∈ axes(A,1)
    #             @inbounds Cᵢⱼ += A[k,i] * B[k,j]
    #         end
    #         C[i,j] = Cᵢⱼ
    #     end
    # end
    AtmulBq = :(for n ∈ axes(C,2), m ∈ axes(C,1)
                Cₘₙ = zero(eltype(C))
                for k ∈ axes(A,1)
                Cₘₙ += A[k,m] * B[k,n]
                end
                C[m,n] = Cₘₙ
                end)
    lsAtmulB = LoopVectorization.loopset(AtmulBq);
    if LoopVectorization.register_count() != 8
        @test LoopVectorization.choose_order(lsAtmulB) == (Symbol[:n,:m,:k], :n, :m, :k, Unumt, Tnumt)
    end
    function AtmulBavx1!(C, A, B)
        @avx for n ∈ axes(C,2), m ∈ axes(C,1)
            Cₘₙ = zero(eltype(C))
            for k ∈ axes(A,1)
                Cₘₙ += A[k,m] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end
    Atq = :(for n ∈ axes(C,2), m ∈ axes(C,1)
            Cₘₙ = zero(eltype(C))
            for k ∈ axes(A,1)
            Cₘₙ += A[k,m] * B[k,n]
            end
            C[m,n] += Cₘₙ * factor
            end);
    atls = LoopVectorization.loopset(Atq);
    # LoopVectorization.operations(atls)
    # LoopVectorization.loopdependencies.(operations(atls))
    # LoopVectorization.reduceddependencies.(operations(atls))
    function AtmulB_avx1!(C, A, B)
        @_avx unroll=(2,2) for n ∈ axes(C,2), m ∈ axes(C,1)
            Cₘₙ = zero(eltype(C))
            for k ∈ axes(A,1)
                Cₘₙ += A[k,m] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end
    function AtmulBavx2!(C, A, B)
        M, N = size(C); K = size(B,1)
        @assert size(C, 1) == size(A, 2)
        @assert size(C, 2) == size(B, 2)
        @assert size(A, 1) == size(B, 1)
        # When the @avx macro is available, this code is faster:
        z = zero(eltype(C))
        @avx unroll=(2,2) for n in axes(C,2), m in axes(C,1)
            Cmn = z
            for k in axes(A,1)
                Cmn += A[k,m] * B[k,n]
            end
            C[m,n] = Cmn
        end
        return C
    end
    function AtmulB_avx2!(C, A, B)
        M, N = size(C); K = size(B,1)
        @assert size(C, 1) == size(A, 2)
        @assert size(C, 2) == size(B, 2)
        @assert size(A, 1) == size(B, 1)
        # When the @avx macro is available, this code is faster:
        z = zero(eltype(C))
        @_avx for n in axes(C,2), m in axes(C,1)
            Cmn = z
            for k in axes(A,1)
                Cmn += A[k,m] * B[k,n]
            end
            C[m,n] = Cmn
        end
        return C
    end
    function rank2AmulB!(C, Aₘ, Aₖ, B)
        @inbounds for m ∈ axes(C,1), n ∈ axes(C,2)
            Cₘₙ = zero(eltype(C))
            @fastmath for k ∈ axes(B,1)
                Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end
    r2ambq = :(for m ∈ axes(C,1), n ∈ axes(C,2)
               Cₘₙ = zero(eltype(C))
               for k ∈ axes(B,1)
               Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
               end
               C[m,n] = Cₘₙ
               end);
    lsr2amb = LoopVectorization.loopset(r2ambq);
    if LoopVectorization.register_count() == 32
        @test LoopVectorization.choose_order(lsr2amb) == ([:n, :m, :k], :m, :n, :m, 3, 7)
    elseif LoopVectorization.register_count() == 16
        # @test LoopVectorization.choose_order(lsr2amb) == ([:m, :n, :k], :m, :n, :m, 1, 6)
        @test LoopVectorization.choose_order(lsr2amb) == ([:m, :n, :k], :m, :n, :m, 2, 4)
    end
    function rank2AmulBavx!(C, Aₘ, Aₖ, B)
        @avx for m ∈ axes(C,1), n ∈ axes(C,2)
            Cₘₙ = zero(eltype(C))
            for k ∈ axes(B,1)
                Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end
    function rank2AmulB_avx!(C, Aₘ, Aₖ, B)
        @_avx for m ∈ axes(C,1), n ∈ axes(C,2)
            Cₘₙ = zero(eltype(C))
            for k ∈ axes(B,1)
                Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end
    function rank2AmulBavx_noinline!(C, Aₘ, Aₖ, B)
        @avx inline=false for m ∈ axes(C,1), n ∈ axes(C,2)
            Cₘₙ = zero(eltype(C))
            for k ∈ axes(B,1)
                Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end

    function mulCAtB_2x2blockavx!(C, A, B)
        M, N = size(C); K = size(B,1)
        @assert size(C, 1) == size(A, 2)
        @assert size(C, 2) == size(B, 2)
        @assert size(A, 1) == size(B, 1)
        T = eltype(C)
	for m ∈ 1:2:(M & -2)
	    m1 = m + 1
	    for n ∈ 1:2:(N & -2)
	    	n1 = n + 1
		C11, C21, C12, C22 = zero(T), zero(T), zero(T), zero(T)
		@avx inline=true for k ∈ 1:K
		    C11 += A[k,m] * B[k,n]
		    C21 += A[k,m1] * B[k,n]
		    C12 += A[k,m] * B[k,n1]
		    C22 += A[k,m1] * B[k,n1]
		end
		C[m,n] = C11
		C[m1,n] = C21
		C[m,n1] = C12
		C[m1,n1] = C22
	    end
            if isodd(N)
	    	C1n = 0.0
	    	C2n = 0.0
	    	@avx inline=true for k ∈ 1:K
	    	    C1n += A[k,m] * B[k,N]
	    	    C2n += A[k,m1] * B[k,N]
	    	end
	    	C[m,N] = C1n
	    	C[m1,N] = C2n
            end
        end
        if isodd(M)
	    for n ∈ 1:2:(N & -2)
	    	n1 = n + 1
		Cm1, Cm2 = zero(T), zero(T)
		@avx inline=true for k ∈ 1:K
		    Cm1 += A[k,M] * B[k,n]
		    Cm2 += A[k,M] * B[k,n1]
		end
		C[M,n] = Cm1
		C[M,n1] = Cm2
	    end
            if isodd(N)
	    	Cmn = 0.0
	    	@avx inline=true for k ∈ 1:K
	    	    Cmn += A[k,M] * B[k,N]
	    	end
	    	C[M,N] = Cmn
            end
        end
        return C
    end
    function mulCAtB_2x2block_avx!(C, A, B)
        M, N = size(C); K = size(B,1)
        @assert size(C, 1) == size(A, 2)
        @assert size(C, 2) == size(B, 2)
        @assert size(A, 1) == size(B, 1)
        T = eltype(C)
	for m ∈ 1:2:(M & -2)
	    m1 = m + 1
	    for n ∈ 1:2:(N & -2)
	    	n1 = n + 1
		C11, C21, C12, C22 = zero(T), zero(T), zero(T), zero(T)
		@_avx for k ∈ 1:K
		    C11 += A[k,m] * B[k,n]
		    C21 += A[k,m1] * B[k,n]
		    C12 += A[k,m] * B[k,n1]
		    C22 += A[k,m1] * B[k,n1]
		end
		C[m,n] = C11
		C[m1,n] = C21
		C[m,n1] = C12
		C[m1,n1] = C22
	    end
            if isodd(N)
	    	C1n = 0.0
	    	C2n = 0.0
	    	@_avx for k ∈ 1:K
	    	    C1n += A[k,m] * B[k,N]
	    	    C2n += A[k,m1] * B[k,N]
	    	end
	    	C[m,N] = C1n
	    	C[m1,N] = C2n
            end
        end
        if isodd(M)
	    for n ∈ 1:2:(N & -2)
	    	n1 = n + 1
		Cm1, Cm2 = zero(T), zero(T)
		@_avx for k ∈ 1:K
		    Cm1 += A[k,M] * B[k,n]
		    Cm2 += A[k,M] * B[k,n1]
		end
		C[M,n] = Cm1
		C[M,n1] = Cm2
	    end
            if isodd(N)
	    	Cmn = 0.0
	    	@_avx for k ∈ 1:K
	    	    Cmn += A[k,M] * B[k,N]
	    	end
	    	C[M,N] = Cmn
            end
        end
        return C
    end

    function mulCAtB_2x2blockavx_noinline!(C, A, B)
        M, N = size(C); K = size(B,1)
        @assert size(C, 1) == size(A, 2)
        @assert size(C, 2) == size(B, 2)
        @assert size(A, 1) == size(B, 1)
        T = eltype(C)
	for m ∈ 1:2:(M & -2)
	    m1 = m + 1
	    for n ∈ 1:2:(N & -2)
	    	n1 = n + 1
		C11, C21, C12, C22 = zero(T), zero(T), zero(T), zero(T)
		@avx inline=false for k ∈ 1:K
		    C11 += A[k,m] * B[k,n]
		    C21 += A[k,m1] * B[k,n]
		    C12 += A[k,m] * B[k,n1]
		    C22 += A[k,m1] * B[k,n1]
		end
		C[m,n] = C11
		C[m1,n] = C21
		C[m,n1] = C12
		C[m1,n1] = C22
	    end
            if isodd(N)
	    	C1n = 0.0
	    	C2n = 0.0
	    	@avx inline=false for k ∈ 1:K
	    	    C1n += A[k,m] * B[k,N]
	    	    C2n += A[k,m1] * B[k,N]
	    	end
	    	C[m,N] = C1n
	    	C[m1,N] = C2n
            end
        end
        if isodd(M)
	    for n ∈ 1:2:(N & -2)
	    	n1 = n + 1
		Cm1, Cm2 = zero(T), zero(T)
		@avx inline=false for k ∈ 1:K
		    Cm1 += A[k,M] * B[k,n]
		    Cm2 += A[k,M] * B[k,n1]
		end
		C[M,n] = Cm1
		C[M,n1] = Cm2
	    end
            if isodd(N)
	    	Cmn = 0.0
	    	@avx inline=false for k ∈ 1:K
	    	    Cmn += A[k,M] * B[k,N]
	    	end
	    	C[M,N] = Cmn
            end
        end
        return C
    end

    # TODO: add fast=false option to `@avx`
    # function gemm_accurate!(C, A, B)
    #     @avx for n in axes(C,2), m in axes(C,1)
    #         Cmn_hi = zero(eltype(C))
    #         Cmn_lo = zero(eltype(C))
    #         for k in axes(B,1)
    #             hiprod = vmul(A[m,k], B[k,n])
    #             loprod = vfmsub(A[m,k], B[k,n], hiprod)
    #             hi_ts = vadd(hiprod, Cmn_hi)
    #             a1_ts = vsub(hi_ts, Cmn_hi)
    #             b1_ts = vsub(hi_ts, a1_ts)
    #             lo_ts = vadd(vsub(hiprod, a1_ts), vsub(Cmn_hi, b1_ts))
    #             thi = vadd(loprod, Cmn_lo)
    #             a1_t = vsub(thi, Cmn_lo)
    #             b1_t = vsub(thi, a1_t)
    #             tlo = vadd(vsub(loprod, a1_t), vsub(Cmn_lo, b1_t))
    #             c1 = vadd(lo_ts, thi)
    #             hi_ths = vadd(hi_ts, c1)
    #             lo_ths = vsub(c1, vsub(hi_ths, hi_ts))
    #             c2 = vadd(tlo, lo_ths)
    #             Cmn_hi = vadd(hi_ths, c2)
    #             Cmn_lo = vsub(c2, vsub(Cmn_hi, hi_ths))
    #         end
    #         C[m,n] = Cmn_hi
    #     end
    # end
    function AB_plus_BA!(du, u, mat)
        @assert size(u, 1) == size(u, 2) == size(mat, 1) == size(mat, 2)
        for i2 in 1:size(u, 2), i1 in 1:size(u, 1)
            for sum_idx in 1:size(u, 1)
                du[i1, i2] += mat[i1, sum_idx] * u[sum_idx, i2] + mat[i2, sum_idx] * u[i1, sum_idx]
            end
        end
        return nothing
    end

    function AB_plus_BA_avx!(du, u, mat)
        @assert size(u, 1) == size(u, 2) == size(mat, 1) == size(mat, 2)
        @avx for i2 in 1:size(u, 2), i1 in 1:size(u, 1)
            for sum_idx in 1:size(u, 1)
                du[i1, i2] += mat[i1, sum_idx] * u[sum_idx, i2] + mat[i2, sum_idx] * u[i1, sum_idx]
            end
        end
        return nothing
    end

    function threegemms!(Ab, Bb, Cb, A, B, C)
        M, N = size(Cb); K = size(B,1)
        @avx for m in 1:M, k in 1:K, n in 1:N
            Ab[m,k] += C[m,n] * B[k,n]
            Bb[k,n] += A[m,k] * C[m,n]
            Cb[m,n] += A[m,k] * B[k,n]
        end
    end
    # M = 77;
    # A = rand(M,M); B = rand(M,M); C = similar(A);
    # mulCAtB_2x2block_avx!(C,A,B)
    # using LoopVectorization
    # mul2x2q = :(for k ∈ 1:K
    # C11 += A[k,m] * B[k,n]
    # C21 += A[k,m1] * B[k,n]
    # C12 += A[k,m] * B[k,n1]
    # C22 += A[k,m1] * B[k,n1]
    # end)
    # lsmul2x2q = LoopVectorization.loopset(mul2x2q)

    lsAtmulBt8 = :(for m ∈ 1:8, n ∈ 1:8
    ΔCₘₙ = zero(eltype(C))
    for k ∈ 1:8
    ΔCₘₙ += A[k,m] * B[n,k]
    end
    C[m,n] = ΔCₘₙ
    end) |> LoopVectorization.loopset;
    if LoopVectorization.register_count() == 32
        # @test LoopVectorization.choose_order(lsAtmulBt8) == ([:n, :m, :k], :m, :n, :m, 1, 8)
        @test LoopVectorization.choose_order(lsAtmulBt8) == ([:n, :m, :k], :k, :n, :m, 1, 8)
    elseif LoopVectorization.register_count() == 16
        # @test LoopVectorization.choose_order(lsAtmulBt8) == ([:n, :m, :k], :m, :n, :m, 2, 4)
        @test LoopVectorization.choose_order(lsAtmulBt8) == ([:n, :m, :k], :n, :m, :n, 2, 4)
    elseif LoopVectorization.register_count() == 8
        @test LoopVectorization.choose_order(lsAtmulBt8) == ([:n, :m, :k], :m, :n, :m, 1, 4)
    end

    struct SizedMatrix{M,N,T} <: DenseMatrix{T}
        data::Matrix{T}
        function SizedMatrix{M,N}(data::Matrix{T}) where {M,N,T}
            @assert (M,N) === size(data)
            new{M,N,T}(data)
        end
    end
    Base.parent(A::SizedMatrix) = A.data
    Base.IndexStyle(::Type{<:SizedMatrix}) = Base.IndexLinear()
    Base.@propagate_inbounds Base.getindex(A::SizedMatrix, i::Int) = getindex(parent(A), i)
    Base.@propagate_inbounds Base.setindex!(A::SizedMatrix, v, i::Int) = setindex!(parent(A), v, i)
    Base.@propagate_inbounds Base.getindex(A::SizedMatrix, i::CartesianIndex) = getindex(parent(A), i + oneunit(i))
    Base.@propagate_inbounds Base.setindex!(A::SizedMatrix, v, i::CartesianIndex) = setindex!(parent(A), v, i + oneunit(i))
    Base.@propagate_inbounds Base.getindex(A::SizedMatrix, i::Int, j::Int) = getindex(parent(A), i+1, j+1)
    Base.@propagate_inbounds Base.setindex!(A::SizedMatrix, v, i::Int, j::Int) = setindex!(parent(A), v, i+1, j+1)
    Base.size(::SizedMatrix{M,N}) where {M,N} = (M,N)
    LoopVectorization.ArrayInterface.size(::SizedMatrix{M,N}) where {M,N} = (LoopVectorization.Static{M}(),LoopVectorization.Static{N}())
    function Base.axes(::SizedMatrix{M,N}) where {M,N}
        (LoopVectorization.CloseOpen(LoopVectorization.Static{M}()), LoopVectorization.CloseOpen(LoopVectorization.Static{N}()))
    end
    function LoopVectorization.ArrayInterface.axes_types(::Type{SizedMatrix{M,N,T}}) where {M,N,T}
        Tuple{LoopVectorization.CloseOpen{LoopVectorization.Static{0},LoopVectorization.Static{M}}, LoopVectorization.CloseOpen{LoopVectorization.Static{0},LoopVectorization.Static{N}}}
    end
    Base.unsafe_convert(::Type{Ptr{T}}, A::SizedMatrix{M,N,T}) where {M,N,T} = pointer(A.data)
    LoopVectorization.ArrayInterface.strides(::SizedMatrix{M}) where {M} = (LoopVectorization.Static{1}(),LoopVectorization.Static{M}())
    LoopVectorization.ArrayInterface.contiguous_axis(::Type{<:SizedMatrix}) = LoopVectorization.One()
    LoopVectorization.ArrayInterface.contiguous_batch_size(::Type{<:SizedMatrix}) = LoopVectorization.Zero()
    LoopVectorization.ArrayInterface.stride_rank(::Type{<:SizedMatrix}) = (LoopVectorization.Static(1), LoopVectorization.Static(2))
    # LoopVectorization.ArrayInterface.offsets(::Type{SizedMatrix{M,N,T}}) where {M,N,T}  = (LoopVectorization.Static{0}(), LoopVectorization.Static{0}())
    LoopVectorization.ArrayInterface.offsets(::SizedMatrix) = (LoopVectorization.Static{0}(), LoopVectorization.Static{0}())
    LoopVectorization.ArrayInterface.dense_dims(::Type{SizedMatrix{M,N,T}}) where {M,N,T} = LoopVectorization.ArrayInterface.dense_dims(Matrix{T})
# struct ZeroInitializedArray{T,N,A<:DenseArray{T,N}} <: DenseArray{T,N}
#     data::A
# end
# Base.size(A::ZeroInitializedArray) = size(A.data)
# Base.length(A::ZeroInitializedArray) = length(A.data)
# Base.axes(A::ZeroInitializedArray, i) = axes(A.data, i)
# @inline Base.getindex(A::ZeroInitializedArray{T}) where {T} = zero(T)
# Base.@propagate_inbounds Base.setindex!(A::ZeroInitializedArray, v, i...) = setindex!(A.data, v, i...)
# function LoopVectorization.VectorizationBase.stridedpointer(A::ZeroInitializedArray)
#     LoopVectorization.VectorizationBase.ZeroInitializedStridedPointer(LoopVectorization.VectorizationBase.stridedpointer(A.data))
# end

    for T ∈ (Float32, Float64, Int32, Int64)
        TC = sizeof(T) == 4 ? Float32 : Float64
        R = T <: Integer ? (T(-1000):T(1000)) : T
        for M ∈ 48:54
            C0 = zeros(TC, M, M); C1 = zeros(TC, M, M);
            A = rand(R, M, M); B = rand(R, M, M);
            AB_plus_BA!(C0, A, B)
            AB_plus_BA_avx!(C1, A, B)
            @test C0 ≈ C1
        end
    # let T = Int32
        # exceeds_time_limit() && break
        @show T, @__LINE__
        # M, K, N = 128, 128, 128;
        N = 69;
        for M ∈ 72:80, K ∈ 72:80
        # M, K, N = 73, 75, 69;
            C = Matrix{TC}(undef, M, N);
            A = rand(R, M, K); B = rand(R, K, N);
            At = copy(A');
            Bt = copy(B');
            C2 = similar(C);
            @time @testset "avx $T dynamc gemm" begin
                AmulB!(C2, A, B)
                AmulBavx1!(C, A, B)
                @test C ≈ C2
                fill!(C, 999.99); AmulBavx1!(C, At', B)
                @test C ≈ C2
                fill!(C, 999.99); AmulBavx2!(C, A, B)
                @test C ≈ C2
                fill!(C, 999.99); AmulBavx2!(C, At', B)
                @test C ≈ C2
                fill!(C, 999.99); AmulBavx2!(C, A, Bt')
                @test C ≈ C2
                fill!(C, 999.99); AmulBavx2!(C, At', Bt')
                @test C ≈ C2
                fill!(C, 999.99); AmulBavx3!(C, A, B)
                @test C ≈ C2
                fill!(C, 999.99); AmulBavx3!(C, At', B)
                @test C ≈ C2
                fill!(C, 0.0); AmuladdBavx!(C, A, B)
                @test C ≈ C2
                AmuladdBavx!(C, At', B)
                @test C ≈ 2C2
                AmuladdBavx!(C, A, B, -1)
                @test C ≈ C2
                AmuladdBavx!(C, At', Bt', -2)
                @test C ≈ -C2
                AmuladdBavx!(C, At', B, 3, 2)
                @test C ≈ C2
                # How much of this can I do before rounding errors are likely to cause test failures?
                # Setting back to zero here...
                AmuladdBavx!(C, At', B, 1, 0)
                @test C ≈ C2
                AmuladdBavx!(C, At', Bt', 2, -1)
                @test C ≈ C2
                # TODO: Reimplement the ZeroInitialized wrappers
                # fill!(C, 9999.999); AmuladdBavx!(ZeroInitializedArray(C), At', Bt')
                # @test C ≈ C2
                fill!(C, 9999.999); AmulB2x2avx!(C, A, B);
                @test C ≈ C2
                fill!(C, 9999.999); AmulB2x2avx!(C, At', B);
                @test C ≈ C2
                fill!(C, 9999.999); AtmulBavx1!(C, At, B)
                @test C ≈ C2
                fill!(C, 9999.999); AtmulBavx1!(C, A', B)
                @test C ≈ C2
                fill!(C, 9999.999); AtmulBavx2!(C, At, B);
                @test C ≈ C2
                fill!(C, 9999.999); AtmulBavx2!(C, A', B);
                @test C ≈ C2
                fill!(C, 9999.999); mulCAtB_2x2blockavx!(C, At, B);
                @test C ≈ C2
                fill!(C, 9999.999); mulCAtB_2x2blockavx!(C, A', B);
                @test C ≈ C2
                fill!(C, 9999.999); mulCAtB_2x2blockavx_noinline!(C, At, B);
                @test C ≈ C2
                fill!(C, 9999.999); mulCAtB_2x2blockavx_noinline!(C, A', B);
                @test C ≈ C2
                if RUN_SLOW_TESTS
                    # fill!(C, 9999.999); gemm_accurate!(C, A, B);
                    # @test C ≈ C2
                    # fill!(C, 9999.999); gemm_accurate!(C, At', B);
                    # @test C ≈ C2
                    # fill!(C, 9999.999); gemm_accurate!(C, A, Bt');
                    # @test C ≈ C2
                    # fill!(C, 9999.999); gemm_accurate!(C, At', Bt');
                    # @test C ≈ C2
                    Ab = zeros(eltype(C), size(A)); Bb = zeros(eltype(C), size(B)); Cb = zero(C);
                    threegemms!(Ab, Bb, Cb, A, B, C)
                    @test Ab ≈ C * B'
                    @test Bb ≈ A' * C
                    @test Cb ≈ A * B
                end
                if iszero(size(A,1) % 8)
                    Abit = A .> 0.5;
                    fill!(C, 9999.999); AmulBavx1!(C, Abit, B);
                    @test C ≈ Abit * B
                end
                if iszero(size(B,1) % 8)
                    Bbit = B .> 0.5;
                    fill!(C, 9999.999); AmulBavx1!(C, A, Bbit);
                    @test C ≈ A * Bbit
                end
            end
            # exceeds_time_limit() && break
            @time @testset "_avx $T dynamic gemm" begin
                AmulB_avx1!(C, A, B)
                @test C ≈ C2
                fill!(C, 999.99); AmulB_avx1!(C, At', B)
                @test C ≈ C2
                fill!(C, 999.99); AmulB_avx2!(C, A, B)
                @test C ≈ C2
                fill!(C, 999.99); AmulB_avx2!(C, At', B)
                @test C ≈ C2
                fill!(C, 999.99); AmulB_avx2!(C, A, Bt')
                @test C ≈ C2
                fill!(C, 999.99); AmulB_avx2!(C, At', Bt')
                @test C ≈ C2
                fill!(C, 999.99); AmulB_avx3!(C, A, B)
                @test C ≈ C2
                fill!(C, 999.99); AmulB_avx3!(C, At', B)
                @test C ≈ C2
                fill!(C, 0.0); AmuladdB_avx!(C, A, B)
                @test C ≈ C2
                AmuladdB_avx!(C, At', B)
                @test C ≈ 2C2
                AmuladdB_avx!(C, A, B, -1)
                @test C ≈ C2
                AmuladdB_avx!(C, At', B, -2)
                @test C ≈ -C2
                fill!(C, 9999.999); AmulB2x2_avx!(C, A, B);
                @test C ≈ C2
                fill!(C, 9999.999); AmulB2x2_avx!(C, At', B);
                @test C ≈ C2
                fill!(C, 9999.999); AtmulB_avx1!(C, At, B);
                @test C ≈ C2
                fill!(C, 9999.999); AtmulB_avx1!(C, A', B);
                @test C ≈ C2
                fill!(C, 9999.999); AtmulB_avx2!(C, At, B);
                @test C ≈ C2
                fill!(C, 9999.999); AtmulB_avx2!(C, A', B);
                @test C ≈ C2
                fill!(C, 9999.999); mulCAtB_2x2block_avx!(C, At, B);
                @test C ≈ C2
                fill!(C, 9999.999); mulCAtB_2x2block_avx!(C, A', B);
                @test C ≈ C2
            end
            # exceeds_time_limit() && break
            if (M,K) === (73,77) # pick a random size, we only want to compile once
                As = SizedMatrix{M,K}(A);
                Ats = SizedMatrix{K,M}(At);
                Bs = SizedMatrix{K,N}(B);
                Bts = SizedMatrix{N,K}(Bt);
                Cs = SizedMatrix{M,N}(C);
                C2z = LoopVectorization.OffsetArray(C2, -1, -1);
                @time @testset "avx $T static gemm" begin
                    # AmulBavx1!(Cs, As, Bs)
                    # @test Cs ≈ C2
                    # fill!(Cs, 999.99); AmulBavx1!(Cs, Ats', Bs)
                    # @test Cs ≈ C2
                    fill!(Cs, 999.99); AmulBavx2!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulBavx2!(Cs, Ats', Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulBavx2!(Cs, As, Bts')
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulBavx2!(Cs, Ats', Bts')
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulBavx3!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulBavx3!(Cs, Ats', Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 0.0); AmuladdBavx!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    AmuladdBavx!(Cs, Ats', Bs)
                    @test Cs ≈ 2C2z
                    AmuladdBavx!(Cs, As, Bs, -1)
                    @test Cs ≈ C2z
                    AmuladdBavx!(Cs, Ats', Bs, -2)
                    @test Cs ≈ -C2z
                    fill!(Cs, 9999.999); AmulB2x2avx!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AmulB2x2avx!(Cs, Ats', Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulBavx1!(Cs, Ats, Bs);
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulBavx1!(Cs, As', Bs);
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulBavx2!(Cs, Ats, Bs);
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulBavx2!(Cs, As', Bs);
                    @test Cs ≈ C2z
                    # fill!(Cs, 9999.999); mulCAtB_2x2blockavx!(Cs, Ats, Bs);
                    # @test Cs ≈ C2
                    # fill!(Cs, 9999.999); mulCAtB_2x2blockavx!(Cs, As', Bs);
                    # @test Cs ≈ C2
                    # fill!(Cs, 9999.999); mulCAtB_2x2blockavx_noinline!(Cs, Ats, Bs);
                    # @test Cs ≈ C2
                    # fill!(Cs, 9999.999); mulCAtB_2x2blockavx_noinline!(Cs, As', Bs);
                    # @test Cs ≈ C2
                end
                # exceeds_time_limit() && break
                @time @testset "_avx $T static gemm" begin
                    # AmulB_avx1!(Cs, As, Bs)
                    # @test Cs ≈ C2
                    # fill!(Cs, 999.99); AmulB_avx1!(Cs, Ats', Bs)
                    # @test Cs ≈ C2
                    fill!(Cs, 999.99); AmulB_avx2!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulB_avx2!(Cs, Ats', Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulB_avx2!(Cs, As, Bts')
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulB_avx2!(Cs, Ats', Bts')
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulB_avx3!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 999.99); AmulB_avx3!(Cs, Ats', Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 0.0); AmuladdB_avx!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    AmuladdB_avx!(Cs, Ats', Bs)
                    @test Cs ≈ 2C2z
                    AmuladdB_avx!(Cs, As, Bs, -1)
                    @test Cs ≈ C2z
                    AmuladdB_avx!(Cs, Ats', Bs, -2)
                    @test Cs ≈ -C2z
                    fill!(Cs, 9999.999); AmulB2x2_avx!(Cs, As, Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AmulB2x2_avx!(Cs, Ats', Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulB_avx1!(Cs, Ats, Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulB_avx1!(Cs, As', Bs)
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulB_avx2!(Cs, Ats, Bs);
                    @test Cs ≈ C2z
                    fill!(Cs, 9999.999); AtmulB_avx2!(Cs, As', Bs);
                    @test Cs ≈ C2z
                    # fill!(Cs, 9999.999); mulCAtB_2x2block_avx!(Cs, Ats, Bs);
                    # @test Cs ≈ C2
                    # fill!(Cs, 9999.999); mulCAtB_2x2block_avx!(Cs, As', Bs);
                    # @test Cs ≈ C2
                end
            end
            # exceeds_time_limit() && break
            @time @testset "$T rank2mul" begin
                Aₘ= rand(R, M, 2); Aₖ = rand(R, 2, K);
                Aₖ′ = copy(Aₖ');
                rank2AmulB!(C2, Aₘ, Aₖ, B)
                rank2AmulBavx!(C, Aₘ, Aₖ, B)
                @test C ≈ C2
                fill!(C, 9999.999); rank2AmulB_avx!(C, Aₘ, Aₖ, B)
                @test C ≈ C2
                fill!(C, 9999.999); rank2AmulBavx_noinline!(C, Aₘ, Aₖ, B)
                @test C ≈ C2
                fill!(C, 9999.999); rank2AmulBavx!(C, Aₘ, Aₖ′', B)
                @test C ≈ C2
                fill!(C, 9999.999); rank2AmulB_avx!(C, Aₘ, Aₖ′', B)
                @test C ≈ C2
                fill!(C, 9999.999); rank2AmulBavx_noinline!(C, Aₘ, Aₖ′', B)
                @test C ≈ C2
            end
        end
    end
end


text:
# ALL TESTS THAT CONCERN UNITCELLS

# include
using LatticePhysics
using Base.Test

# begin the Unitcell testblock
SVG_plotting_testset = @testset "SVG Plotting of lattices" begin

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT LATTICE PLOTTING FUNCTIONS FOR
#   PLOTTING TO SVG IMAGES
#
################################################################################

# 2D
@testset "2D Plotting" begin

    # get a test unitcell
    unitcell = getUnitcellHoneycomb(4)
    # get a lattice
    lattice = getLattice(unitcell, [-10, -5])

    # TODO naive method plotLattice

    # Test naive plotting
    @test typeof(plotLattice2D(lattice)) == String
    # Naive plotting and Inkscape pdf export
    @test typeof(plotLattice2D(lattice, inkscape_export_pdf=true)) == String
    # Naive plotting and printing options
    @test typeof(plotLattice2D(lattice, print_used_options=true)) == String

    # custom canvas border size and size
    @test typeof(plotLattice2D(lattice, border_percentage=0.02, size_long_side=2000)) == String

    # periodic connections
    @test typeof(plotLattice2D(lattice, visualize_periodic=true)) == String
    # certain bond thickness
    @test typeof(plotLattice2D(lattice, bond_thickness=-0.2)) == String
    @test typeof(plotLattice2D(lattice, bond_thickness=10.0)) == String
    # certain bond colors
    @test typeof(plotLattice2D(lattice, colorcode_bonds_automation="GREY")) == String
    @test typeof(plotLattice2D(lattice, colorcode_bonds_automation="COLOR")) == String
    @test typeof(plotLattice2D(lattice, colorcode_bonds=Dict("tx"=>[155,0,0], "ty"=>[0,155,0], "tz"=>[0,0,155]))) == String

    # custom site colors
    @test typeof(plotLattice2D(lattice, colorcode_sites=Dict(1=>[155,0,50], 2=>[0,155,50]))) == String
    # custom site labels
    @test typeof(plotLattice2D(lattice, site_labels="POSITION INDEX")) == String
    @test typeof(plotLattice2D(lattice, site_labels="LATTICE INDEX")) == String
    # custom site radius
    @test typeof(plotLattice2D(lattice, site_radius=40.0)) == String
    @test typeof(plotLattice2D(lattice, site_radius=-0.3)) == String
    # cuistom site border
    @test typeof(plotLattice2D(lattice, site_border_width_percentage=0.2)) == String

# end the testset here
end;

# 3D
@testset "3D Plotting" begin

    # get a test unitcell
    unitcell = getUnitcell_10_3_a(4)
    # get a lattice
    lattice = getLattice(unitcell, [-2, -2, -2])

# end the testset here
end;

# 2D Plaquettes
@testset "2D Plaquette Plotting" begin

    # get a test unitcell
    unitcell = getUnitcellHoneycomb(4)
    # get a lattice
    lattice = getLatticeInBox(unitcell, [10.0, 5.0])
    # get the plaquettes
    plaquettes = getPlaquettesOfSite(lattice, 3, 6)
    # get the plaquette values
    plaquette_values = [rand()>0.5 ? 1.0 : -1.0 for p in plaquettes]

    # PLOTTING RELATED (SAME AS 2D PLOTTING)

    # Test naive plotting
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values)) == String
    # Naive plotting and Inkscape pdf export
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, inkscape_export_pdf=true)) == String
    # Naive plotting and printing options
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, print_used_options=true)) == String

    # custom canvas border size and size
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, border_percentage=0.02, size_long_side=2000)) == String

    # periodic connections
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, visualize_periodic=true)) == String
    # certain bond thickness
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, bond_thickness=-0.2)) == String
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, bond_thickness=10.0)) == String
    # certain bond colors
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, colorcode_bonds_automation="GREY")) == String
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, colorcode_bonds_automation="COLOR")) == String
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, colorcode_bonds=Dict("tx"=>[155,0,0], "ty"=>[0,155,0], "tz"=>[0,0,155]))) == String

    # custom site colors
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, colorcode_sites=Dict(1=>[155,0,50], 2=>[0,155,50]))) == String
    # custom site labels
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, site_labels="POSITION INDEX")) == String
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, site_labels="LATTICE INDEX")) == String
    # custom site radius
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, site_radius=40.0)) == String
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, site_radius=-0.3)) == String
    # cuistom site border
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, site_border_width_percentage=0.2)) == String

    # PLAQUETTE RELATED

    # cuistom colorcode
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, colorcode_plaquettes=Dict(-1.0=>[0,125,0], 1.0=>[255,0,125]))) == String
    # cuistom opacity
    @test typeof(plotPlaquettes2D(lattice, plaquettes, plaquette_values, opacity_plaquettes=0.1)) == String

# end the testset here
end;

# end the SVG plotting testset here
end;


text:
#using Test

#=
Test rpn operations:
=#

@testset "RPN" begin
end


text:
module RPRMakie

using Makie
using RadeonProRender
using GeometryBasics
using Colors
const RPR = RadeonProRender

include("scene.jl")
include("lines.jl")
include("meshes.jl")

end


text:
@testset "trailer.jl" begin

    @testset "StateObject{Int}()" begin
        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Int}(3, trailer)

        @test reversibleInt.value == 3
        @test reversibleInt.trailer == trailer
    end

    @testset "StateObject{Bool}()" begin
        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Bool}(true, trailer)

        @test reversibleInt.value == true
        @test reversibleInt.trailer == trailer
    end

    @testset "trail!()" begin
        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Int}(3, trailer)

        SeaPearl.trail!(reversibleInt)

        @test length(trailer.current) == 1

        se = first(trailer.current)

        @test se.object == reversibleInt
        @test se.value == 3
    end

    @testset "setValue!()" begin
        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Int}(3, trailer)

        # Check when same value
        SeaPearl.setValue!(reversibleInt, 3)
        @test length(trailer.current) == 0
        @test reversibleInt.value == 3

        # With a different value
        SeaPearl.setValue!(reversibleInt, 5)
        @test length(trailer.current) == 1
        @test reversibleInt.value == 5

        se = first(trailer.current)
        @test se.value == 3

    end

    @testset "saveState!()" begin
        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Int}(3, trailer)

        current = trailer.current

        SeaPearl.setValue!(reversibleInt, 5)
        SeaPearl.saveState!(trailer)

        @test first(trailer.prior) == current
        @test isempty(trailer.current)
    end

    @testset "restoreState!()" begin
        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Int}(3, trailer)

        current = trailer.current

        SeaPearl.setValue!(reversibleInt, 5)
        SeaPearl.saveState!(trailer)
        SeaPearl.setValue!(reversibleInt, 8)

        @test reversibleInt.value == 8

        SeaPearl.restoreState!(trailer)

        @test reversibleInt.value == 5

        SeaPearl.restoreState!(trailer)

        @test reversibleInt.value == 3
    end

    @testset "withNewState!()" begin

        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Int}(3, trailer)

        SeaPearl.withNewState!(trailer) do
            SeaPearl.setValue!(reversibleInt, 5)

            @test reversibleInt.value == 5
        end

        @test reversibleInt.value == 3

    end

    @testset "restoreInitialState!()" begin
        trailer = SeaPearl.Trailer()
        reversibleInt = SeaPearl.StateObject{Int}(3, trailer)

        SeaPearl.saveState!(trailer)
        SeaPearl.setValue!(reversibleInt, 4)
        SeaPearl.saveState!(trailer)
        SeaPearl.setValue!(reversibleInt, 5)
        SeaPearl.saveState!(trailer)
        SeaPearl.setValue!(reversibleInt, 6)
        SeaPearl.saveState!(trailer)
        SeaPearl.setValue!(reversibleInt, 7)

        SeaPearl.restoreInitialState!(trailer)

        @test reversibleInt.value == 3
    end

end


text:
# saturated addition for positive numbers
function addsaturated(a::T,b::T) where {T}
    r = a+b
    r<a || r<b ? typemax(T) : r # check for overflow
end
addsaturated(a::T,b::T) where {T<:AbstractFloat} = a+b # Inf+a = Inf, so no overflow check is needed

# Finds the shortest paths between all pairs of nodes in an directed graph, specified by a matrix of pairwise distances.
# If there are no edges between two nodes, the distance should be typemax(T) (which is Inf for Float64).
function floydwarshall!(D::AbstractMatrix{T}) where {T<:Real}
    @assert ndims(D)==2
    @assert size(D,1)==size(D,2)
    @assert all(D.>=0)
    N = size(D,1)

    for k=1:N
        for j=1:N
            for i=1:N
                #D[j,i] = min(D[j,i], D[j,k]+D[k,i])
                D[j,i] = min(D[j,i], addsaturated(D[j,k],D[k,i]))
            end
        end
    end

    # This is not needed since we assume the matrix has only positive entries
    # for i=1:N
    #     @assert D[i,i]>=0 "The graph has negative cycles."
    # end
    D
end
floydwarshall(D::AbstractMatrix{T}) where {T<:Real} = floydwarshall!(copy(D))


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 39cf2802-71fa-11eb-2ba9-aba632614050
md"""
# machine learning
[project placeholder]
"""

# ╔═╡ Cell order:
# ╟─39cf2802-71fa-11eb-2ba9-aba632614050

max_stars_repo_path: Projects/Project 2 - machine learning.jl