text:
module GHPGHX
# -------
export ghp_model
export get_ghpghx_results_for_reopt

include("ghpghx_inputs.jl")

include("ghpghx_results.jl")

"""
    ghp_model(d)

Process the inputs and run the `size_borefield(p)` function.

Return the processed set of inputs and the results from `size_borefield(p)`.

"""

function ghp_model(d)
    # Load parameters and process from input Dict, to create instance of type InputsStruct
    p = InputsProcess(d)

    # Run size_borefield for each attempt to SOLVE
    results = size_borefield(p)

    return results, p
end

"""
    size_borefield(p)

Determine the size and dispatch of the heat pump and ground heat exchanger to serve the heating and cooling loads.

Return the mutated results struct.

"""
function size_borefield(p)
    # Declare and initialize arrays which get passed and mutated by GHX model (different length arrays used for different models)
    INFO = zeros(Int32, 15)  # Used for initialization and incrementing the number of times the GHX model is called by timestep
    if p.ghx_model == "DST"
        PAR = zeros(Float64, 112)  # Only changed by assign_PAR! each sizing iteration for ResultsStruct variables
        XIN = zeros(Float64, 5)  # Key inputs passed from GHP calcs to the GHX model (T_in, Mdot_in, T_amb)
        OUT = zeros(Float64, 29)  # Outputs from the GHX model to go into the GHP calcs - most notably OUT[1] = Tout_GHX
    end
    if p.ghx_model == "TESS"
        PAR = zeros(Float64, 32)  # Only changed by assign_PAR! each sizing iteration for ResultsStruct variables
        XIN = zeros(Float64, 4)  # Key inputs passed from GHP calcs to the GHX model (T_in, Mdot_in, T_amb)
        OUT = zeros(Float64, 12)  # Outputs from the GHX model to go into the GHP calcs - most notably OUT[1] = Tout_GHX
    end

    # Define TimeArray; TimeArray[1] varies with time, and TimeArray[4] varies for TESS GHX model with variable time step per hour
    Time0 = 0.0
    Time = Time0
    Tfinal = 8760 * p.simulation_years
    Delt = 1.0 / p.dst_ghx_timesteps_per_hour  # This gets overwritten for TESS GHX model for each hour
    TimeArray = [Time, Time0, Tfinal, Delt]

    # Initial sizing guess and initialize ResultsStruct
    X_Now::Array{Float64, 1} = zeros(p.max_sizing_iterations + 1)
    FX_Now::Array{Float64, 1} = zeros(p.max_sizing_iterations + 1)
    N_Bores::Array{Float64, 1} = zeros(p.max_sizing_iterations + 1)

    # Initial size guess for first sizing iteration
    size_iter = 1
    X_Now[size_iter] = p.X_init

    r = ResultsStruct(X_Now = X_Now, FX_Now = FX_Now, N_Bores = N_Bores,
                        total_hours = p.simulation_years * 8760)
    init_sizing!(r, p, size_iter)  # Only modifies r, not p

    #Set the parameters for the ground model
    assign_PAR!(p, r, PAR, size_iter)

    # Instantiate TESS GHX storage/save variables
    #tess_save = TESSStruct()

    # Initialize GHX model with 2 calls
    #init_ghx_calls_2x!(p, TimeArray, XIN, OUT, PAR, INFO, tess_save)
    init_ghx_calls_2x!(p, TimeArray, XIN, OUT, PAR, INFO)

    ###   LOOP THROUGH FOR EACH SIZING ITERATION UNTIL "break" or reaching max iterations   ###
    while size_iter <= p.max_sizing_iterations
        println("size_iter = ", size_iter)
        # Initialize/reset and update variables (already done for first iteration above)
        if size_iter > 1
            r = ResultsStruct(X_Now = r.X_Now, FX_Now = r.FX_Now, N_Bores = r.N_Bores,
                                total_hours = p.simulation_years * 8760)
            init_sizing!(r, p, size_iter)  # Only modifies r, not p
            assign_PAR!(p, r, PAR, size_iter)
            init_ghx_calls_2x!(p, TimeArray, XIN, OUT, PAR, INFO)
        end
        println("N_bores = ", r.N_Bores[size_iter])
        for year in 1:p.simulation_years
            for hr in 1:8760  # TESS model has hourly-based timesteps per hour
                Q_Heat = p.HeatingThermalLoadKW[hr] * 3600.0  # Convert kW to kJ/hr
                Q_Cool = p.CoolingThermalLoadKW[hr] * 3600.0  # Convert kW to kJ/hr
                # Find time step per hour based on flow rate
                Tons_HeatPump_H = Q_Heat / 1.055 / 12000.0
                Tons_HeatPump_C = Q_Cool / 1.055 / 12000.0
                GPM_GHX = max(p.fMin_VSP_GHXPump * p.GPM_GHXPump, (Tons_HeatPump_H + Tons_HeatPump_C) * p.GPMperTon_WSHP)
                Mdot_GHX = GPM_GHX * 60.0 / 264.172 * p.Rho_GHXFluid
                Mdot_Circuit = max(Mdot_GHX / r.N_Bores[size_iter] * p.N_Series, 0.001)
                Circuit_Time = max(r.Mass_Circuit / Mdot_Circuit, 0.001)

                # TESS GHX model, Model = 3, uses variable time step per hour; DST uses fixed p.dst_ghx_timesteps_per_hour
                if p.ghx_model == "TESS"
                   N_Steps_ThisHour = max(p.tess_ghx_minimum_timesteps_per_hour, convert(Int64, floor(1.0 / Circuit_Time)) + 1)
                else
                   N_Steps_ThisHour = p.dst_ghx_timesteps_per_hour
                end

                for step in 1:N_Steps_ThisHour
                    Delt = 1 / N_Steps_ThisHour
                    if year == 1 && hr == 1 && step == 1
                        TimeArray[1] = 0.0  # Reset time to 0; especially important for initialization in GHX model
                    else
                        TimeArray[1] += Delt
                    end
                    TimeArray[4] = Delt
                    INFO[7] = 0
                    INFO[8] += 1
                    Tout_GHX = OUT[1]  # [C] This is always OUT[1] for ANY GHX model, but if that changes need to update
                    EWT_F = Tout_GHX * 1.8 + 32.0  # [F]
                    # Find COP_Heat (Col 2) and COP_Cool (Col 3) based on temperature (Col 1) map; use interpolation between data points
                    if EWT_F <= p.HeatPumpCOPMap[1, 1]
                        COP_Heat = p.HeatPumpCOPMap[1, 2]
                        COP_Cool = p.HeatPumpCOPMap[1, 3]
                    elseif EWT_F > p.HeatPumpCOPMap[end, 1]
                        COP_Heat = p.HeatPumpCOPMap[end, 2]
                        COP_Cool = p.HeatPumpCOPMap[end, 3]
                    else
                        # Loop over all temperature points in p.HeatPumpCOPMap, but break out of loop if it is found
                        for (index, tmp) in enumerate(p.HeatPumpCOPMap[2:end, 1])  # Omit first and last temp checks, done above
                            if EWT_F > p.HeatPumpCOPMap[index, 1] && EWT_F <= tmp  # "index" starts at 1, not 2, which is i-1 of tmp index in HeatPumpCOPMap
                                slope_heat = (p.HeatPumpCOPMap[index+1, 2] - p.HeatPumpCOPMap[index, 2]) / (tmp - p.HeatPumpCOPMap[index, 1])
                                COP_Heat = p.HeatPumpCOPMap[index, 2] + (tmp - EWT_F) * slope_heat
                                slope_cool = (p.HeatPumpCOPMap[index+1, 3] - p.HeatPumpCOPMap[index, 3]) / (tmp - p.HeatPumpCOPMap[index, 1])
                                COP_Cool = p.HeatPumpCOPMap[index, 3] + (tmp - EWT_F) * slope_cool
                                break
                            end
                        end
                    end

                    # Override COP_Heat/COP_Cool calculation if Q_Heat/Q_Cool is negative
                    if Q_Heat <= 0.0
                        COP_Heat = 3.0
                    end
                    if Q_Cool <= 0.0
                        COP_Cool = 3.0
                    end

                    # Calculate the heat pump outputs
                    P_WSHP_H = Q_Heat / COP_Heat
                    P_WSHP_C = Q_Cool / COP_Cool

                    Q_Absorbed = Q_Heat - P_WSHP_H
                    Q_Rejected = Q_Cool + P_WSHP_C

                    Tons_HeatPump_H = Q_Heat / 1.055 / 12000.0
                    Tons_HeatPump_C = Q_Cool / 1.055 / 12000.0

                    if (Tons_HeatPump_H + Tons_HeatPump_C) > 0.0
                        GPM_GHX = max(p.fMin_VSP_GHXPump * p.GPM_GHXPump,
                                        (Tons_HeatPump_H + Tons_HeatPump_C) * p.GPMperTon_WSHP)
                    else
                        GPM_GHX = 0.0
                    end

                    Mdot_GHX = GPM_GHX * 60.0 / 264.172 * p.Rho_GHXFluid

                    if Mdot_GHX > 0.0
                        Tout_HeatPumps_GHX = Tout_GHX + (Q_Rejected - Q_Absorbed) / Mdot_GHX / p.Cp_GHXFluid
                    else
                        Tout_HeatPumps_GHX = Tout_GHX
                    end

                    # Calculate the pump outputs
                    if GPM_GHX > 0.0
                        f_GHXPump = GPM_GHX / p.GPM_GHXPump
                    else
                        f_GHXPump = 0.0
                    end

                    P_GHXPump = p.Prated_GHXPump * (f_GHXPump^p.Exponent_GHXPump)
                    Qf_GHXPump = P_GHXPump * 0.3
                    if Mdot_GHX > 0.0
                        Tout_Pump = Tout_HeatPumps_GHX + Qf_GHXPump / Mdot_GHX / p.Cp_GHXFluid
                    else
                        Tout_Pump = Tout_HeatPumps_GHX
                    end

                    # Call the GHX model each iteration
                    ErrorFound = [false]
                    if p.ghx_model == "DST"
                        XIN[1] = Tout_Pump
                        XIN[2] = Mdot_GHX
                        XIN[3] = p.AmbientTemperature[hr]
                        XIN[4] = p.AmbientTemperature[hr]
                        XIN[5] = 1.0
                        # Call the ground model
                        ccall((:type557_, "/opt/julia_src/dst.so"), Cvoid,
                        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
                        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)
                    elseif p.ghx_model == "TESS"
                        XIN[1] = Tout_Pump
                        XIN[2] = Mdot_GHX
                        XIN[3] = p.AmbientTemperature[hr]
                        XIN[4] = p.AmbientTemperature[hr]
                        ccall((:type1373_, "/opt/julia_src/GHX-Model-REopt/tess.so"), Cvoid,
                        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
                        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)
                    end

                    # Calculate the outputs
                    # This might be an inefficient calc to do every timestep (could just do at the end if not using before)
                    r.Tmax_GHX = max(r.Tmax_GHX, OUT[1])
                    r.Tmin_GHX = min(r.Tmin_GHX, OUT[1])

                    r.Power_GHXPump += P_GHXPump * Delt
                    r.Qfluid_GHXPump += Qf_GHXPump * Delt
                    r.Power_WSHP_C += P_WSHP_C * Delt
                    r.Power_WSHP_H += P_WSHP_H * Delt
                    r.LoadMet_WSHP_C += Q_Cool * Delt
                    r.LoadMet_WSHP_H += Q_Heat * Delt
                    r.Q_Rejected_Total += Q_Rejected * Delt
                    r.Q_Absorbed_Total += Q_Absorbed * Delt
                    r.Q_GHX_Net += (Q_Rejected - Q_Absorbed) * Delt

                    if p.ghx_model == "DST"
                        r.Q_GHX_In += OUT[4] * Delt
                        r.Q_GHX_Top += OUT[5] * Delt
                        r.Q_GHX_Sides += OUT[6] * Delt
                        r.Q_GHX_Bottom += OUT[7] * Delt
                        r.Q_GHX_Stored += OUT[8] * Delt
                    elseif p.ghx_model == "TESS"
                        r.Q_GHX_In += OUT[3] * Delt
                        r.Q_GHX_Top += OUT[5] * Delt
                        r.Q_GHX_Sides += OUT[8] * Delt
                        r.Q_GHX_Bottom += OUT[10] * Delt
                        r.Q_GHX_Stored += OUT[3] * Delt
                        r.Qf_GHX_Stored += OUT[12] * Delt
                    end

                    # EWT is an AVERAGE temperature over the hour
                    r.EWT[8760*(year-1)+hr] += OUT[1]*Delt

                    r.P_GHXPump_Hourly[8760*(year-1)+hr] += P_GHXPump * Delt / 3600.0
                    r.P_WSHPc_Hourly[8760*(year-1)+hr] += P_WSHP_C * Delt / 3600.0
                    r.P_WSHPh_Hourly[8760*(year-1)+hr] += P_WSHP_H * Delt / 3600.0
                    r.Qc_Hourly[8760*(year-1)+hr] += Q_Cool * Delt / 3600.0
                    r.Qh_Hourly[8760*(year-1)+hr] += Q_Heat * Delt / 3600.0

                    # Call the ground heat exchanger model to clean up as the timestep is complete (assign Ti=Tf)
                    INFO[13] = 1
                    if p.ghx_model == "DST"
                        ccall((:type557_, "/opt/julia_src/dst.so"), Cvoid,
                        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
                        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)
                    elseif p.ghx_model == "TESS"
                        ccall((:type1373_, "/opt/julia_src/GHX-Model-REopt/tess.so"), Cvoid,
                        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
                        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)
                    end
                    INFO[13] = 0
                end
            end
            println("EFT after year ", year ," = ", round(OUT[1] * 1.8 + 32.0, digits=1)," F")
        end
        # Use a Newton's method to calculate the bore size to meet the design goals
        # FX => error of temperature limits
        r.N_Bores_Final = r.N_Bores[size_iter]
        FX_Cooling = p.Tmax_Sizing - r.Tmax_GHX
        FX_Heating = r.Tmin_GHX - p.Tmin_Sizing
        # If tolerance is met for either heating or cooling, break the while loop and accept size solution
        if (abs(FX_Cooling) < p.solver_eft_tolerance) && (abs(FX_Heating) < p.solver_eft_tolerance)
            r.FX_Now[size_iter] = min(FX_Cooling, FX_Heating)
            break
        elseif (abs(FX_Cooling) < p.solver_eft_tolerance) && (FX_Heating > 0.0)
            r.FX_Now[size_iter] = FX_Cooling
            break
        elseif (abs(FX_Heating) < p.solver_eft_tolerance) && (FX_Cooling > 0.0)
            r.FX_Now[size_iter] = FX_Heating
            break
        end
        # Check the first iteration if it solved and if not, give it a next best guess
        if size_iter == 1
            # X_Now[size_iter=1] is already defined by X_init
            # FX_Now[size_iter=1] is NOT already defined (removed "X/FX_Previous" from TESS)
            if (FX_Cooling > 0.0) && (FX_Heating > 0.0)
                r.X_Now[size_iter+1] = r.X_Now[size_iter] / 2.0
                r.FX_Now[size_iter] = max(FX_Cooling, FX_Heating)
            elseif FX_Cooling > 0.0
                r.X_Now[size_iter+1] = 2.0 * r.X_Now[size_iter]
                r.FX_Now[size_iter] = FX_Heating
            elseif FX_Heating > 0.0
                r.X_Now[size_iter+1] = 2.0 * r.X_Now[size_iter]
                r.FX_Now[size_iter] = FX_Cooling
            else
                r.X_Now[size_iter+1] = 2.0 * r.X_Now[size_iter]
                r.FX_Now[size_iter] = min(FX_Cooling, FX_Heating)
            end
            size_iter += 1
            continue
        elseif size_iter > 1
            if (FX_Cooling > p.solver_eft_tolerance) || (FX_Heating > p.solver_eft_tolerance)
                if FX_Heating < p.solver_eft_tolerance
                    r.FX_Now[size_iter] = FX_Heating
                elseif FX_Cooling < p.solver_eft_tolerance
                    r.FX_Now[size_iter] = FX_Cooling
                else  # Both are greater than tolerance
                    r.FX_Now[size_iter] = min(FX_Cooling, FX_Heating)
                end
            else
                r.FX_Now[size_iter] = max(FX_Cooling, FX_Heating)
            end
            X_New = r.X_Now[size_iter] - r.FX_Now[size_iter] *
                    (r.X_Now[size_iter] - r.X_Now[size_iter-1]) /
                    (r.FX_Now[size_iter] - r.FX_Now[size_iter-1])
            if X_New <= 0.0
                X_New = min(r.X_Now[size_iter], r.X_Now[size_iter-1]) / 2.0
            end
            X_New = min(r.X_Now[size_iter] * 2.0, X_New)
            X_New = max(r.X_Now[size_iter] / 2.0, X_New)
            r.X_Now[size_iter+1] = X_New

            #!Check to see if the program is stuck in a tolerance situation hunting a solution between two adjacent # bore sizes
            if abs(r.N_Bores[size_iter] - r.N_Bores[size_iter-1]) == 1
                if (r.FX_Now[size_iter] < 0) && (r.FX_Now[size_iter-1] > 0)
                    break
                elseif (r.FX_Now[size_iter] > 0) && (r.FX_Now[size_iter-1] < 0)
                    break
                end
            end

            size_iter += 1
            continue
        end
    end
    # Record actual final EWT error; if the loop exits by hitting max_sizing_iterations, record previous error
    if r.FX_Now[size_iter] == 0
        r.FX_Final = r.FX_Now[size_iter-1]
    else
        r.FX_Final = r.FX_Now[size_iter]
    end
    return r
end

"""
    assign_PAR!(p, r, PAR, size_iter)

Assign the `PAR` parameter array for the GHX compiled Fortran subroutine.

"""
function assign_PAR!(p, r, PAR, size_iter)
    if p.ghx_model == "DST"
        PAR[1] = r.StorageVolume
        PAR[2] = r.Length_Boreholes + p.Depth_Header
        PAR[3] = p.Depth_Header
        PAR[4] = r.N_Bores[size_iter]
        PAR[5] = p.Radius_Bores
        PAR[6] = p.N_Series
        PAR[7] = p.N_Radial
        PAR[8] = p.N_Vertical
        PAR[9] = p.K_Soil
        PAR[10] = p.RhoCp_Soil
        PAR[11] = -1.0
        PAR[12] = p.Ro_Pipe
        PAR[13] = p.Ri_Pipe
        PAR[14] = p.Center_Center_Distance / 2.0  # Was Center_Center_HalfDistance
        PAR[15] = p.K_Grout
        PAR[16] = p.K_Pipe
        PAR[17] = p.K_Grout
        PAR[18] = 0.0
        PAR[19] = r.Mdot_Borehole
        PAR[20] = p.T_Ground
        PAR[21] = -1.0
        PAR[22] = p.Cp_GHXFluid
        PAR[23] = p.Rho_GHXFluid
        PAR[24] = 0.0
        PAR[25] = 0.0
        PAR[26] = 0.0
        PAR[27] = 0.04 * 3.6
        PAR[28] = 25.0
        PAR[29] = 100.0
        PAR[30] = p.T_Ground
        PAR[31] = 0.0
        PAR[32] = 0.0
        PAR[33] = p.T_Ground
        PAR[34] = p.T_Ground
        PAR[35] = 0.
        PAR[36] = p.T_Ground
        PAR[37] = p.Tamp_Ground
        PAR[38] = p.DayMin_DST
        PAR[39] = 1.0
        PAR[40] = p.K_Soil
        PAR[41] = p.RhoCp_Soil
        PAR[42] = r.Depth_DST
        PAR[43] = 0.0
        PAR[44] = 0.0
    elseif p.ghx_model == "TESS"
        PAR[1] = r.N_Bores[size_iter]
        PAR[2] = 1.0
        PAR[3] = r.Length_Boreholes
        PAR[4] = p.BoreSpacing
        PAR[5] = p.SpacingType
        PAR[6] = p.Radius_Bores * 2.0
        PAR[7] = p.K_Grout
        PAR[8] = p.Center_Center_Distance
        PAR[9] = p.Ri_Pipe * 2.0
        PAR[10] = p.Ro_Pipe * 2.0
        PAR[11] = p.K_Pipe
        PAR[12] = p.Depth_Header
        PAR[13] = p.Cp_GHXFluid
        PAR[14] = p.Rho_GHXFluid
        PAR[15] = p.K_GHXFluid
        PAR[16] = p.Mu_GHXFluid
        PAR[17] = 0.0
        PAR[18] = 0.0
        PAR[19] = 0.0
        PAR[20] = p.K_Soil
        PAR[21] = p.Rho_Soil
        PAR[22] = p.Cp_Soil
        PAR[23] = p.T_Ground
        PAR[24] = p.Tamp_Ground
        PAR[25] = p.DayMin_Surface
        PAR[26] = 0.0
        PAR[27] = r.Radius_Field  #!Far-field distance
        PAR[28] = p.Depth_Bores / 2.0  #!Deep earth distance
        PAR[29] = 0.1
        PAR[30] = 2.0
        PAR[31] = 0.02
        PAR[32] = 2.0
    end
end

"""
    init_ghx_calls_2x!(p, TimeArray, XIN, OUT, PAR, INFO)

Run the initial two (2) GHX Fortran subroutines using initialization parameters.

"""
function init_ghx_calls_2x!(p, TimeArray, XIN, OUT, PAR, INFO)
    # Call the GHX model two times for initialization
    INFO[1] = 1
    INFO[7] = -1
    INFO[8] = 2
    ErrorFound = [false]
    if p.ghx_model == "DST"
        INFO[2] = 557
        XIN[1] = p.T_Ground
        XIN[2] = 0.0
        XIN[3] = 20.0
        XIN[4] = 20.0
        XIN[5] = 1.0

        # First time
        ccall((:type557_, "/opt/julia_src/dst.so"), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)

        INFO[7] = 0
        INFO[8] = 3

        # Second time
        ccall((:type557_, "/opt/julia_src/dst.so"), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)
    elseif p.ghx_model == "TESS"
        INFO[2] = 1373
        XIN[1] = p.T_Ground
        XIN[2] = 0.0
        XIN[3] = 20.0
        XIN[4] = 20.0

        # First time
        ccall((:type1373_, "/opt/julia_src/GHX-Model-REopt/tess.so"), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)

        INFO[7] = 0
        INFO[8] = 3

        # Second time
        ccall((:type1373_, "/opt/julia_src/GHX-Model-REopt/tess.so"), Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}),
        TimeArray, XIN, OUT, PAR, INFO, ErrorFound)
    end
end

# -------
end # module GHPGHX


text:
"""
ParametersOrthophoto

# Constructors
```jldoctest
ParametersOrthophoto(
	txtpotreedirs::String,
	outputimage::String,
	bbin::Union{String,AABB},
	GSD::Float64,
	PO::String,
	altitude::Union{Float64,Nothing},
	thickness::Union{Float64,Nothing},
	ucs::Union{String,Matrix{Float64}},
	BGcolor::Array{Float64,1},
	pc::Bool
	)
```
Input description:
 - txtpotreedirs: path to Potree project or a text file containing a list of Potree projects
 - outputimage: Output filename for image
 - bbin: Path to json file or AABB
 - GSD: Ground sampling distance
 - PO: Projection plane and sight direction. Option: "XY+","XY-","XZ+","XZ-","YZ+","YZ-"
 - quote: Distance of plane from the origin
 - thickness: Thickness of plane
 - ucs: User coordinates system
 - BGcolor: Background color image
 - pc: If true return LAS file of point cloud

# Fields
```jldoctest
PO::String
outputimage::String
outputfile::String
potreedirs::Array{String,1}
model::LAR
coordsystemmatrix::Array{Float64,2}
RGBtensor::Array{Float64,3}
rasterquote::Array{Float64,2}
GSD::Float64
refX::Float64
refY::Float64
pc::Bool
ucs::Matrix
tightBB::AABB
mainHeader::LasIO.LasHeader
numPointsProcessed::Int64
numNodes::Int64
numFilesProcessed::Int64
stream_tmp::Union{Nothing,IOStream}
```

Fields description:
 - PO: Projection plane and sight direction. Option: "XY+","XY-","XZ+","XZ-","YZ+","YZ-"
 - outputimage: Output filename for image
 - outputfile: Output filename for point cloud
 - potreedirs: Potree projects
 - model: Cuboidal LAR model
 - coordsystemmatrix: Rotation matrix related to PO
 - RGBtensor: Image tensor
 - rasterquote: Matrix of heigth (z-buffer)
 - GSD: Ground sampling distance
 - refX: x coordinate of landmark
 - refY: y coordinate of landmark
 - pc: if true extract point cloud
 - ucs: user coordinate system
 - tightBB: segmented point cloud AABB
 - mainHeader: header of point cloud file las
 - numPointsProcessed: number of points processed
 - numNodes: number of files in potree
 - numFilesProcessed: number of files processed
 - stream_tmp: temporary files for laspoints

"""
mutable struct ParametersOrthophoto
    PO::String
    outputimage::String
    outputfile::String
    potreedirs::Array{String,1}
    model::LAR
    coordsystemmatrix::Array{Float64,2}
    RGBtensor::Array{Float64,3}
    rasterquote::Array{Float64,2}
    GSD::Float64
    refX::Float64
    refY::Float64
    pc::Bool
    ucs::Matrix
    tightBB::AABB
    mainHeader::LasIO.LasHeader
    numPointsProcessed::Int64
    numNodes::Int64
    numFilesProcessed::Int64
    stream_tmp::Union{Nothing,IOStream}

    function ParametersOrthophoto(
        txtpotreedirs::String,
        outputimage::String,
        bbin::Union{String,AABB},
        GSD::Float64,
        PO::String,
        altitude::Union{Float64,Nothing},
        thickness::Union{Float64,Nothing},
        ucs::Union{String,Matrix{Float64}},
        BGcolor::Array{Float64,1},
        pc::Bool,
        epsg::Union{Nothing,Integer},
    )

        # check validity
        @assert length(PO) == 3 "orthoprojectionimage: $PO not valid view"

        numPointsProcessed = 0
        numNodes = 0
        numFilesProcessed = 0
        stream_tmp = nothing

        outputfile = splitext(outputimage)[1] * ".las"

        potreedirs = get_potree_dirs(txtpotreedirs)

        if typeof(ucs) == Matrix{Float64}
            coordsystemmatrix = PO2matrix(PO, ucs)
        else
            ucs = FileManager.ucs2matrix(ucs)
            coordsystemmatrix = PO2matrix(PO, ucs)
        end

        model = getmodel(bbin)
        aabb = AABB(model[1])

        if !isnothing(altitude) && !isnothing(thickness)
            directionview = PO[3]
            if directionview == '-'
                altitude = -altitude
            end
            origin = Common.inv(ucs)[1:3, 4]
            model = getmodel(
                Common.inv(coordsystemmatrix),
                altitude,
                thickness,
                aabb;
                new_origin = origin,
            )
            aabb = AABB(model[1])
        end

        RGBtensor, rasterquote, refX, refY =
            init_raster_array(coordsystemmatrix, GSD, model, BGcolor)

        #per l'header devo creare il nuovo AABB dato dal nuovo orientamento.
        new_verts_BB = Common.apply_matrix(ucs, model[1])
        aabb = AABB(new_verts_BB)
        mainHeader = FileManager.newHeader(
            aabb,
            "ORTHOPHOTO",
            FileManager.SIZE_DATARECORD,
        )

        if !isnothing(epsg)
            FileManager.LasIO.epsg_code!(mainHeader, epsg)
        end

        return new(
            PO,
            outputimage,
            outputfile,
            potreedirs,
            model,
            coordsystemmatrix,
            RGBtensor,
            rasterquote,
            GSD,
            refX,
            refY,
            pc,
            ucs,
            AABB(-Inf, Inf, -Inf, Inf, -Inf, Inf),
            mainHeader,
            numPointsProcessed,
            numNodes,
            numFilesProcessed,
            stream_tmp,
        )
    end

end

"""
	Point

Point information.

# Constructors
```jldoctest
Point(lasPoint::FileManager.LasPoint, header::FileManager.LasHeader)::Point
```

# Fields
```jldoctest
position::Vector{Float64}
color::Vector{UInt16}
normal::Vector{Float64}
intensity::UInt16
classification::Char
returnNumber::Char
numberOfReturns::Char
pointSourceID::UInt16
gpsTime::Float64
```
"""
mutable struct Point
    position::Vector{Float64}
    color::Vector{LasIO.FixedPointNumbers.N0f16}
    normal::Vector{Float64}
    intensity::UInt16
    classification::Char
    returnNumber::Char
    numberOfReturns::Char
    pointSourceID::UInt16
    gpsTime::Float64

    function Point(
        lasPoint::FileManager.LasPoint,
        header::FileManager.LasHeader,
    )::Point
        position = FileManager.xyz(lasPoint, header)
        color = [lasPoint.red, lasPoint.green, lasPoint.blue]
        normal = Float64[]
        intensity = lasPoint.intensity
        classification = lasPoint.raw_classification
        returnNumber = 0
        numberOfReturns = 0
        pointSourceID = lasPoint.pt_src_id
        gpsTime = 0.0
        return new(
            position,
            color,
            normal,
            intensity,
            classification,
            returnNumber,
            numberOfReturns,
            pointSourceID,
            gpsTime,
        )
    end

end

function Base.show(io::IO, point::Point)
    println(io, "position: $(point.position)")
    println(io, "color: $(point.color)")
end

ext: jl


text:
using Pandas
using IteratorInterfaceExtensions
using TableTraits
using DataValues
using Test

@testset "TableTraits" begin

table_array = [(a=1, b="John", c=3.2), (a=2, b="Sally", c=5.8)]

df = DataFrame(table_array)

@test collect(columns(df)) == ["a", "b", "c"]
@test values(df[:a]) == [1,2]
@test values(df[:c]) == [3.2, 5.8]
@test [df[:b][i] for i in 1:2] == ["John", "Sally"]

@test TableTraits.isiterabletable(df) == true

it = IteratorInterfaceExtensions.getiterator(df)

@test eltype(it) == NamedTuple{(:a,:b,:c),Tuple{Int,String,Float64}}

it_collected = collect(it)

@test eltype(it_collected) == NamedTuple{(:a,:b,:c),Tuple{Int,String,Float64}}
@test length(it_collected) == 2
@test it_collected[1] == (a=1, b="John", c=3.2)
@test it_collected[2] == (a=2, b="Sally", c=5.8)

@test TableTraits.supports_get_columns_copy_using_missing(df) == true
cols = TableTraits.get_columns_copy_using_missing(df)
@test cols == (a=[1,2], b=["John", "Sally"], c=[3.2, 5.8])

table_array2 = [(a=1, b=DataValue("John"), c=3.2), (a=2, b=DataValue("Sally"), c=5.8)]

@test_throws ArgumentError DataFrame(table_array2)

table_array3 = [(a=DataValue{Int}(), b="John", c=DataValue(3.2)), (a=DataValue(2), b="Sally", c=DataValue{Float64}())]

df3 = DataFrame(table_array3)

it3_collected = collect(IteratorInterfaceExtensions.getiterator(df3))

@test length(it3_collected) == 2
@test isnan(it3_collected[1].a)
@test it3_collected[1].b == "John"
@test it3_collected[1].c == 3.2
@test it3_collected[2].a == 2
@test it3_collected[2].b == "Sally"
@test isnan(it3_collected[2].c)

cols3 = TableTraits.get_columns_copy_using_missing(df3)
@test isequal(cols3, (a=[NaN,2.], b=["John", "Sally"], c=[3.2, NaN]))

end


text:
using LightGraphs
using LightGraphs.SimpleGraphs
using LinearAlgebra

function normalized_laplacian(g::SimpleGraph)
    L = float(laplacian_matrix(g))
    for i in vertices(g)
        if L[i,i] != 0
            L[i,:] /= L[i,i]
        end
    end
    return L
end

function normalized_laplacian(g::SimpleDiGraph)
    L = float(laplacian_matrix(g))
    for i in vertices(g)
        if L[i,i] != 0
            L[i,:] /= L[i,i]
        end
    end
    return L
end

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
module Compologicircuits

#include("Circuits.jl")
include("AugmentedCircuit.jl")
#include("SAT.jl")
#include("Clc.jl")

export CircuitDom, Circuit, show_diagram, apply, iNOT, iAND, iOR, Impl, Circuits,
Augmented_Circuit, aiAND, aiNOT, aiOR, AImpl

# Make some circuits

#=XOR = @program Circuits (x::B, y::B) begin
    xnandy = NOT(AND(x,y))
    xory = OR(x,y)
    return AND(xnandy,xory)
end

XOR_expr = to_hom_expr(FreeCartesianCategory, XOR)

iXOR = Impl(XOR_expr)

# Print XOR truth table

print_xor() = for i in Iterators.product([false, true], [false, true])
    input = collect(i)
    output = apply(iXOR, input)
    println("XOR on $input = $output")
end

full_adder = @program Circuits (x::B, y::B, cin::B) begin
    z = XOR(x,y)
    s = XOR(z,cin)
    cout = OR(AND(z,cin),AND(x,y))
    return (s, cout)
end

AND = to_wiring_diagram(Circuits[:AND])
NOT = to_wiring_diagram(Circuits[:NOT])
OR = to_wiring_diagram(Circuits[:OR])
total_full_adder = ocompose(full_adder, [XOR, XOR, AND, AND, OR])
iFA = Impl(to_hom_expr(FreeCartesianCategory, total_full_adder))

# Print Full Adder truth table

print_fa() = for i in Iterators.product([false, true], [false, true], [false, true])
    input = collect(i)
    output = apply(iFA, input)
    println("FA on $input = $output")
end

four_bit_adder = @program Circuits (cin::B, x1::B, x2::B, x3::B, x4::B, y1::B, y2::B, y3::B, y4::B) begin
    s1, cout1 = FA(x4,y4,cin)
    s2, cout2 = FA(x3,y3,cout1)
    s3, cout3 = FA(x2,y2,cout2)
    s4, cout = FA(x1,y1,cout3)
    return (cout, s4, s3, s2, s1)
end

total_four_bit_adder = ocompose(four_bit_adder, [total_full_adder, total_full_adder, total_full_adder, total_full_adder])
tfba_expr = to_hom_expr(FreeCartesianCategory, total_four_bit_adder)
i4bitadder = Impl(tfba_expr)

sum = apply(i4bitadder,
    Bool[0,      #carry in
        0,1,1,1, #x
        0,1,1,0] #y
)=#

end


text:
function test_tmatrix_deriv(scenario)
    @unpack ea, tx, bfield, species = scenario

    M = LMP.susceptibility(LMPParams().topheight, tx.frequency, bfield, species)

    for i = 1:4
        for j = 1:4
            Tfcn(θ) = (ea = EigenAngle(θ); T = LMP.tmatrix(ea, M)[i,j])
            dTref = FiniteDiff.finite_difference_derivative(Tfcn, θs, Val{:central})
            dT(θ) = (ea = EigenAngle(θ); T = LMP.dtmatrix(ea, M)[i,j])

            @test maxabsdiff(dT.(θs), dTref) < 1e-7
        end
    end
end

@testset "TMatrix.jl" begin
    @info "Testing TMatrix"

    for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario,
        multiplespecies_scenario)
        test_tmatrix_deriv(scn)
    end
end


text:
using FTCTests
const FTC = FaultTolerantControl
using JLD2
using Transducers
using UnPack
using Flux
using Flux.Data: DataLoader
using LinearAlgebra

function faults_to_effectiveness(faults; dim_input=6, _fault_time=0.0)
    # actuator faults
    actuator_faults = faults |> Filter(fault -> typeof(fault) <: AbstractActuatorFault) |> collect
    actuator_fault_times = actuator_faults |> Map(fault -> fault.time) |> collect
    @assert actuator_fault_times == _fault_time*ones(size(actuator_fault_times)...)  # fault time: assumed to be zero
    affect_actuator! = FTC.Affect!(AbstractActuatorFault[actuator_faults...])
    # effectiveness matrix
    Λ = ones(dim_input) |> diagm
    affect_actuator!(Λ, _fault_time)
    λ = diag(Λ)
end

function preprocess(file_path::String; cf=PositionAngularVelocityCostFunctional(), verbose=true)
    res = try
        jldopen(file_path, "r") do jld2
            df = jld2["df"]
            traj_des = jld2["traj_des"]
            faults = jld2["faults"]
            @assert length(faults) == 1  # currently, only single fault is considered
            # desired trajectory parameter
            _θ = vcat(traj_des.θ...)  # concatenated control points
            # actuator fault (effectiveness)
            λ = faults_to_effectiveness(faults)
            # initial condition
            x0 = df.sol[1].plant.state
            # cost
            ts = df.time
            poss = df.sol |> Map(datum -> datum.plant.state.p) |> collect
            poss_des = ts |> Map(traj_des) |> collect
            e_ps = poss .- poss_des
            if verbose
                @warn("TODO: there is no desired angular velocity")
            end
            e_ωs = ts |> Map(t -> zeros(3)) |> collect
            J = cost(cf, ts, e_ps, e_ωs)
            (; x0=x0, λ=λ, _θ=_θ, J=J)
        end
    catch
        @warn("Fail to load a file; is the file loadable via JLD2?")
        missing
    end
end

function preprocess(file_paths::Vector{String}; cf=PositionAngularVelocityCostFunctional(), verbose=true)
    data = file_paths |> Map(path -> preprocess(path; cf=cf, verbose=verbose)) |> collect
    data_filtered = filter(x -> typeof(x) != Missing, data)
end

function training_test(Ĵ, data_train, data_test, epochs)
    _data_train = make_a_trainable(data_train)
    _data_test = make_a_trainable(data_test)
    loss(d) = Flux.Losses.mse(Ĵ(d.feature), d.J)
    opt = ADAM(1e-3)
    ps = Flux.params(Ĵ)
    dataloader = DataLoader(_data_train; batchsize=32, shuffle=true, partial=false)
    println("Training $(epochs) epoch...")
    for epoch in 0:epochs
        println("epoch: $(epoch)/$(epochs)")
        if epoch != 0
            for d in dataloader
                train_loss, back = Flux.Zygote.pullback(() -> loss(d), ps)
                gs = back(one(train_loss))
                Flux.update!(opt, ps, gs)
            end
        end
        println("train loss: $(loss(_data_train)), test loss: $(loss(_data_test))")
    end
end

function make_a_trainable(data)
    # feature
    features = data |> Map(datum -> vcat(
                                         datum.x0,
                                         datum.λ,
                                         datum._θ,
                                        )) |> collect
    # output
    Js = data |> Map(datum -> datum.J) |> collect
    _data = (;
             feature = hcat(features...),
             J = hcat(Js...),
            )  # for Flux
end

function main(epochs; dir_log="data/hovering/adaptive", seed=2021)
    Random.seed!(seed)
    # load data
    file_paths = readdir(dir_log; join=true)
    @time data = preprocess(file_paths; verbose=false)
    # construct approximator
    n_nt = (;
            x=18,
            λ=6,
            θ=3,
           )  # dimensions
    n = sum(n_nt)  # total dim (feature dimension)
    n_h = 256  # hidden layer nodes
    Ĵ = Chain(
              Dense(n, n_h, leakyrelu),
              Dense(n_h, n_h, leakyrelu),
              Dense(n_h, 1, leakyrelu),
             )
    data_train, data_test = partitionTrainTest(data, 0.9)  # 90:10
    training_test(Ĵ, data_train, data_test, epochs)
end


text:
using ProblemParser
using Utils

using Chain
import IterTools

struct Off0Matrix
    rules::Vector{Bool}
    parent::Matrix{Bool}
    at_infinity::Bool
end
Off0Matrix(rules, matrix) = Off0Matrix(rules, matrix, false)

const file = find_input(@__FILE__)
tobool = Map(e -> e == '#')
const input = Off0Matrix(parse(FirstRest(Blocks(), tobool, Rectangular(tobool)), slurp(file))...)

Base.size(matrix::Off0Matrix) = size(matrix.parent)

Base.checkbounds(::Type{Bool}, matrix::Off0Matrix, ci::CartesianIndex{2}) = checkbounds(Bool, matrix.parent, ci)

function Base.getindex(matrix::Off0Matrix, inds)
    checkbounds(Bool, matrix, inds) ? getindex(matrix.parent, inds) : matrix.at_infinity
end

numfrom3x3(matrix::Off0Matrix, num::Int) = matrix.rules[num + 1]
numfrom3x3(matrix::Off0Matrix, bits::AbstractArray{Bool}) = numfrom3x3(matrix, foldl((l,r) -> (l << 1) | Int(r), bits; init=0))

neighbors = CartesianIndices((-1:1, -1:1)) |> permutedims

function Base.iterate(matrix::Off0Matrix, state=matrix)
    border = 1
    output = Matrix{Bool}(undef, (size(state) .+ 2border)...)
    for ci in CartesianIndices(output)
        output[ci] = @chain ci begin
            [state[i - CartesianIndex(1,1)] for i in (neighbors .+ [_])]
            numfrom3x3(state, _)
        end
    end
    # Also calculate new value at infinity.
    output, Off0Matrix(state.rules, output, numfrom3x3(state, Iterators.repeated(state.at_infinity, 9) |> collect))
end

Base.IteratorSize(::Off0Matrix) = Base.IsInfinite()

solve(input, n) = IterTools.nth(input, n) |> count

part1 = Base.Fix2(solve, 2)

assertequal(part1(input), 5400)

part2 = Base.Fix2(solve, 50)

assertequal(part2(input), 18989)

max_stars_repo_name: CmdQ/AoC


text:
using Particles
using StatsBase
using StatsPlots

nx = 16
nv = 12

n = 10000
x = randn(n)
v = randn(n)

xmax = ceil(maximum(abs.(x)))
x .+= xmax
x ./= 2*xmax

vmax = ceil(maximum(abs.(v)))

p = PoissonSolver{eltype(x)}(nx)
solve!(p, x)

xgrid = collect(0:p.Δx:1)
vgrid = LinRange(-vmax, +vmax, nv)

hx = fit(Histogram, x, xgrid)
hv = fit(Histogram, v, vgrid)

plot(hx; legend=nothing, xlabel="x", ylabel="n")
savefig("poisson_histogram_x.png")

plot(hv; legend=nothing, xlabel="v", ylabel="n")
savefig("poisson_histogram_v.png")

plot(p.xgrid[1:end-1], p.ρ; legend=nothing, xlabel="x", ylabel="ρ(x)")
savefig("poisson_ρ.png")

plot(p.xgrid[1:end-1], p.ϕ; legend=nothing, xlabel="x", ylabel="ϕ(x)")
savefig("poisson_ϕ.png")

plot(p.xgrid, x -> eval_field(p, x); legend=nothing, xlabel="x", ylabel="E(x)")
savefig("poisson_E.png")

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
function arg_type(arg, ismethod)
    if ismethod
        if hasbinding(arg)
            if bindingof(arg) isa Binding && bindingof(arg).type !== nothing
                type = bindingof(arg).type
                if type isa Binding && type.val isa SymbolServer.DataTypeStore
                    type = type.val
                end
                return type
            end
        end
    else
        if hasref(arg)
            if refof(arg) isa Binding && refof(arg).type !== nothing
                type = refof(arg).type
                if type isa Binding && type.val isa SymbolServer.DataTypeStore
                    type = type.val
                end
                return type
            end
        elseif headof(arg) === :STRING
            return CoreTypes.String
        elseif headof(arg) === :CHAR
            return CoreTypes.Char
        elseif headof(arg) === :FLOAT
            return CoreTypes.Float64
        elseif headof(arg) === :INT
            return CoreTypes.Int
        elseif headof(arg) === :HEXINT
            if length(arg.val) < 5
                return CoreTypes.UInt8
            elseif length(arg.val) < 7
                return CoreTypes.UInt16
            elseif length(arg.val) < 11
                return CoreTypes.UInt32
            else
                return CoreTypes.UInt64
            end
        elseif headof(arg) === :TRUE || headof(arg) === :FALSE
            return CoreTypes.Bool
        elseif isquotedsymbol(arg)
            return SymbolServer.stdlibs[:Core][:Symbol]
        end
    end
    # VarRef(VarRef(nothing, :Core), :Any)
    CoreTypes.Any
end

isquotedsymbol(x) = x isa EXPR && x.head === :quotenode && length(x.args) == 1 && x.args[1].head === :IDENTIFIER && hastrivia(x)

function call_arg_types(call::EXPR, ismethod)
    types, kws = [], []
    call.args === nothing && return types, kws
    if length(call.args) > 1 && headof(call.args[2]) === :parameters
        for i = 1:length(call.args[2].args)
            push!(kws, call.args[2].args[i].args[1])
        end
        for i = 3:length(call.args)
            push!(types, arg_type(call.args[i], ismethod))
        end
    else
        for i = 2:length(call.args)
            push!(types, arg_type(call.args[i], ismethod))
        end
    end
    types, kws
end

function method_arg_types(call::EXPR)
    types, opts, kws = [], [], []
    call.args === nothing && return types, opts, kws
    if length(call.args) > 1 && headof(call.args[2]) === :parameters
        for i = 1:length(call.args[2].args)
            push!(kws, call.args[2].args[i].args[1])
        end
        for i = 3:length(call.args)
            if CSTParser.iskwarg(call.args[i])
                push!(opts, arg_type(call.args[i].args[1], true))
            else
                push!(types, arg_type(call.args[i], true))
            end
        end
    else
        for i = 2:length(call.args)
            if CSTParser.iskwarg(call.args[i])
                push!(opts, arg_type(call.args[i].args[1], true))
            else
                push!(types, arg_type(call.args[i], true))
            end
        end
    end
    types, opts, kws
end

function find_methods(x::EXPR, store)
    possibles = []
    if iscall(x)
        length(x.args) === 0 && return possibles
        func_ref = refof_call_func(x)
        func_ref === nothing && return possibles
        args, kws = call_arg_types(x, false)
        if func_ref isa Binding && func_ref.val isa SymbolServer.FunctionStore ||
            func_ref isa Binding && func_ref.val isa SymbolServer.DataTypeStore
            func_ref = func_ref.val
        end
        if func_ref isa SymbolServer.FunctionStore || func_ref isa SymbolServer.DataTypeStore
            for method in func_ref.methods
                if match_method(args, kws, method, store)
                    push!(possibles, method)
                end
            end
        elseif func_ref isa Binding
            if (CoreTypes.isfunction(func_ref.type) || CoreTypes.isdatatype(func_ref.type)) && func_ref.val isa EXPR
                for method in func_ref.refs
                    method = get_method(method)
                    if method !== nothing
                        if method isa SymbolServer.FunctionStore
                            for method1 in method.methods
                                if match_method(args, kws, method1, store)
                                    push!(possibles, method1)
                                end
                            end
                        elseif match_method(args, kws, method, store)
                            push!(possibles, method)
                        end
                    end
                end
            elseif (method = method_of_callable_datatype(func_ref)) !== nothing
                if match_method(args, kws, method, store)
                    push!(possibles, method)
                end
            end
        end
    end
    possibles
end

function match_method(args::Vector{Any}, kws::Vector{Any}, method::SymbolServer.MethodStore, store)
    !isempty(kws) && isempty(method.kws) && return false
    nmargs = length(method.sig)
    varargval = nothing
    if nmargs > 0 && last(method.sig)[2] isa SymbolServer.FakeTypeofVararg
        if length(args) == nmargs - 1
            nmargs -= 1
            # vararg can be zero length
        elseif length(args) >= nmargs
            # set aside the type param of the Vararg for later use
            varargval = last(method.sig)[2].T
        end
    end
    if length(args) == nmargs
        for i in 1:length(args)
            if varargval !== nothing && i >= nmargs
                !_issubtype(args[i], varargval, store) && !_issubtype(varargval, args[i], store) && return false
            else
                !_issubtype(args[i], method.sig[i][2], store) && !_issubtype(method.sig[i][2], args[i], store) && return false
            end

        end
        return true
    end
    return false
end

function match_method(args::Vector{Any}, kws::Vector{Any}, method::EXPR, store)
    margs, mkws = [], []
    vararg = false
    if CSTParser.defines_struct(method)
        for i in 1:length(method.args[3].args)
            arg = method.args[3].args[i]
            if defines_function(arg)
                # Hit an inner constructor so forget about the default one.
                for arg in method.args[3].args
                    if defines_function(arg)
                        !match_method(args, kws, arg, store) && return false
                    end
                end
                return true
            end
            push!(margs, arg_type(arg, true))
        end
    else
        sig = CSTParser.rem_decl(CSTParser.get_sig(method))
        margs, mopts, mkws = method_arg_types(sig)
        # vararg
        if length(sig.args) > 0
            if CSTParser.issplat(last(sig.args))
                vararg = true
            end
        end
    end
    !isempty(kws) && isempty(mkws) && return false

    if length(margs) < length(args)
        for i in 1:min(length(mopts), length(args) - length(margs))
            push!(margs, mopts[i])
        end
        if vararg
            for _ in 1:length(args) - length(margs)
                push!(margs, CoreTypes.Any)
            end
        end
    end

    if length(args) == length(margs) || (vararg && length(args) == length(margs) - 1)
        for i in 1:length(args)
            !_issubtype(args[i], margs[i], store) && !_issubtype(margs[i], args[i], store) && return false
        end
        return true
    end
    return false
end

function refof_call_func(x)
    if isidentifier(first(x.args)) && hasref(first(x.args))
        return refof(first(x.args))
    elseif is_getfield_w_quotenode(x.args[1]) && (rhs = rhs_of_getfield(x.args[1])) !== nothing && hasref(rhs)
        return refof(rhs)
    else
        return
    end
end

function is_sig_of_method(sig::EXPR, method = maybe_get_parent_fexpr(sig, defines_function))
    method !== nothing && sig == CSTParser.get_sig(method)
end

function method_of_callable_datatype(b::Binding)
    if b.type isa Binding && b.type.type === CoreTypes.DataType
        for ref in b.type.refs
            if ref isa EXPR && ref.parent isa EXPR && isdeclaration(ref.parent) && is_in_fexpr(ref.parent, x -> x.parent isa EXPR && x.parent.head === :call && x == x.parent.args[1] && is_in_funcdef(x.parent))
                return get_parent_fexpr(ref, defines_function)
            end
        end
    end
end


text:
function to_roman(number::Integer)

end


text:
# Test of Dijkstra's algorithm for shorest paths

using Graphs
using Base.Test

# g1: the example in CLRS (2nd Ed.)
g1 = simple_inclist(5)

g1_wedges = [
    (1, 2, 10.),
    (1, 3, 5.),
    (2, 3, 2.),
    (3, 2, 3.),
    (2, 4, 1.),
    (3, 5, 2.),
    (4, 5, 4.),
    (5, 4, 6.),
    (5, 1, 7.),
    (3, 4, 9.) ]

ne = length(g1_wedges)
eweights1 = zeros(ne)
for i = 1 : ne
    we = g1_wedges[i]
    add_edge!(g1, we[1], we[2])
    eweights1[i] = we[3]
end

@assert num_vertices(g1) == 5
@assert num_edges(g1) == 10

s1 = dijkstra_shortest_paths(g1, eweights1, 1)

@test s1.parents == [1, 3, 1, 2, 3]
@test s1.dists == [0., 8., 5., 9., 7.]
@test s1.colormap == [2, 2, 2, 2, 2]

s1 = dijkstra_shortest_paths(g1, eweights1, [1])

@test s1.parents == [1, 3, 1, 2, 3]
@test s1.dists == [0., 8., 5., 9., 7.]
@test s1.colormap == [2, 2, 2, 2, 2]

g1ex = graph([i for i in 1:5], ExEdge{Int}[], is_directed=true)

#todo -- this should be doable in a comprehension, I can't get the types to work
for (i,v) in enumerate(g1_wedges)
    ed = ExEdge(i, v[1], v[2])
    ed.attributes["length"] = v[3]
    add_edge!(g1ex, ed)
end

edgel = AttributeEdgePropertyInspector{Float64}("length")
s1ex = dijkstra_shortest_paths(g1ex, edgel, [1])

@test s1ex.parents == [1, 3, 1, 2, 3]
@test s1ex.dists == [0., 8., 5., 9., 7.]
@test s1ex.colormap == [2, 2, 2, 2, 2]

# Check early termination

type EndWhenNode <: AbstractDijkstraVisitor
  n::Int
end

function Graphs.include_vertex!(visitor::EndWhenNode, u, v, d)
  v != visitor.n
end

s1b = dijkstra_shortest_paths(g1, eweights1, [1], visitor=EndWhenNode(5))

@test s1b.parents == [1, 3, 1, 3, 3]
@test s1b.dists == [0.0, 8.0, 5.0, 14.0, 7.0]
@test s1b.colormap == [2, 1, 2, 1, 2]

# g2: the example in Wikipedia
g2 = simple_inclist(6, is_directed=false)

g2_wedges = [
    (5, 6, 9.),
    (5, 4, 6.),
    (6, 3, 2.),
    (4, 3, 11.),
    (6, 1, 14.),
    (3, 1, 9.),
    (3, 2, 10.),
    (4, 2, 15.),
    (1, 2, 7.) ]

ne = length(g2_wedges)
eweights2 = zeros(ne)
for i = 1 : ne
    we = g2_wedges[i]
    add_edge!(g2, we[1], we[2])
    eweights2[i] = we[3]
end

@assert num_vertices(g2) == 6
@assert num_edges(g2) == 9

s2 = dijkstra_shortest_paths(g2, eweights2, 1)

@test s2.parents == [1, 1, 1, 3, 6, 3]
@test s2.dists == [0., 7., 9., 20., 20., 11.]
@test s2.colormap == [2, 2, 2, 2, 2, 2]

s2 = dijkstra_shortest_paths(g2, eweights2, [1])

@test s2.parents == [1, 1, 1, 3, 6, 3]
@test s2.dists == [0., 7., 9., 20., 20., 11.]
@test s2.colormap == [2, 2, 2, 2, 2, 2]


text:
export MPI_rank, MPI_size

MPI_rank = 0
MPI_size = 1
MPI_COMM_WORLD = nothing

mpirank(comm) = 0
mpisize(comm) = 1

function distribute(size::Int, comm_size::Int, rank::Int)
    return 1, size
end

function Allreduce(data, mpi_sum, comm)
    return data
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
using DataFrames

function zaf(
        mat::Material,
        e0::Float64,
        toa::Float64=deg2rad(40.0);
        mc::Type{<:MatrixCorrection} = XPP,
        fc::Type{<:FluorescenceCorrection}=ReedFluorescence,
        coatings::Dict{String,Film}=Dict{String,Film}(),
        stds::Dict{Element,Material}=Dict{Element,Material}(),
        minweight=0.01)
    # @info "Using the $mc Z⋅A correction and the $fc F correction."
    stdfor(elm) = haskey(stds, elm) ? stds[elm] : pure(elm)
    coatz(mat) = get(coatings,name(mat),missing)
    @assert toa>=0.0 && toa <= π/2 "The take-off angle is out-of-range [0, π/2]"
    @assert e0 > 1.0e3 "The beam energy less than 1 keV."
    res = DataFrame(Material=String[], Standard=String[], Line=CharXRay[], weight=Float64[], E₀=Float64[], Z=Float64[], A=Float64[], #
                    F=Float64[], c=Float64[], ZAFc=Float64[], k=Float64[])
    for elm in keys(mat)
        cxrs = characteristic(elm, alltransitions, 0.01, e0)
        std = stdfor(elm)
        for (sh, cxrs) in splitbyshell(cxrs)
            zaf = zafcorrection(mc, fc, Coating, mat, std, sh, e0, unkCoating=coatz(mat), stdCoating=coatz(std))
            for cxr in filter(cxr->weight(cxr)>=minweight, cxrs)
                row = [
                    mat.name, std.name, cxr, weight(cxr), e0, #
                    Z(zaf...),
                    A(zaf..., cxr, toa, toa),
                    F(zaf..., cxr, toa, toa),
                    NeXLMatrixCorrection.coating(zaf..., cxr, toa, toa),
                    ZAFc(zaf..., cxr, toa, toa),
                    k(zaf..., cxr, toa, toa)
                ]
                push!(res, row)
            end
        end
    end
    return res
end

function NeXLMatrixCorrection.zaf(ir::IterationResult; minweight=0.01)
    stds = Dict( kr.element=>kr.standard for kr in ir.kratios)
    spec = ir.label.spectrum
    NeXLMatrixCorrection.zaf(ir.comp, spec[:BeamEnergy], spec[:TakeOffAngle];
        mc=ir.iterate.mctype, fc=ir.iterate.fctype, stds=stds, minweight=minweight)
end

function NeXLUncertainties.asa(
    ::Type{DataFrame},
    krs::AbstractVector{KRatio},
    withComputedKs::Bool,
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection}=ReedFluorescence
)
    fm, sm, cm = Union{Float64,Missing}, Union{String,Missing}, Union{Material,Missing}
    xrays, mease0, meastoa, meascomp = String[], fm[], fm[], cm[]
    refe0, reftoa, refcomp, krv, dkrv, cks, ratio = fm[], fm[], cm[], Float64[], Float64[], fm[], fm[]
    for kr in krs
        push!(xrays, repr(kr.xrays))
        meas = kr.unkProps
        push!(mease0, get(meas, :BeamEnergy, missing))
        push!(meastoa, get(meas, :TakeOffAngle, missing))
        push!(meascomp, get(meas, :Composition, missing))
        ref = kr.stdProps
        push!(refe0, get(ref, :BeamEnergy, missing))
        push!(reftoa, get(ref, :TakeOffAngle, missing))
        push!(refcomp, get(ref, :Composition, missing))
        push!(krv, value(kr.kratio))
        push!(dkrv, σ(kr.kratio))
        if withComputedKs
            elm = element(kr.xrays[1])
            prim = brightest(kr.xrays)
            if any(ismissing.((meascomp[end], mease0[end], meastoa[end], refcomp[end], refe0[end], reftoa[end]))) ||
               (NeXLCore.nonneg(meascomp[end], elm) < 1.0e-6) ||
               (NeXLCore.nonneg(refcomp[end], elm) < 1.0e-6) || (energy(inner(prim))>0.95*min(mease0[end],refe0[end]))
                push!(cks, missing)
                push!(ratio, missing)
            else
                zs = zafcorrection(mc, fc, Coating, meascomp[end], kr.xrays, mease0[end])
                zr = zafcorrection(mc, fc, Coating, refcomp[end], kr.xrays, refe0[end])
                k =
                    gZAFc(zs, zr, meastoa[end], reftoa[end]) * NeXLCore.nonneg(meascomp[end], elm) /
                    NeXLCore.nonneg(refcomp[end], elm)
                push!(cks, k)
                push!(ratio, value(kr.kratio) / k)
            end
        end
    end
    res = DataFrame(
        Xrays = xrays,
        E0meas = mease0,
        TOAmeas = meastoa,
        Cmeas = meascomp,
        E0ref = refe0,
        TOAref = reftoa,
        Cref = refcomp,
        K = krv,
        ΔK = dkrv,
    )
    if withComputedKs
        res[:, :Kxpp] = cks
        res[:, :Ratio] = ratio
    end
    return res
end


text:
import Mill: catobs, MaybeHotMatrix
"""
	ExtractCategorical(s::Entry)
	ExtractCategorical(s::UnitRange)
	ExtractCategorical(s::Vector)

Converts a single item to a one-hot encoded vector.
Converts array of items into matrix of one-hot encoded columns.
There is always alocated an extra element for a unknown value.
If passed `missing`, returns column of missing values.

# Examples

```jldoctest
julia> e = ExtractCategorical(2:4);

julia> e([2,3,1,4]).data
4×4 Mill.MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}:
 1  0  0  0
 0  1  0  0
 0  0  0  1
 0  0  1  0

julia> e([1,missing,5]).data
4×3 Mill.MaybeHotMatrix{Union{Missing, Int64},Array{Union{Missing, Int64},1},Int64,Union{Missing, Bool}}:
 false  missing  false
 false  missing  false
 false  missing  false
  true  missing   true

julia> e(4).data
4×1 Mill.MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}:
 0
 0
 1
 0

julia> e(missing).data
4×1 Mill.MaybeHotMatrix{Missing,Array{Missing,1},Int64,Missing}:
 missing
 missing
 missing
 missing
```
"""
struct ExtractCategorical{V,I} <: AbstractExtractor
	keyvalemap::Dict{V,I}
	n::Int
end

ExtractCategorical(s::Entry) = ExtractCategorical(collect(keys(s.counts)))
ExtractCategorical(s::UnitRange) = ExtractCategorical(collect(s))
function ExtractCategorical(ks::Vector)
	if isempty(ks)
		@warn "Skipping initializing empty categorical variable does not make much sense to me"
		return nothing
	end
	ks = sort(unique(ks))
	ExtractCategorical(Dict(zip(ks, 1:length(ks))), length(ks) +1)
end

function (s::ExtractCategorical{V,I})(v::V) where {V,I}
    x = MaybeHotMatrix([get(s.keyvalemap, v, s.n)], s.n)
    ArrayNode(x)
end

function (s::ExtractCategorical{V,I})(vs::Vector{V}) where {V,I}
	x = MaybeHotMatrix([get(s.keyvalemap, v, s.n) for v in vs], s.n)
	ArrayNode(x)
end

# following 2 methods are to let us extract float from int extractor and vice versa
function (s::ExtractCategorical{U,I})(v::V) where {U<:Number,V<:Number,I}
    x = MaybeHotMatrix([get(s.keyvalemap, v, s.n)], s.n)
    ArrayNode(x)
end

function (s::ExtractCategorical{U,I})(vs::Vector{V}) where {U<:Number,V<:Number,I}
	x = MaybeHotMatrix([get(s.keyvalemap, v, s.n) for v in vs], s.n)
	ArrayNode(x)
end

function (s::ExtractCategorical{V,I})(vs::Vector{<:Union{V, Missing, Nothing}}) where {V,I}
	x = MaybeHotMatrix([ismissing(v) || isnothing(v) ? missing : get(s.keyvalemap, v, s.n) for v in vs], s.n)
	ArrayNode(x)
end

# todo: dodělat missingy, všchny nothing předělat na missing a pořádně to otestovat
(s::ExtractCategorical)(::V) where {V<:Union{Missing, Nothing}} = ArrayNode(MaybeHotMatrix([missing], s.n))
(s::ExtractCategorical)(::ExtractEmpty) = ArrayNode(MaybeHotMatrix(Vector{Int}(), s.n))
(s::ExtractCategorical)(v) = s(missing)

Base.reduce(::typeof(catobs), a::Vector{S}) where {S<:Flux.OneHotMatrix} = _catobs(a[:])
catobs(a::Flux.OneHotMatrix...) = _catobs(collect(a))
_catobs(a::AbstractArray{<:Flux.OneHotMatrix}) = Flux.OneHotMatrix(a[1].height,reduce(vcat, [i.data for i in a]))

Base.hash(e::ExtractCategorical, h::UInt) = hash((e.keyvalemap, e.n), h)
Base.:(==)(e1::ExtractCategorical, e2::ExtractCategorical) = e1.keyvalemap == e2.keyvalemap && e1.n === e2.n


text:
module Test_dlyap

using LinearAlgebra
using MatrixEquations
using Test

@testset "Testing discrete Lyapunov equation solvers" begin

n = 10
Ty = Float64

@testset "Discrete Lyapunov equations" begin

reltol = eps(float(1000))
a = -2+im; b = 2; @time x = lyapd(a,b)
@test abs(a*x*a'-x+b) < reltol

reltol = eps(float(1000f0))
a = 1f0-2f0im; b = 2f0; @time x = lyapd(a,b)
@test abs(a*x*a'-x+b)  < reltol

try
   x = lyapd(ones(1,1),ones(1,1))
   @test false
 catch
   @test true
 end

 try
   x = lyapd([0 1;  -1 0 ],ones(2,2))
   @test false
 catch
   @test true
 end

  try
   x = lyapd([0.5 1;  0 2],ones(2,2))
   @test false
 catch
   @test true
 end

 try
   x = lyapd([1 1 1 1;  -1 1 0 1; 0 0 0.5 0.5; 0 0 -0.5 0.5],ones(4,4))
   @test false
 catch
   @test true
 end

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
c = rand(Ty,n,n)+im*rand(Ty,n,n)
qc = c'*c
Qr = real(qc)
Ty == Float64 ? reltol = eps(float(10000)) : reltol = eps(10000*n*one(Ty))

@time x = lyapd(ac,qc);
@test norm(ac*x*ac'-x+qc)/norm(x)/max(1.,norm(ac)^2) < reltol

α = 3+im; # α = 1; #SingularException
@time x = lyapd(α*I,qc);
@test norm(α*x*α'-x+qc)/norm(x) < reltol

try
   α = 1; #SingularException
   @time x = lyapd(α,qc);
catch
   @test true
end

α = 3. +im;
@time x = lyapd(α,qc);
@test norm(α*x*α'-x+qc)/norm(x) < reltol

@time x = lyapd(ac',qc);
@test norm(ac'*x*ac-x+qc)/norm(x)/max(1.,norm(ac)^2) < reltol

@time x = lyapd(ar,Qr)
@test norm(ar*x*ar'-x+Qr)/norm(x)/max(1.,norm(ar)^2) < reltol

@time x = lyapd(ar',Qr)
@test norm(ar'*x*ar-x+Qr)/norm(x)/max(1.,norm(ar)^2) < reltol

@time x = lyapd(ac,Qr);
@test norm(ac*x*ac'-x+Qr)/norm(x)/max(1.,norm(ac)^2)  < reltol

@time x = lyapd(ac',Qr);
@test norm(ac'*x*ac-x+Qr)/norm(x)/max(1.,norm(ac)^2)  < reltol

@time x = lyapd(ar,qc)
@test norm(ar*x*ar'-x+qc)/norm(x)/max(1.,norm(ar)^2) < reltol

@time x = lyapd(ar',qc)
@test norm(ar'*x*ar-x+qc)/norm(x)/max(1.,norm(ar)^2) < reltol
end
end

@testset "Discrete generalized Lyapunov equations" begin

reltol = eps(float(10000))
ac = -2+im; ec = 4+im; b = 2; @time x = lyapd(ac,ec,b)
@test abs(ac*x*ac'-ec*x*ec'+b) < reltol

try
   x = lyapd(ones(1,1),ones(1,1),ones(1,1))
   @test false
catch
   @test true
end

try
   x = lyapd([0 2;  -2 0 ],[2 0; 0 2],ones(2,2))
   @test false
catch
   @test true
end

try
   x = lyapd([2 1;  0 2],[2 0; 0 2],ones(2,2))
   @test false
catch
   @test true
end

try
   x = lyapd(2*[1 1 1 1;  -1 1 0 1; 0 0 0.5 0.5; 0 0 -0.5 0.5],2*Matrix{Float64}(I,4,4),ones(4,4))
   @test false
catch
   @test true
end

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
er = rand(Ty,n,n)
ec = er+im*rand(Ty,n,n)

c = rand(Ty,n,n)+im*rand(Ty,n,n)
qc = c'*c
Qr = real(qc)
Ty == Float64 ? reltol = eps(float(1000)) : reltol = eps(1000*n*one(Ty))

@time x = lyapd(ac,ec,qc);
@test norm(ac*x*ac'-ec*x*ec'+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

β = 3
@time x = lyapd(ac,β*I,qc);
@test norm(ac*x*ac'-β*x*β'+qc)/norm(x)/norm(ac) < reltol

β = (1+im);
@time x = lyapd(ac,β,qc);
@test norm(ac*x*ac'-β*x*β'+qc)/norm(x)/norm(ac) < reltol

α = 2+3im
@time x = lyapd(α*I,ec,qc);
@test norm(α*x*α'-ec*x*ec'+qc)/norm(x)/norm(ec) < reltol

α = 2+3im
@time x = lyapd(α,ec,qc);
@test norm(α*x*α'-ec*x*ec'+qc)/norm(x)/norm(ec) < reltol

try
   @time x = lyapd(0*ac,0*ec,qc);
   @test false
catch
   @test true
end

α = 2+3im; β = (1+im);
@time x = lyapd(α*I,β*I,qc);
@test norm(α*x*α'-β*x*β'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

try
   α = 1; β = 1; #SingularException
   @time x = lyapd(α,β,qc);
   @test false
catch
   @test true
end

α = 2+3im; β = (1+im);
@time x = lyapd(α,β,qc);
@test norm(α*x*α'-β*x*β'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

@time x = lyapd(ac',ec',qc);
@test norm(ac'*x*ac-ec'*x*ec+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

β = (1+im);
@time x = lyapd(ac',β*I,qc);
@test norm(ac'*x*ac-β*x*β'+qc)/norm(x)/norm(ac) < reltol

α = 2+3im
@time x = lyapd(α*I,ec',qc);
@test norm(α*x*α'-ec'*x*ec+qc)/norm(x)/norm(ec) < reltol

@time x = lyapd(ac',ec,qc);
@test norm(ac'*x*ac-ec*x*ec'+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac,ec',qc);
@test norm(ac*x*ac'-ec'*x*ec+qc)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ar,er,Qr);
@test norm(ar*x*ar'-er*x*er'+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er',Qr);
@test norm(ar'*x*ar-er'*x*er+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er,Qr);
@test norm(ar'*x*ar-er*x*er'+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar,er',Qr);
@test norm(ar*x*ar'-er'*x*er+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ac,ec,Qr);
@test norm(ac*x*ac'-ec*x*ec'+Qr)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac',ec',Qr);
@test norm(ac'*x*ac-ec'*x*ec+Qr)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac',ec,Qr);
@test norm(ac'*x*ac-ec*x*ec'+Qr)/norm(x)/max(norm(ac)^2,norm(ec)^2) < reltol

@time x = lyapd(ac,ec',Qr);
@test norm(ac*x*ac'-ec'*x*ec+Qr)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar,er,qc);
@test norm(ar*x*ar'-er*x*er'+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er',qc);
@test norm(ar'*x*ar-er'*x*er+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar',er,qc);
@test norm(ar'*x*ar-er*x*er'+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol

@time x = lyapd(ar,er',qc);
@test norm(ar*x*ar'-er'*x*er+qc)/norm(x)/max(norm(ar)^2,norm(er)^2) < reltol
end
end

@testset "Discrete Lyapunov equations - Schur form" begin

for Ty in (Float64, Float32)

ar = rand(Ty,n,n);
ac = rand(Ty,n,n)+im*rand(Ty,n,n);
er = rand(Ty,n,n);
ec = er+im*rand(Ty,n,n);
as, es = schur(ar,er);
acs, ecs = schur(ac,ec);

c = rand(Ty,n,n)+im*rand(Ty,n,n);
qc = c'*c;
Qr = real(qc);
Ty == Float64 ? reltol = eps(float(1000)) : reltol = eps(1000*n*one(Ty))

x = copy(Qr)
@time lyapds!(as,x);
@test norm(as*x*as'+Qr-x)/norm(x)/max(norm(as)^2,1.) < reltol

x = copy(Qr)
@time lyapds!(as,x,adj=true);
@test norm(as'*x*as+Qr-x)/norm(x)/max(norm(as)^2,1.) < reltol

x = copy(qc)
@time lyapds!(acs,ecs,x);
@test norm(acs*x*acs'+qc-ecs*x*ecs')/norm(x)/max(norm(acs)^2,norm(ecs)^2) < reltol

x = copy(qc)
@time lyapds!(acs,I,x);
@test norm(acs*x*acs'+qc-x)/norm(x)/max(norm(acs)^2,norm(ecs)^2) < reltol

x = copy(qc)
@time lyapds!(acs,ecs,x,adj=true);
@test norm(acs'*x*acs+qc-ecs'*x*ecs)/norm(x)/max(norm(acs)^2,norm(ecs)^2) < reltol

x = copy(Qr)
@time lyapds!(as,es,x);
@test norm(as*x*as'+Qr-es*x*es')/norm(x)/max(norm(as)^2,norm(es)^2) < reltol

x = copy(Qr)
@time lyapds!(as,es,x,adj=true);
@test norm(as'*x*as+Qr-es'*x*es)/norm(x)/max(norm(as)^2,norm(es)^2) < reltol

end
end

end

end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
using ForwardDiff,SparseArrays

include(string(pwd(),"/src/gensys.jl"))
include(string(pwd(),"/misc/aux_matrix.jl"))

function dF1_theta(A,dGamma_0,dGamma_1,dGamma_2)
    m = size(A,1)
    return kron(A'*I,I(m))*dGamma_0 - kron(A'^2,I(m))*dGamma_1 - dGamma_2
end

function dF1_A(A,Gamma_0,Gamma_1)
    m = size(A,1)
    b1 = Gamma_0 - Gamma_1*A
    return kron(I(m),b1) - kron(A'*I,Gamma_1)
end

function dA_theta(A,Gamma_0,Gamma_1,dGamma_0,dGamma_1,dGamma_2)
    dvecA = dF1_A(A,Gamma_0,Gamma_1)
    dTheta = dF1_theta(A,dGamma_0,dGamma_1,dGamma_2)
    return -dvecA\Array(dTheta)
end

function dB_theta(Gamma_0,Gamma_1,Gamma_3,A,dGamma_0,dGamma_1,dGamma_2,dGamma_3)
    dA = dA_theta(A,Gamma_0,Gamma_1,dGamma_0,dGamma_1,dGamma_2)
    dA[abs.(dA) .< 10*eps()] .= 0
    dA = sparse(dA)

    m = size(A,1)
    n = size(Gamma_3,2)

    b0 = Array(Gamma_0 - Gamma_1*A)
    b0[abs.(b0) .< 10*eps()] .= 0
    b0 = pinv(b0)

    b1 = kron((b0*Gamma_3)', b0)
    b2 = dGamma_0 - kron(A',I(m))*dGamma_1 - kron(I(m),Gamma_1)*dA

    return -b1*b2 + kron(I(n),b0)*dGamma_3
end

## These are the gensys matrices

#order to par
#1 alfa
#2 beta
#3 epsilon
#4 theta
#5 sig
#6 sigma: this is the std dev of the innovation
#7 phi
#8 phi_pi
#9 phi_y
#10 rho_v

GAMMA_0_foo(par) = [par[2]  0    0  0;
           1    par[5]  0  0;
           0    0    0  0;
           0    0    0  1]

function GAMMA_1_foo(par)
    THETA = (1-par[1])/(1-par[1]+par[1]*par[3])
    lamb = (1-par[4])*(1-par[2]*par[4])/par[4]*THETA
    kappa = lamb*(par[5]+(par[7]+par[1])/(1-par[1]))
     res = [ 1       -kappa  0   0;
            0        par[5]    1   0;
           -par[8]  -par[9]  1  -1;
            0        0      0   par[10]]
end

PSI_foo(par) = [0; 0; 0; par[6]]

PI_foo(par) = [par[2]  0;
      1    par[5];
      0    0;
      0    0]

## These are the matrices for the same model written as iskrev (2008)

#order to par
#1 alfa
#2 beta
#3 epsilon
#4 theta
#5 sig
#6 sigma: this is the std dev of the innovation
#7 phi
#8 phi_pi
#9 phi_y
#10 rho_v

function Gamma_0(par)
    THETA = (1-par[1])/(1-par[1]+par[1]*par[3])
    lamb = (1-par[4])*(1-par[2]*par[4])/par[4]*THETA
    kappa = lamb*(par[5]+(par[7]+par[1])/(1-par[1]))
    res = [1 -kappa 0 0;
            0       1   1/par[5] 0;
            -par[8] -par[9] 1  -1;
            0           0    0   1
            ]
    return res
end

Gamma_1(par) = [par[2] 0 0 0;
                1/par[5] 1 0 0;
                0        0 0 0;
                0        0 0 0]

Gamma_2(par) =[0 0 0 0;
               0 0 0 0;
               0 0 0 0;
               0 0 0 par[10]]

Gamma_3(par) = [0;
                0;
                0;
               par[6]]

## Diff for the model
# Now it will be a #estrut par x ## reduce form par
# In which the reduced form depends of 3 matrices, in the following order
# A, that governs the auto regressive part of the reduce form: first 16 cols
# C that govens the observation matrix in the Kalman Filter (usually should be all zeros): 4 columns
# Omega, that is BB' in which B governs the shock transmission: remaining 10 columns (it is a vech!)

function diff_mod(par;spar=true)
    Gamma0 = Gamma_0(par)
    Gamma1 = Gamma_1(par)
    Gamma2 = Gamma_2(par)
    Gamma3 = Gamma_3(par)

    G0 = GAMMA_0_foo(par)
    G1 = GAMMA_1_foo(par)
    Psi = PSI_foo(par)
    Pi = PI_foo(par)

    dGamma0 = ForwardDiff.jacobian(Gamma_0,par)
    dGamma1 = ForwardDiff.jacobian(Gamma_1,par)
    dGamma2 = ForwardDiff.jacobian(Gamma_2,par)
    dGamma3 = ForwardDiff.jacobian(Gamma_3,par)

    gen_sol = gensys(G0,G1,Psi,Pi, verbose = false)
    A = gen_sol.Theta1

    if spar == true #use sparse arrays
        dGamma0 = sparse(dGamma0)
        dGamma1 = sparse(dGamma1)
        dGamma2 = sparse(dGamma2)
        dGamma3 = sparse(dGamma3)

        Gamma0 = sparse(Gamma0)
        Gamma1 = sparse(Gamma1)
        Gamma2 = sparse(Gamma2)
        Gamma3 = sparse(Gamma3)

        A[abs.(A) .< 10*eps()] .= 0
        A = sparse(A)
    end

    dA = dA_theta(A,Gamma0,Gamma1,dGamma0,dGamma1,dGamma2)
    dB = dB_theta(Gamma0,Gamma1,Gamma3,A,dGamma0,dGamma1,dGamma2,dGamma3)

    if spar == true
        dA = sparse(dA)
    end

    return dA,dB
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using DrWatson
@quickactivate :Continuation

using Revise
using SparseArrays, LinearAlgebra, DiffEqOperators, Setfield, Parameters
using BifurcationKit
using Plots
const BK = BifurcationKit

N = 1000
Δx = 0.1
p = @dict ϵ=1.16 ν=1
params = GenericParams(N, Δx, p)

# initial condition
x = collect((1:N) * Δx);
x0 = x[500]
σ = 1
sol0 = @. cos(x) * exp(-((x - x0) / (2 * σ))^2)

plot(x, sol0)
# ok

# newton!
p = @set p[:ϵ] -= 0.1
X2 = newton(X2, RHS_SHE1D, Jacobian_SHE1D, p, N, Δx)
plot(x, X2)

S = read_state(params)

xs = getXs(S)
heatmap(xs, xs, S.u)

S = @set S.flags.useFFT = false

V = RHS_SHE(S);
heatmap(xs, xs, V)
L = sum(abs.(V)) / params.N^2
m = maximum(abs.(V))

X = Array(vec(transpose(S.u)))
Y = RHS_SHE(X, S.gp.p, N, Δx)

J = Jacobian_SHE(X, S.gp.p, N, Δx)

X2 = newton(X2, RHS_SHE, Jacobian_SHE, S.gp.p, N, Δx)

S = @set S.gp.p[:ϵ] = 1.175
heatmap(xs, xs, transpose(reshape(X2, (N, N))))


text:
# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Evaluation Criteria][Evaluation Criteria:1]]
function generateForwardPrediction(prob::LearningProblem, startTimeInd, endTimeInd)
    dt = prob.tsteps[2] - prob.tsteps[1];
    tspan = (prob.tsteps[startTimeInd], prob.tsteps[endTimeInd]);

    W = WienerProcess(0.0, prob.validationData[:,1,:]);
    sdeProb = SDEProblem(prob.fullSystemDrift_,
			 prob.modelDiff_,
			 prob.validationData[:,1,1],
			 tspan,
			 p = prob.params,
			 noise = W,
			 noise_rate_prototype=zeros(18,18));
    predicted = Array(solve(sdeProb,
			    u0 = prob.validationData[:,1,:],
			    EM(),
			    tstops = prob.tsteps,
			    p=prob.params,
			    saveat = prob.tsteps));
    predicted = permutedims(predicted, [1,3,2]);

    return predicted;
end

function createRollup(prob::LearningProblem; startTimeInd = 1, endTimeInd = nothing)
   numDims, numTsteps, numTrajs = size(prob.validationData);
   if (endTimeInd == nothing)
      endTimeInd = numTsteps;
   end
   results = Dict{String, Vector{Any}};
   predicted = generateForwardPrediction(prob, startTimeInd, endTimeInd);

   #--------------Prediction Plots--------------#
   #skewness/flatness of velocity gradient tensor
   stats_predicted = skewnessKurtosisVGT(predicted[10:18,endTimeInd,:]);
   stats_pltPredicted = scatter(stats_predicted);
   results = merge!(results, Dict("skew/kurt_predicted"=>[stats_pltPredicted]));

   #conditional mean tangents
    cmt_pltPredicted = qrConditionalMeanTangents(predicted[10:18,[1;2],:], dt = 1e-1);
    plot!(cmt_pltPredicted, xlims=(-1.0,1.0), ylims=(-1.0,1.0));
   results = merge!(results, Dict("qrCMT_predicted"=>[cmt_pltPredicted]));

   #evolved q-r pdfs
   qrPDF_predicted = qrHistogram(predicted[10:18,endTimeInd,:],nbins=30);
   qrPDF_predictedPlot = plot(qrPDF_predicted);
   results = merge!(results, Dict("qrPDF_predicted"=>[qrPDF_predictedPlot]));

   #pdfs of cosine angle between vorticity and strain-rate eigenvectors
   alignmentHists_predicted = strainVorticityAlignmentPDFs(predicted[10:18,endTimeInd,:]);
   results = merge!(results, Dict("alignmentPDFs_predicted"=>[plot(alignmentHists_predicted[i]) for i in 1:length(alignmentHists_predicted)]))

   #pdfs of s^*
   sStarHist_predicted = sStarPDF(predicted[10:18,endTimeInd,:]);
   results = merge!(results, Dict("sStarPDF_predicted"=>[plot(sStarHist_predicted)]));

   #pdf of enstrophy
   enstrophy_predicted = normalizedEnstrophyPDF(predicted[10:18,endTimeInd,:]);
   results = merge!(results, Dict("normalizedEnstrophyPDF_predicted"=>[plot(enstrophy_predicted)]));

   #pdf of dissipation
   dissipation_predicted = normalizedDissipationPDF(predicted[10:18,endTimeInd,:]);
   results = merge!(results, Dict("normalizedDissipationPDF_predicted"=>[plot(dissipation_predicted)]));

   #loss
   results = merge!(results, Dict("loss"=>[plot(prob.lossArray)]));

   #--------------Ground Truth Plots--------------#
   #skewness/flatness of velocity gradient tensor
   stats_groundTruth = skewnessKurtosisVGT(prob.validationData[10:18,endTimeInd,:]);
   stats_pltGroundTruth = scatter(stats_groundTruth);
   results = merge!(results, Dict("skew/kurt_groundTruth"=>[stats_pltGroundTruth]));

   #conditional mean tangents
   cmt_pltGroundTruth = qrConditionalMeanTangents(prob.validationData[10:18,[1;2],:], dt = 1e-1, plotToOverlay=results["qrCMT_predicted"][1]);
    results = merge!(results, Dict("qrCMT_groundTruth"=>[cmt_pltGroundTruth]));

   #evolved q-r pdfs
   qrPDF_initial = qrHistogram(prob.validationData[10:18,startTimeInd,:],nbins=30);
   qrPDF_initialPlot = plot(qrPDF_initial);
   results = merge!(results, Dict("qrPDF_initial"=>[qrPDF_initialPlot]));

   qrPDF_groundTruth = qrHistogram(prob.validationData[10:18,endTimeInd,:]);
   qrPDF_groundTruthPlot = plot(qrPDF_groundTruth);
   results = merge!(results, Dict("qrPDF_groundTruth"=>[qrPDF_groundTruthPlot]));

   #pdfs of cosine angle between vorticity and strain-rate eigenvectors
   alignmentHists_groundTruth = strainVorticityAlignmentPDFs(prob.validationData[10:18,endTimeInd,:]);
   results = merge!(results, Dict("alignmentPDFs_groundTruth"=>[plot(alignmentHists_groundTruth[i]) for i in 1:length(alignmentHists_groundTruth)]))

   #pdfs of s^*
   sStarHist_groundTruth = sStarPDF(prob.validationData[10:18,endTimeInd,:]);
   results = merge!(results, Dict("sStarPDF_groundTruth"=>[plot(sStarHist_groundTruth)]));

   #pdf of enstrophy
   enstrophy_groundTruth = normalizedEnstrophyPDF(prob.validationData[10:18,endTimeInd,:]);
   results = merge!(results, Dict("normalizedEnstrophyPDF_groundTruth"=>[plot(enstrophy_groundTruth)]));

   #pdf of dissipation
   dissipation_groundTruth = normalizedDissipationPDF(prob.validationData[10:18,endTimeInd,:]);
   results = merge!(results, Dict("normalizedDissipationPDF_groundTruth"=>[plot(dissipation_groundTruth)]));

   return results;
end
# Evaluation Criteria:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Skewness & Kurtosis of predicted velocity gradient tensor][Skewness & Kurtosis of predicted velocity gradient tensor:1]]
function vorticityElements(Aij::Array{Float64,2})
   numDims, numTrajs = size(Aij);
   @assert numDims==9;

   arr = zeros(3,numTrajs);
   for i in 1:numTrajs
      M = reshape(Aij[:,i],(3,3));
      W = 0.5(M-M');
      arr[:,i] .= [W[3,2]; W[1,3]; W[2,1]];
   end

   return reshape(arr, prod(size(arr)));
end
function onDiagonalElements(Aij::Array{Float64,2})
   numDims, numTrajs = size(Aij);
   @assert numDims==9;

   arr = zeros(3,numTrajs);
   for i in 1:numTrajs
      M = reshape(Aij[:,i],(3,3));
      arr[:,i] .= [M[1,1]; M[2,2]; M[3,3]];
   end
   return reshape(arr, prod(size(arr)));
end
function skewnessKurtosisVGT(Aij::Array{Float64,2})
   onDiagonal = onDiagonalElements(Aij);
   offDiagonal = vorticityElements(Aij);

   return Dict("onDiagonalSkewness"=>skewness(onDiagonal),
	       "onDiagonalKurtosis"=>kurtosis(onDiagonal),
	       "offDiagonalSkewness"=>skewness(offDiagonal),
	       "offDiagonalKurtosis"=>kurtosis(offDiagonal));
end
# Skewness & Kurtosis of predicted velocity gradient tensor:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Conditional mean tangents][Conditional mean tangents:1]]
# Aij needs to be of size (9, 2, numTrajs), representing 2 timesteps (start/end)
function qrConditionalMeanTangents(Aij::Array{Float64,3};
				   nbins_r=15,
				   nbins_q=15,
				   dt=3e-4,
				   Rrange=range(-sqrt(3)/9, sqrt(3)/9, length=nbins_r),
				   Qrange=range(-0.5, 0.5, length=nbins_q),
				   label = nothing,
				   pltToOverlay = nothing)

   numDims, numTsteps, numTrajs = size(Aij);
   @assert numDims == 9;
   rmid = Array(midpoints(Rrange));
   qmid = Array(midpoints(Qrange));
   seperatix_r = range(Rrange[1], Rrange[end], length=1000);
   seperatix_q = -((27.0/4.0).*(seperatix_r.^2)).^(1/3);

   function calc_rqFromSample(Marr)
      @assert length(Marr) == 9;
      local M = reshape(Marr, (3,3));
      local b = M/norm(M);
      local r = -tr(b^3)/3.0;
      local q = -tr(b^2)/2.0;
      return [r; q];
   end

   function inRange(M, R_bounds, Q_bounds)
      local r,q = calc_rqFromSample(M);
      return (R_bounds[1] <= r < R_bounds[2]) & (Q_bounds[1] <= q < Q_bounds[2]);
   end

   function find(M_vals, R_bounds, Q_bounds)
      numTrajs = size(M_vals)[2];
      trueIdx = [inRange(M_vals[:,i], R_bounds, Q_bounds) for i in 1:numTrajs];
      return findall(trueIdx);
   end

   r = zeros(length(rmid), length(qmid));
   q = zeros(length(rmid), length(qmid));
   delta = zeros(length(rmid), length(qmid), 2);

   for i in 1:length(rmid)
      for j in 1:length(qmid)
	 r[i,j] = rmid[i];
	 q[i,j] = qmid[j];
      end
   end

   Threads.@threads for i in 1:length(rmid)
      for j in 1:length(qmid)
	 samples = find(Aij[:,1,:],
			(Rrange[i], Rrange[i+1]),
			(Qrange[j], Qrange[j+1]));
	 for l in 1:length(samples)
	    end_rq = calc_rqFromSample(Aij[:,end,samples[l]]);
	    start_rq = calc_rqFromSample(Aij[:,1,samples[l]]);
	    delta[i,j,:] .+= (end_rq .- start_rq)/dt;
	 end
	 delta[i,j,:] /= length(samples);
      end
   end

   r = reshape(r, prod(size(r)));
   q = reshape(q, prod(size(q)));
   delta = Flux.flatten(delta);

   if (pltToOverlay == nothing)
      plt = plot(seperatix_r, seperatix_q, label="seperatix");
      quiver!(plt, r,q, quiver=(delta[:,1], delta[:,2]), label=label);
      return plt;
   else
      quiver!(pltToOverlay, r,q, quiver=(delta[:,1], delta[:,2]),label=label);
      return pltToOverlay;
   end
end
# Conditional mean tangents:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Evolved $q\text{-}r$ PDFs][Evolved $q\text{-}r$ PDFs:1]]
# Aij has size = (9,numTrajs)
function qrHistogram(Aij::Array{Float64,2}; nbins = 15)

    numDims, numTrajs = size(Aij);
    @assert numDims == 9;

    r = zeros(numTrajs);
    q = zeros(numTrajs);

    Threads.@threads for i in 1:numTrajs
	M = reshape(Aij[:,i],(3,3));
	b = M/norm(M);
	q[i] = -tr(b^2)/2.0;
	r[i] = -tr(b^3)/3.0;
    end

    h = fit(Histogram, (r, q), nbins = nbins);
    return h;
end
# Evolved $q\text{-}r$ PDFs:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*PDFs of cosine of angle between vorticity and strain-rate eigenvectors][PDFs of cosine of angle between vorticity and strain-rate eigenvectors:1]]
#Aij is of size (9,numTrajs)
function strainVorticityAlignmentPDFs(Aij::Array{Float64,2}; nbins = 15)
   numDims, numTrajs = size(Aij);
   @assert numDims == 9;

   stats_e1 = zeros(numTrajs);
   stats_e2 = zeros(numTrajs);
   stats_e3 = zeros(numTrajs);

   function costheta(e_i, omega)
      return (dot(e_i,omega))/(norm(e_i)*norm(omega));
   end

   Threads.@threads for i in 1:numTrajs
      M = reshape(Aij[:,i], (3,3));
      S = 0.5*(M+M');
      W = 0.5*(M-M');
      omega = 2*[ W[3,2]; W[1,3]; W[2,1] ];
      vecs = eigen(S).vectors;
      e1 = vecs[:,3];
      e2 = vecs[:,2];
      e3 = vecs[:,1];
      stats_e1[i] = costheta(e1,omega);
      stats_e2[i] = costheta(e2,omega);
      stats_e3[i] = costheta(e3,omega);
   end
   h1 = fit(Histogram, stats_e1, nbins=nbins);
   h2 = fit(Histogram, stats_e2, nbins=nbins);
   h3 = fit(Histogram, stats_e3, nbins=nbins);

   return [h1, h2, h3];
end
# PDFs of cosine of angle between vorticity and strain-rate eigenvectors:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*PDF of $s^*$][PDF of $s^*$:1]]
#Aij is of size (9,numTrajs)
function sStarPDF(Aij::Array{Float64,2}; nbins = 15)
   numDims, numTrajs = size(Aij);
   @assert numDims == 9;
   stats_sStar = zeros(numTrajs);

   Threads.@threads for i in 1:numTrajs
      M = reshape(Aij[:,i], (3,3));
      S = 0.5*(M+M');
      values = eigen(S).values;
      alpha = values[3];
      beta  = values[2];
      gamma = values[1];
      stats_sStar[i] = -(3*sqrt(6)*alpha*beta*gamma)/((alpha^2+beta^2+gamma^2)^(3/2));
   end
   h = fit(Histogram, stats_sStar, nbins=nbins);

   return h;
end
# PDF of $s^*$:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*PDF of enstrophy][PDF of enstrophy:1]]
function normalizedEnstrophyPDF(Aij::Array{Float64,2}; nbins::Int = 15)
   numDims, numTrajs = size(Aij);
   @assert numDims == 9;
   stats = [sum(abs2,Aij[:,i]) for i in 1:numTrajs];
   meanEnstrophy = mean(stats);
   stats /= meanEnstrophy;
   h = fit(Histogram, stats, nbins=nbins);
   return h;
end
# PDF of enstrophy:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*PDF of dissipation][PDF of dissipation:1]]
#Aij is of size (9,numTrajs)
function normalizedDissipationPDF(Aij::Array{Float64,2}; nu::Float64 = 1.1277e-3, nbins::Int =15)
   numDims, numTrajs = size(Aij);
   @assert numDims == 9;
   stats = zeros(numTrajs);

   Threads.@threads for i in 1:numTrajs
      M = reshape(Aij[:,i], (3,3));
      S = 0.5*(M+M');
      stats[i] = nu*(norm(S))^2
   end
   meanDissipation = mean(stats);
   stats /= meanDissipation;
   h = fit(Histogram, stats, nbins=nbins);

   return h;
end
# PDF of dissipation:1 ends here

# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Loss][Loss:1]]
function plotLoss(prob::LearningProblem)
   p = plot(prob.lossArray, yaxis=log);
   return p;
end
# Loss:1 ends here

max_stars_repo_name: cmhyett/MLCoarseGrainedVGT

max_issues_repo_name: cmhyett/MLCoarseGrainedVGT

max_forks_repo_name: cmhyett/MLCoarseGrainedVGT


text:
import FITSIO

function get_fits_members(fits_name)
include("/Users/alex/CODES/Julia_Custom_jl/alex/src/remove_char.jl")

# fits_name = "/Users/alex/Desktop/make_pros/old/completeness_test_results.fits"

  # open the fits_file
fits_data=FITSIO.fits_open_data(fits_name)

  # get the header information
fits_total_header=FITSIO.fits_hdr2str(fits_data)
  # split the header into individual strings
fits_header_items = split(fits_total_header)

  # search for NAXIS2, TFIELDS and GCOUNT in the header
  # in order to get the rowsxcolumns of the data
tfields = 0.0
naxis2 = 0.0
gcount = 0.0
  for i in 1:size(fits_header_items,1)
      # get TFIELDS
       if fits_header_items[i] == "TFIELDS"
            tfields = float(fits_header_items[i + 2])
      end

      # get NAXIS2
      if fits_header_items[i] == "NAXIS2"
        naxis2 = float(fits_header_items[i + 2])
      end

      # get gcount
      if fits_header_items[i] == "GCOUNT"
        gcount = float(fits_header_items[i + 2])
      end

  end

##########################################################################
#### get and create each TFIELD string names                  ############
#### which are the structure member names                     ############
#### to be used in the definition of the return type          ############
##########################################################################
                                                              ############
  # ttype_names will hold the name of each ttype              ############
ttype_names=Array(String, int64(tfields))                     ############
  # ttype_array_call will hold the  TTYPE id'                 ############
ttype_array_call=Array(String, int64(tfields))                ############
                                                              ############
    # populate ttype_array_call with the ttype identifiers    ############
  for i in 1:tfields                                          ############
    ttype_array_call[i] = string("TTYPE", string(int64(i)))   ############
  end                                                         ############
                                                              ############
    # search for all TTYPE's as populated in ttype_array_call       ######
    # and pass them to ttype_names                            ############
  for i in 1:size(fits_header_items,1)                        ############
    # i cycles through all individual header items            ############
      for x in 1:tfields                                      ############
        # x cycles trhough ttype_array_call to search for each TTYPE   ###
          if fits_header_items[i] == ttype_array_call[x]      ############
            ttype_names[x]=fits_header_items[i+2]             ############
          end                                                 ############
      end                                                     ############
  end                                                         ############
                                                              ############
##########################################################################

##########################################################################
#### get and create each TFORM string names                   ############
#### which are the structure member names                     ############
#### to be used in the definition of the return type          ############
##########################################################################
                                                              ############
  # tform_names will hold the name of each tform              ############
tform_names=Array(String, int64(tfields))                     ############
  # tform_array_call will hold the  tform id'                 ############
tform_array_call=Array(String, int64(tfields))                ############
                                                              ############
    # populate tform_array_call with the tform identifiers    ############
  for i in 1:tfields                                          ############
    tform_array_call[i] = string("TFORM", string(int64(i)))   ############
  end                                                         ############
                                                              ############
    # search for all tform's as populated in tform_array_call       ######
    # and pass them to tform_names                            ############
  for i in 1:size(fits_header_items,1)                        ############
    # i cycles through all individual header items            ############
      for x in 1:tfields                                      ############
        # x cycles trhough tform_array_call to search for each tform   ###
          if fits_header_items[i] == tform_array_call[x]      ############
            tform_names[x]=fits_header_items[i+2]             ############
          end                                                 ############
      end                                                     ############
  end                                                         ############
                                                              ############
##########################################################################

##########################################################################
####                                                          ############
####    the problem now is that tform_names and ttype_names   ############
####    could be "polluted" by "'"                            ############
####                                                          ############
##########################################################################
                                                              ############
tform_names_final=Array(String, int64(tfields))               ############
ttype_names_final=Array(String, int64(tfields))               ############
                                                              ############
                                                              ############
  for i in 1:tfields                                          ############
    tform_names_final[i] = remove_char(tform_names[i], "'")   ############
    ttype_names_final[i] = remove_char(ttype_names[i], "'")   ############
  end                                                         ############
                                                              ############
##########################################################################

##########################################################################
####                                                          ############
####     create a custom Dict to return                       ############
####     as an object                                         ############
####                                                          ############
##########################################################################
                                                              ############
return_dictionary = Dict{Symbol,Any}()                        ############
return_dictionary[:tform_names] = tform_names_final           ############
return_dictionary[:ttype_names] = ttype_names_final           ############
return_dictionary[:tfields] = tfields                         ############
return_dictionary[:naxis2] = naxis2                           ############
return_dictionary[:gcount] = gcount                           ############
                                                              ############
return return_dictionary                                      ############
##########################################################################

end


text:
""" A diagnostic type associated with FourierFlows.Problem types """
mutable struct Diagnostic{T} <: AbstractDiagnostic
  calc::Function
  prob::Problem
  num::Int
  data::Array{T,1}
  time::Array{Float64,1}
  step::Array{Int,1}
  value::T
  count::Int
  freq::Int
end

""" Constructor for the ProblemDiagnostic type. """
function Diagnostic(calc::Function, prob::FourierFlows.Problem; freq=1, nsteps=1, num=ceil(Int, (nsteps+1)/freq))
  value = calc(prob)
  T = typeof(value)
  data = Array{T}(undef, num)
  time = Array{Float64}(undef, num)
  step = Array{Int}(undef, num)
  data[1] = value
  time[1] = prob.t
  step[1] = prob.step
  Diagnostic{T}(calc, prob, num, data, time, step, value, 1, freq)
end

getindex(d::Diagnostic, inds...) = getindex(d.data, inds...)

"""
    resize!(diag, newnum)

Resize the Diagnostic data and time arrays to length newnum.
"""
function resize!(diag::AbstractDiagnostic, newnum::Int)
  resize!(diag.data, newnum)
  resize!(diag.time, newnum)
  resize!(diag.step, newnum)
  nothing
end

"""
    update!(diag)

Update diag with its current value.
"""
function update!(diag::AbstractDiagnostic)
  diag.data[diag.count] = diag.calc(diag.prob)
  diag.time[diag.count] = diag.prob.t
  diag.step[diag.count] = diag.prob.step
  diag.value = diag.data[diag.count]
  nothing
end

function update!(diags::AbstractArray)
  for diag in diags
    update!(diag)
  end
  nothing
end

"""
    increment!(diag)
    increment!(diags)

Increment the Diagnostic diag, or an array of Diagnostics diags.
"""
function increment!(diag::AbstractDiagnostic)
  if diag.count < diag.num
    diag.data[diag.count+1] = diag.calc(diag.prob)
    diag.time[diag.count+1] = diag.prob.t
    diag.step[diag.count+1] = diag.prob.step
  else
    push!(diag.data, diag.calc(diag.prob))
    push!(diag.time, diag.prob.t)
    push!(diag.step, diag.prob.step)
    diag.num += 1
  end
  diag.count += 1
  diag.value = diag.data[diag.count]
  nothing
end

function increment!(diags::AbstractArray)
  for d in diags
    increment!(d)
  end
  nothing
end