text:
function inlet_broadband(settings::PyObject, theta, M_tip, tsqem, M_d_fan::Float64, RSS_fan::Float64)

    T = eltype(M_tip)
    # Fan inlet broadband noise component:
    if settings.fan_BB_method == "original"
        # Tip Mach-dependent term (F1 of Eqn 4 in report, Figure 4A):
        if M_d_fan <= 1
            if M_tip <= 0.9
                F1IB = 58.5
            else
                F1IB = 58.5 - 20 * log10(M_tip / 0.9)
            end
        else
            if M_tip <= 0.9
                F1IB = 58.5 + 20 * log10(M_d_fan)
            else
                F1IB = 58.5 + 20 * log10(M_d_fan) - 20 * log10(M_tip / 0.9)
            end
        end

        # Rotor-stator correction term (F2 of Eqn 4, Figure 6B):
        if settings.fan_id == false
            F2IB = -5 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan <= 100
                F2IB = -5 * log10(RSS_fan / 300)
            else
                F2IB = -5 * log10(100 / 300)  # This is set to a constant 2.3856
            end
        end
    elseif settings.fan_BB_method == "allied_signal"
        # Tip Mach-dependent term (F1 of Eqn 4 in report, Figure 4A, modified by AlliedSignal):
        if M_d_fan <= 1
            if M_tip <= 0.9
                F1IB = 55.5
            else
                F1IB = 55.5 - 20 * log10(M_tip / 0.9)
            end
        else
            if M_tip <= 0.9
                F1IB = 55.5 + 20 * log10(M_d_fan)
            else
                F1IB = 55.5 + 20 * log10(M_d_fan) - 20 * log10(M_tip / 0.9)
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 4, Figure 6B):
        if settings.fan_id == false
            F2IB = -5 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan <= 100
                F2IB = -5 * log10(RSS_fan / 300)
            else
                F2IB = -5 * log10(100. / 300)  # This is set to a constant 2.3856
            end
        end
    elseif settings.fan_BB_method == "geae"
        # Tip Mach-dependent term (F1 of Eqn 4 in report, Figure 4A, modified by GE):
        if M_d_fan <= 1
            if M_tip <= 0.9
                F1IB = 58.5
            else
                F1IB = 58.5 - 50 * log10(M_tip / 0.9)
            end
        else
            if M_tip <= 0.9
                F1IB = 58.5 + 20 * log10(M_d_fan)
            else
                F1IB = 58.5 + 20 * log10(M_d_fan) - 50 * log10(M_tip / 0.9)
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 4, Figure 6B):
        F2IB = zeros(T, 1)
    elseif settings.fan_BB_method == "kresja"
        # Tip Mach-dependent term (F1, of Eqn 4 in report, Figure 4A, modified by Krejsa):
        if M_d_fan <= 1
            if M_tip < 0.72
                F1IB = 34 + 20 * log10(1. / 1.245)
            else
                F1IB = 34 - 43 * (M_tip - 0.72) + 20 * nlog10(1. / 1.245)
            end
        else
            if M_tip < 0.72
                F1IB = 34 + 20 * log10(M_d_fan / 1.245)
            else
                F1IB = 34 - 43 * (M_tip - 0.72) + 20 * log10(M_d_fan/ 1.245)
            end
        end

        # Rotor-stator spacing correction term (F2, of Eqn 4, Figure 6B):
        if settings.fan_id == false
            F2IB = -5 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan <= 100
                F2IB = -5 * log10(RSS_fan / 300)
            else
                F2IB = -5 * log10(100 / 300)  # This is set to a constant 2.3856
            end
        end
    else
        throw(DomainError("Invalid fan_BB_method specified. Specify: original / allied_signal / geae / kresja."))
    end

    # Theta correction term (F3 of Eqn 4, Figure 7A):
    if settings.fan_BB_method == "kresja"
        # Krejsa method:
        THET7A = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190]
        FIG7A = [-0.5, -1, -1.25, -1.41, -1.4, -2.2, -4.5, -8.5, -13, -18.5, -24, -30, -36, -42, -48, -54, -60,-66, -73, -66]
        f_F3IB = LinearInterpolation(THET7A, FIG7A)
        F3IB = f_F3IB(theta)
    elseif settings.fan_BB_method in ["original", "allied_signal", "geae"]
        # All other methods:
        THET7A = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 180, 250]
        FIG7A = [-2, -1, 0, 0, 0, -2, -4.5, -7.5, -11, -15, -19.5, -25, -63.5, -25]
        f_F3IB = LinearInterpolation(THET7A, FIG7A)
        F3IB = f_F3IB(theta)
    else
        throw(DomainError("Invalid fan_BB_method specified. Specify: original / allied_signal / geae / kresja."))
    end
    # Component value:
    spl_IB = tsqem .+ F1IB .+ F2IB .+ F3IB

    return spl_IB
end
function discharge_broadband(settings::PyObject, theta, M_tip, tsqem, M_d_fan::Float64, RSS_fan::Float64)

    T = eltype(M_tip)
    # Fan discharge broadband noise component:
    if settings.fan_BB_method == "original"
        # Tip Mach-dependent term (F1 of Eqn 10 in report, Figure 4B):
        if M_d_fan <= 1
            if M_tip <= 1
                F1DB = 60
            else
                F1DB = 60 - 20 * log10(M_tip / 1)
            end
        else
            if M_tip <= 1
                F1DB = 60 + 20 * log10(M_d_fan)
            else
                F1DB = 60 + 20 * log10(M_d_fan) - 20 * log10(M_tip / 1)
            end
        end

        # Rotor-stator correction term (F2 of Eqn 4, Figure 6B):
        if settings.fan_id == false
            F2DB = -5 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan <= 100
                F2DB = -5 * log10(RSS_fan / 300)
            else
                F2DB = -5 * log10(100 / 300)  # This is set to a constant 2.3856
            end
        end
    elseif settings.fan_BB_method == "allied_signal"
        # Tip Mach-dependent term (F1 of Eqn 10 in report, Figure 4B, modified by AlliedSignal):
        if M_d_fan <= 1
            if M_tip <= 1
                F1DB = 58
            else
                F1DB = 58 - 20 * log10(M_tip / 1)
            end
        else
            if M_tip <= 1
                F1DB = 58 + 20 * log10(M_d_fan)
            else
                F1DB = 58 + 20 * log10(M_d_fan) - 20 * log10(M_tip / 1)
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 10, Figure 6B, modified by AlliedSignal):
        if settings.fan_id == false
            F2DB = -5 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan <= 100
                F2DB = -5 * log10(RSS_fan / 300)
            else
                F2DB = -5 * log10(100 / 300)  # This is set to a constant 2.3856
            end
        end
    elseif settings.fan_BB_method == "geae"
        # Tip Mach-dependent term (F1 of Eqn 10 in report, Figure 4B, modified by GE):
        if M_d_fan <= 1
            if M_tip <= 1
                F1DB = 63
            else
                F1DB = 63 - 30 * log10(M_tip / 1)
            end
        else
            if M_tip <= 1
                F1DB = 63 + 20 * log10(M_d_fan)
            else
                F1DB = 63 + 20 * log10(M_d_fan) - 30 * log10(M_tip / 1)
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 10, Figure 6B, modified by GE):
        F2DB = -5 * log10(RSS_fan / 300)
    elseif settings.fan_BB_method == "kresja"
        # Tip Mach-dependent term (F1, of Eqn 10 in report, Figure 4B, modified by Krejsa):
        if M_d_fan <= 1
            # If M_tip < 0.65 Then
            #   F1DBBkrejsa = 34 + 20 * np.log10(1 / 1.245)
            # Else
            F1DB = 34 - 17 * (M_tip - 0.65) + 20 * log10(1 / 1.245)
        else
            # If M_tip < 0.65 Then
            #   F1DBBkrejsa = 34 + 20 * np.log10(M_d_fan / 1.245)
            # Else
            F1DB = 34 - 17 * (M_tip - 0.65) + 20 * log10(M_d_fan / 1.245)
        end

        # Rotor-stator spacing correction term (F2, of Eqn 10, Figure 6B):
        if settings.fan_id == false
            F2DB = -5 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan <= 100
                F2DB = -5 * log10(RSS_fan / 300)
            else
                F2DB = -5 * log10(100 / 300)  # This is set to a constant 2.3856
            end
        end
    else
        throw(DomainError("Invalid fan_BB_method specified. Specify: original / allied_signal / geae / kresja."))
    end

    # Theta correction term (F3 of Eqn 10, Figure 7B):
    if settings.fan_BB_method == "allied_signal"
        THET7B = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
        FIG7B = [0, -29.5, -26, -22.5, -19, -15.5, -12, -8.5, -5, -3.5, -2.5, -2, -1.3, 0, -3, -7, -11, -15, -20]
        f_F3DB = LinearInterpolation(THET7B, FIG7B)
        F3DB = f_F3DB(theta)
    elseif settings.fan_BB_method == "kresja"
        THET7B = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190]
        FIG7B = [-30, -25, -20.8, -19.5, -18.4, -16.7, -14.5, -12, -9.6, -6.9, -4.5, -1.8, -0.3, 0.5, 0.7, -1.9,-4.5, -9, -15, -9]
        f_F3DB = LinearInterpolation(THET7B, FIG7B)
        F3DB = f_F3DB(theta)
    elseif settings.fan_BB_method in ["original", "geae"]
        THET7B = [0, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190]
        FIG7B = [-41.6, -15.8, -11.5, -8, -5, -2.7, -1.2, -0.3, 0, -2, -6, -10, -15, -20, -15]
        f_F3DB = LinearInterpolation(THET7B, FIG7B)
        F3DB = f_F3DB(theta)
    else
        throw(DomainError("Invalid fan_BB_method specified. Specify: original / allied_signal / geae / kresja."))
    end

    # Added noise factor if there are inlet guide vanes present:
    if settings.fan_igv == true
        CDB = 3 * ones(T, 1)
    else
        CDB = zeros(T, 1)
    end

    # Component value:
    spl_DB = tsqem .+ F1DB .+ F2DB .+ F3DB .+ CDB

    return spl_DB
end
function inlet_tones(settings::PyObject, theta, M_tip, tsqem, M_d_fan::Float64, RSS_fan::Float64)

    T = eltype(M_tip)
    # Fan inlet discrete tone noise component:
    if settings.fan_RS_method == "original"  # Original method:
        # Tip Mach-dependent term (F1 of Eqn 6 in report, Figure 10A):
        if M_d_fan <= 1
            if M_tip <= 0.72
                F1TI = 60.5
            else
                F1TIA = 60.5 + 50 * log10(M_tip / 0.72)
                F1TIB = 59.5 + 80 * log10(1. / M_tip)
                if F1TIA < F1TIB
                    F1TI = F1TIA
                else
                    F1TI = F1TIB
                end
            end
        else
            if M_tip <= 0.72
                F1TI = 60.5 + 20 * log10(M_d_fan)
            else
                F1TIA = 60.5 + 20 * log10(M_d_fan) + 50 * log10(M_tip / 0.72)
                # Note the 1975 version of NASA TMX-71763 writes it this way:
                # F1TIB = 59.5 + 20 * log10(M_d_fan) + 80 * log10(M_d_fan / M_tip)
                # But the 1979 version of NASA TMX-71763 writes it this way:
                F1TIB = 59.5 + 80 * log10(M_d_fan / M_tip)
                if F1TIA < F1TIB
                    F1TI = F1TIA
                else
                    F1TI = F1TIB
                end
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 6, Figure 12):
        if settings.fan_id == false
            F2TI = -10 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan < 100
                F2TI = -10 * log10(RSS_fan / 300)
            else
                F2TI = -10 * log10(100 / 300)  # This is set to a constant 4.7712
            end
        end
    elseif settings.fan_RS_method == "allied_signal"
        # Tip Mach-dependent term (F1 of Eqn 6 in report, Figure 10A, modified by AlliedSignal):
        if M_d_fan <= 1
            if M_tip <= 0.72
                F1TI = 54.5
            else
                F1TIA = 54.5 + 50 * log10(M_tip / 0.72)
                F1TIB = 53.5 + 80 * log10(1. / M_tip)
                if F1TIA < F1TIB
                    F1TI = F1TIA
                else
                    F1TI = F1TIB
                end
            end
        else
            if M_tip <= 0.72
                F1TI = 54.5 + 20 * log10(M_d_fan)
            else
                F1TIA = 54.5 + 20 * log10(M_d_fan) + 50 * log10(M_tip / 0.72)
                F1TIB = 53.5 + 80 * log10(M_d_fan / M_tip)
                if F1TIA < F1TIB
                    F1TI = F1TIA
                else
                    F1TI = F1TIB
                end
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 6, Figure 12, modified by AlliedSignal):
        if settings.fan_id == false
            F2TI = -10 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan < 100
                F2TI = -10 * log10(RSS_fan / 300)
            else
                F2TI = -10 * log10(100. / 300)  # This is set to a constant 4.7712
            end
        end
    elseif settings.fan_RS_method == "geae"
        # Tip Mach-dependent term (F1 of Eqn 6 in report, Figure 10A, modified by GE):
        if M_d_fan <= 1
            if M_tip <= 0.72
                F1TI = 60.5
            else
                F1TIA = 60.5 + 50 * log10(M_tip / 0.72)
                F1TIB = 64.5 + 80 * log10(1. / M_tip)
                if F1TIA < F1TIB
                    F1TI = F1TIA
                else
                    F1TI = F1TIB
                end
            end
        else
            if M_tip <= 0.72
                F1TI = 60.5 + 20 * log10(M_d_fan)
            else
                F1TIA = 60.5 + 20 * log10(M_d_fan) + 50 * log10(M_tip / 0.72)
                F1TIB = 64.5 + 80 * log10(M_d_fan) - 80 * log10(M_tip)
                if F1TIA < F1TIB
                    F1TI = F1TIA
                else
                    F1TI = F1TIB
                end
            end
        end
        # Rotor-stator spacing correction term (F2 of Eqn 6, Figure 12, modified to zero by GE):
        F2TI = zeros(T, 1)
    elseif settings.fan_RS_method == "kresja"
        # Tip Mach-dependent term (F1 of Eqn 6 in report, Figure 10A, modified by Krejsa):
        if M_d_fan <= 1
            F1TI = 42 - 20 * M_tip + 20 * log10(1. / 1.245)
        else
            F1TI = 42 - 20 * M_tip + 20 * log10(M_d_fan / 1.245)
        end

        # Rotor-stator spacing correction term (F2 of Eqn 6, Figure 12):
        if settings.fan_id == false
            F2TI = -10 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan < 100
                F2TI = -10 * log10(RSS_fan / 300)
            else
                F2TI = -10 * log10(100 / 300)  # This is set to a constant 4.7712
            end
        end
    else
        throw(DomainError("Invalid fan_RS_method specified. Specify: original / allied_signal / geae / kresja."))
    end

    # Theta correction term (F3 of Eqn 6, Figure 13A):
    if settings.fan_RS_method == "allied_signal"
        THT13A = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
        FIG13A = [-3, -1.5, -1.5, -1.5, -1.5, -2, -3, -4, -6, -9, -12.5, -16, -19.5, -23, -26.5, -30, -33.5, -37,-40.5]
        f_F3TI = LinearInterpolation(THT13A, FIG13A)
        F3TI = f_F3TI(theta)
    elseif settings.fan_RS_method == "kresja"
        THT13A = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
        FIG13A = [-3, -1.5, 0, 0, 0, -1.2, -3.5, -6.8, -10.5, -15.5, -19, -25, -32, -40, -49, -59, -70, -80, -90]
        f_F3TI = LinearInterpolation(THT13A, FIG13A)
        F3TI = f_F3TI(theta)
    elseif settings.fan_RS_method in ["original", "geae"]  # For original and GE large fan methods:
        THT13A = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 180, 260]
        FIG13A = [-3, -1.5, 0, 0, 0, -1.2, -3.5, -6.8, -10.5, -14.5, -19, -55, -19]
        f_F3TI = LinearInterpolation(THT13A, FIG13A)
        F3TI = f_F3TI(theta)
    else
        throw(DomainError("Invalid fan_RS_method specified. Specify: original / allied_signal / geae / kresja."))
    end

    # Component value:
    tonlv_I = tsqem .+ F1TI .+ F2TI .+ F3TI

    return tonlv_I
end
function discharge_tones(settings::PyObject, theta, M_tip, tsqem, M_d_fan::Float64, RSS_fan::Float64)

    T = eltype(M_tip)
    # Fan discharge discrete tone noise component:
    if settings.fan_RS_method == "original"  # Original method:
        # Tip Mach-dependent term (F1 of Eqn 12 in report, Figure 10B):
        if M_d_fan <= 1
            if M_tip <= 1
                F1TD = 63
            else
                F1TD = 63 - 20 * log10(M_tip / 1)
            end
        else
            if M_tip <= 1
                F1TD = 63 + 20 * log10(M_d_fan)
            else
                F1TD = 63 + 20 * log10(M_d_fan) - 20 * log10(M_tip / 1)
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 12, Figure 12):
        if settings.fan_id == false
            F2TD = -10 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan < 100
                F2TD = -10 * log10(RSS_fan / 300)
            else
                F2TD = -10 * log10(100 / 300)  # This is set to a constant 4.7712
            end
        end
    elseif settings.fan_RS_method == "allied_signal"
        # Tip Mach-dependent term (F1 of Eqn 12 in report, Figure 10B, modified by AlliedSignal):
        if M_d_fan <= 1
            if M_tip <= 1
                F1TD = 59
            else
                F1TD = 59 - 20 * log10(M_tip / 1)
            end
        else
            if M_tip <= 1
                F1TD = 59 + 20 * log10(M_d_fan)
            else
                F1TD = 59 + 20 * log10(M_d_fan) - 20 * log10(M_tip / 1)
            end
        end

        # Rotor-stator spacing correction term (F2 of Eqn 12, Figure 12, modified by AlliedSignal):
        if settings.fan_id == false
            F2TD = -10 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan < 100
                F2TD = -10 * log10(RSS_fan / 300)
            else
                F2TD = -10 * log10(100 / 300)  # This is set to a constant 4.7712
            end
        end
    elseif settings.fan_RS_method == "geae"
        # Tip Mach-dependent term (F1 of Eqn 12 in report, Figure 10B, modified by GE):
        if M_d_fan <= 1
            if M_tip <= 1
                F1TD = 63
            else
                F1TD = 63 - 20 * log10(M_tip / 1)
            end
        else
            if M_tip <= 1
                F1TD = 63 + 20 * log10(M_d_fan)
            else
                F1TD = 63 + 20 * log10(M_d_fan) - 20 * log10(M_tip / 1)
            end
        end
        # Rotor-stator spacing correction term (F2 of Eqn 12, Figure 12, modified by GE):
        F2TD = -10 * log10(RSS_fan / 300)
    elseif settings.fan_RS_method == "kresja"
        # Tip Mach-dependent term (F1 of Eqn 12 in report, Figure 10B, modified by Krejsa):
        if M_d_fan <= 1
            F1TD = 46 - 20 * M_tip + 20 * log10(1 / 1.245)
        else
            F1TD = 46 - 20 * M_tip + 20 * log10(M_d_fan / 1.245)
        end
        # Rotor-stator spacing correction term (F2 of Eqn 12, Figure 12, modified by Krejsa):
        if settings.fan_id == false
            F2TD = -10 * log10(RSS_fan / 300)  # If no distortion
        else
            if RSS_fan < 100
                F2TD = -10 * log10(RSS_fan / 300)
            else
                F2TD = -10 * log10(100 / 300)  # This is set to a constant 4.7712
            end
        end
    else
        throw(DomainError("Invalid fan_RS_method specified. Specify: original / allied_signal / geae / kresja."))
    end

    # Theta correction term (F3 of Eqn 6, Figure 13B):
    if settings.fan_RS_method == "allied_signal"
        THT13B = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
        FIG13B = [-34, -30, -26, -22, -18, -14, -10.5, -6.5, -4, -1, 0, 0, 0, 0, -1, -3.5, -7, -11, -16]
        f_F3TD = LinearInterpolation(THT13B, FIG13B)
        F3TD = f_F3TD(theta)
    elseif settings.fan_RS_method == "kresja"
        THT13B = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
        FIG13B = [-50, -41, -33, -26, -20.6, -17.9, -14.7, -11.2, -9.3, -7.1, -4.7, -2, 0, 0.8, 1, -1.6, -4.2, -9,-15]
        f_F3TD = LinearInterpolation(THT13B, FIG13B)
        F3TD = f_F3TD(theta)
    elseif settings.fan_RS_method in ["original", "geae"]
        THT13B = [0, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190]
        FIG13B = [-39, -15, -11, -8, -5, -3, -1, 0, 0, -2, -5.5, -9, -13, -18, -13]
        f_F3TD = LinearInterpolation(THT13B, FIG13B)
        F3TD = f_F3TD(theta)
    else
        throw(DomainError("Invalid fan_RS_method specified. Specify: original / allied_signal / geae / kresja."))
    end

    # Added noise factor if there are inlet guide vanes:
    if settings.fan_igv == true
        CDT = 6 * ones(T, 1)
    else
        CDT = zeros(T, 1)
    end

    # Component value:
    tonlv_X = tsqem .+ F1TD .+ F2TD .+ F3TD .+ CDT

    return tonlv_X
end
function combination_tones(settings::PyObject, freq::Array{Float64, 1}, theta, M_tip, bpf, tsqem)

    # Combination tone (multiple pure tone or buzzsaw) calculations:
    # Note the original Heidmann reference states that MPTs should be computed if
    # the tangential tip speed is supersonic, but the ANOPP implementation states MPTs
    # should be computed if the relative tip speed is supersonic.  The ANOPP implementation
    # is used here, i.e., if M_tip >= 1.0, MPTs are computed.

    # Initialize solution matrices
    T = eltype(M_tip)
    dcp = zeros(T, settings.N_f)

    if M_tip >= 1
        if settings.fan_RS_method == "original"
            # Theta correction term (F2 of Eqn 8, Figure 16), original method:
            THT16 = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 180, 270]
            FIG16 = [-9.5, -8.5, -7, -5, -2, 0, 0, -3.5, -7.5, -9, -13.5, -9]
            f_F2CT = LinearInterpolation(THT16, FIG16)
            F2CT = f_F2CT(theta)

            # Spectrum slopes, original method:
            SL3 = [0, -30, -50, -30]
            SL4 = [0, 30, 50, 50]
            YINT3 = [0, -9.0309, -30.103, -27.0927]
            YINT4 = [0, 9.0309, 30.103, 45.1545]
        elseif settings.fan_RS_method == "allied_signal"
            # Theta correction term (F2 of Eqn 8, Figure 16), small fan method:
            THT16 = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 270]
            FIG16 = [-5.5, -4.5, -3, -1.5, 0, 0, 0, 0, -2.5, -5, -6, -6.9, -7.9, -8.8, -9.8, -10.7, -11.7, -12.6,-13.6, -6]
            f_F2CT = LinearInterpolation(THT16, FIG16)
            F2CT = f_F2CT(theta)
            # Spectrum slopes, small fan method:
            SL3 = [0, -15, -50, -30]
            SL4 = [0, 30, 50, 50]
            YINT3 = [0, -4.51545, -30.103, -27.0927]
            YINT4 = [0, 9.0309, 30.103, 45.1545]
        elseif settings.fan_RS_method == "geae"
            # Theta correction term (F2 of Eqn 8, Figure 16), large fan method:
            THT16 = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 180, 270]
            FIG16 = [-9.5, -8.5, -7, -5, -2, 0, 0, -3.5, -7.5, -9, -13.5, -9]
            f_F2CT = LinearInterpolation(THT16, FIG16)
            F2CT = f_F2CT(theta)
            # Spectrum slopes, GE large fan method:
            SL3 = [0, -30, -50, -30]
            SL4 = [0, 30, 50, 50]
            YINT3 = [0, -9.0309, -30.103, -27.0927]
            YINT4 = [0, 9.0309, 30.103, 45.1545]
        elseif settings.fan_RS_method == "kresja"
            # Theta correction term (F2 of Eqn 8, Figure 16), Krejsa method:
            THT16 = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190]
            FIG16 = [-28, -23, -18, -13, -8, -3, 0, -1.3, -2.6, -3.9, -5.2, -6.5, -7.9, -9.4, -11, -12.7, -14.5,-16.4, -18.4]
            f_F2CT = LinearInterpolation(THT16, FIG16)
            F2CT = f_F2CT(theta)
            # Spectrum slopes, Krejsa method:
            SL3 = [0, -20, -30, -20]
            SL4 = [0, 20, 30, 30]
            YINT3 = [0, -6.0206, -18.0618, -18.0618]
            YINT4 = [0, 6.0206, 18.0618, 27.0927]
        else
            throw(DomainError("Invalid fan_RS_method specified. Specify: original / allied_signal / geae / kresja."))
        end

        # Noise adjustment (reduction) if there are inlet guide vanes, for all methods:
        if settings.fan_igv == true
            CCT = -5 * ones(T, 1)
        else
            CCT = zeros(T, 1)
        end

        # Loop through the three sub-bpf terms:
        # K = 1; 1/2 bpf term
        # K = 2; 1/4 bpf term
        # K = 3; 1/8 bpf term
        for K in range(1, 3, step=1)
            # Tip Mach-dependent term (F1 of Eqn 8 in Heidmann report, Figure 15A):
            if settings.fan_RS_method == "original"
                # Original tip Mach number-dependent term of multiple pure tone noise
                FIG15MT = [0, 1.14, 1.25, 1.61]
                SL1 = [0, 785.68, 391.81, 199.2]
                SL2 = [0, -49.62, -50.06, -49.89]
                YINT1 = [0, 30, 30, 30]
                YINT2 = [0, 79.44, 83.57, 81.52]

                if M_tip < FIG15MT[K]
                    F1CT = SL1[K] * log10(M_tip) + YINT1[K]
                else
                    F1CT = SL2[K] * log10(M_tip) + YINT2[K]
                end
            elseif settings.fan_RS_method == "allied_signal"
                # Revised by AlliedSignal:  tip Mach number-dependent term of multiple pure tone noise.
                # Note that a 20log10 independent variable distribution is specified.
                FIG15MT = [0, 1.135, 1.135, 1.413]
                SL1 = [0, 5.9036, 5.54769, 2.43439]
                SL2 = [0, -0.632839, -0.632839, -1.030931]
                YINT1 = [0, 50, 50, 40]
                YINT2 = [0, 57.1896, 56.7981, 50.4058]

                if M_tip < FIG15MT[K]
                    F1CT = SL1[K] * 20 * log10(M_tip) + YINT1[K]
                else
                    F1CT = SL2[K] * 20 * log10(M_tip) + YINT2[K]
                end
            elseif settings.fan_RS_method == "geae"
                # Revised by GE: tip Mach number-dependent term of multiple pure tone noise
                FIG15MT = [0, 1.14, 1.25, 1.61]
                SL1 = [0, 746.8608, 398.3077, 118.9406]
                SL2 = [0, -278.96, -284.64, -43.52]
                YINT1 = [0, 30, 30, 36]
                YINT2 = [0, 88.37, 96.18, 69.6]

                if M_tip < FIG15MT[K]
                    F1CT = SL1[K] * log10(M_tip) + YINT1[K]
                else
                    F1CT = SL2[K] * log10(M_tip) + YINT2[K]
                end
            elseif settings.fan_RS_method == "kresja"
                FIG15MT = [0, 1.146, 1.322, 1.61]
                SL1 = [0, 785.68, 391.81, 199.2]
                SL2 = [0, -49.62, -50.06, -49.89]
                YINT1 = [0, -18, -15, -12]
                YINT2 = [0, 31.44, 38.57, 39.52]

                if M_tip < FIG15MT[K]
                    F1CT = SL1[K] * log10(M_tip) + YINT1[K]
                else
                    F1CT = SL2[K] * log10(M_tip) + YINT2[K]
                end
            else
                throw(DomainError("Invalid fan_RS_method specified. Specify: original / allied_signal / geae / kresja."))
            end

            CTLC = tsqem + F1CT + F2CT + CCT

            # Frequency-dependent term (F3 of Eqn 9, Figure 14):
            FK = 2^K

            # Cycle through frequencies and make assignments:
            for j in range(1, settings.N_f, step=1)
                FQFB = freq[j] / bpf

                if FQFB <= 1 / FK
                    # For frequencies less than the subharmonic:
                    F3CT = SL4[K] * log10(FQFB) + YINT4[K]
                else
                    # For frequencies greater than the subharmonic:
                    F3CT = SL3[K] * log10(FQFB) + YINT3[K]
                end
                # Be sure to add the three sub-bpf components together at each frequency:
                dcp[j] = dcp[j] + 10^(0.1 * (CTLC + F3CT))
            end
        end
    end

    return dcp
end
function calculate_harmonics(settings::PyObject, freq::Array{Float64, 1}, theta, tonlv_I, tonlv_X, i_cut, M_tip, bpf, comp::String)

    # Assign discrete interaction tones at bpf and harmonics to proper bins (see figures 8 and 9):
    # Initialize solution matrices
    T = eltype(M_tip)
    dp = zeros(T, settings.N_f)
    dpx = zeros(T, settings.N_f)

    nfi = 1
    for ih in range(1, settings.n_harmonics, step=1)

        # Determine the tone fall-off rates per harmonic (harm_i and harm_x):
        if settings.fan_RS_method == "original"
            if settings.fan_igv == false
                # For fans without inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    else
                        # For cut-on fans, harmonics:
                        harm_i = 3 * (ih - 1) * ones(T, 1)
                        harm_x = 3 * (ih - 1) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    else
                        # For cut-off fans, harmonics:
                        harm_i = 3 * (ih - 1) * ones(T, 1)
                        harm_x = 3 * (ih - 1) * ones(T, 1)
                    end
                end
            elseif settings.fan_igv == true
                # For fans with inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    else
                        # For cut-on fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    else
                        # For cut-off fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                end
            end
        elseif settings.fan_RS_method == "allied_signal"
            if settings.fan_igv == false
                # For fans without inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    elseif ih == 2
                        # For cut-on fans, second harmonic:
                        harm_i = 9.2 * ones(T, 1)
                        harm_x = 9.2 * ones(T, 1)
                    else
                        # For cut-on fans, upper harmonics:
                        harm_i = (3 * ih + 1.8) * ones(T, 1)
                        harm_x = (3 * ih + 1.8) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    elseif ih == 2
                        # For cut-off fans, second harmonic:
                        harm_i = 9.2 * ones(T, 1)
                        harm_x = 9.2 * ones(T, 1)
                    else
                        # For cut-off fans, upper harmonics:
                        harm_i = (3 * ih + 1.8) * ones(T, 1)
                        harm_x = (3 * ih + 1.8) * ones(T, 1)
                    end
                end
            elseif settings.fan_igv == true
                # For fans with inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    else
                        # For cut-on fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    else
                        # For cut-off fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                end
            end
        elseif settings.fan_RS_method == "geae"
            if settings.fan_igv == false
                # For fans without inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    else
                        # For cut-on fans, harmonics:
                        if M_tip < 1.15
                            harm_i = 6 * (ih - 1) * ones(T, 1)
                        else
                            harm_i = 9 * (ih - 1) * ones(T, 1)
                        end
                        harm_x = 3 * (ih - 1) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    else
                        # For cut-off fans, harmonics:
                        if M_tip < 1.15
                            harm_i = 6 * (ih - 1) * ones(T, 1)
                        else
                            harm_i = 9 * (ih - 1) * ones(T, 1)
                        harm_x = 3 * (ih - 1) * ones(T, 1)
                        end
                    end
                end
            elseif settings.fan_igv == true
                # For fans with inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    else
                        # For cut-on fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    else
                        # For cut-off fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                end
            end
        elseif settings.fan_RS_method == "kresha"
            if settings.fan_igv == false
                # For fans without inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    else
                        # For cut-on fans, harmonics:
                        harm_i = 3 * (ih - 1) * ones(T, 1)
                        harm_x = 3 * (ih - 1) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    else
                        # For cut-off fans, harmonics:
                        harm_i = 3 * (ih - 1) * ones(T, 1)
                        harm_x = 3 * (ih - 1) * ones(T, 1)
                    end
                end
            elseif settings.fan_igv == true
                # For fans with inlet guide vanes:
                if i_cut == [0]
                    # For cut-on fans, fundamental:
                    if ih == 1
                        harm_i = zeros(T, 1)
                        harm_x = zeros(T, 1)
                    else
                        # For cut-on fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                elseif i_cut == [1]
                    # For cut-off fans, fundamental:
                    if ih == 1
                        harm_i = 8 * ones(T, 1)
                        harm_x = 8 * ones(T, 1)
                    else
                        # For cut-off fans, harmonics:
                        harm_i = 3 * (ih + 1) * ones(T, 1)
                        harm_x = 3 * (ih + 1) * ones(T, 1)
                    end
                end
            end
        else
            throw(DomainError("Invalid fan_RS_method specified. Specify: original / allied_signal / geae / kresja."))
        end
        # Calculate TCS and distor
        if (settings.fan_id == true) && (settings.fan_RS_method != "geae")
            # Assign the increment to the fundamental tone along with a 10 dB per harmonic order fall-off
            # for cases with inlet flow distortion (see figure 9):
            distor = 10^(0.1 * tonlv_I - ih + 1)
            TCS = zeros(T, 1)
        elseif (settings.fan_id == true) && (settings.fan_RS_method == "geae")
            # Compute suppression factors for GE#s "Flight cleanup Turbulent Control Structure."
            # Approach or takeoff values to be applied to inlet discrete interaction tones
            # at bpf and 2bpf.  Accounts for observed in-flight tendencies.
            TCSTHA = [0, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
            TCSAT1 = [0, 4.8, 4.8, 5.5, 5.5, 5.3, 5.3, 5.1, 4.4, 3.9, 2.6, 2.3, 1.8, 2.1, 1.7, 1.7, 2.6, 3.5, 3.5,3.5]
            TCSAT2 = [0, 5.8, 5.8, 3.8, 5.3, 6.4, 3.5, 3, 2.1, 2.1, 1.1, 1.4, 0.9, 0.7, 0.7, 0.4, 0.6, 0.8, 0.8,0.8]
            TCSAA1 = [0, 5.6, 5.6, 5.8, 4.7, 4.6, 4.9, 5.1, 2.9, 3.2, 1.6, 1.6, 1.8, 2.1, 2.4, 2.2, 2, 2.8, 2.8,2.8]
            TCSAA2 = [0, 5.4, 5.4, 4.3, 3.4, 4.1, 2, 2.9, 1.6, 1.3, 1.5, 1.1, 1.4, 1.5, 1, 1.8, 1.6, 1.6, 1.6, 1.6]

            if settings.ge_flight_cleanup == "takeoff"
                # Apply takeoff values:
                if ih == 1
                    f_TCS = LinearInterpolation(TCSTHA, TCSAT1)
                    TCS = f_TCS(theta)
                elseif ih == 2
                    f_TCS = LinearInterpolation(TCSTHA, TCSAT2)
                    TCS = f_TCS(theta)
                else
                    TCS = zeros(T, 1)
                end
            elseif settings.ge_flight_cleanup == "approach"
                # Apply approach values:
                if ih == 1
                    f_TCS = LinearInterpolation(TCSTHA, TCSAA1)
                    TCS = f_TCS(theta)
                elseif ih == 2
                    f_TCS = LinearInterpolation(TCSTHA, TCSAA2)
                    TCS = f_TCS(theta)
                else
                    TCS = zeros(T, 1)
                end
            elseif settings.ge_flight_cleanup == "none"
                # Apply zero values (i.e., fan inlet flow is distorted):
                TCS = zeros(T, 1)
            end
                # Inlet distortion effects are always included in basic inlet tone model of the GE method.
            # Flight cleanup levels are then subtracted from the inlet tones if the flow is not distorted.
            # The flight cleanup levels are set to zero if the flow is distorted.
            # Use the same increment as the original method and the same 10 dB per harmonic fall-off rate.
            distor = 10 .^(0.1 * (tonlv_I .- TCS) - ih + 1)
        else
            distor = zeros(T, 1)
            TCS = zeros(T, 1)
        end

        # Calculate tone power
        if comp == "fan_inlet" # or comp == "inlet RS":
            tonpwr_i = 10 .^(0.1 * (tonlv_I .- harm_i .- TCS)) .+ distor
        else
            tonpwr_i = zeros(T, 1)
        end
        if comp == "fan_discharge" # or comp == "discharge RS" or comp == "total":
            tonpwr_x = 10 .^(0.1 * (tonlv_X .- harm_x))
        else
            tonpwr_x = zeros(T, 1)
        end

        # Compute filter bandwidths:
        filbw = 1  # Fraction of filter bandwidth with gain of unity (default to unity)
        F1 = 0.78250188 + 0.10874906 * filbw
        F2 = 1 - 0.10874906 * filbw
        F3 = 1 + 0.12201845 * filbw
        F4 = 1.2440369 - 0.12201845 * filbw

        # Cycle through frequencies and assign tones to 1/3rd octave bins:
        ll = NaN
        for l in range(nfi, settings.N_f, step=1)
            Frat = bpf * ih / freq[l]
            FR = 1
            if Frat .< F1
                break
            elseif Frat .> F4
                ll = l
                continue
            #elseif Frat .> F3
                #FR = (F4 .- Frat) ./ (F4 .- F3)
            #elseif Frat .< F2
                #FR = (Frat .- F1) ./ (F2 .- F1)
            end

            dp[l] = dp[l] .+ tonpwr_i[1] .* FR
            dpx[l] = dpx[l] .+ tonpwr_x[1] .* FR
            nfi = ll

            #continue
        end
    end

    return dp, dpx
end
function calculate_cutoff(M_tip_tan, B_fan::Int64, V_fan::Int64)

    # Vane/blade ratio parameter:
    vane_blade_ratio = 1 - V_fan / B_fan
    if vane_blade_ratio == 0
        vane_blade_ratio = 1e-6
    end

    # Cutoff parameter:
    cutoff = abs(M_tip_tan / vane_blade_ratio)

    # if the cutoff parameter is less than 1.05 and the tip Mach is less than
    # unity, the fan is cut off (i.e., the tones are reduced in magnitude):
    if (cutoff < 1.05) && (M_tip_tan < 1)
        i_cut = ones(eltype(cutoff), (1, ))
    elseif (cutoff < 1.05) && (M_tip_tan >= 1)
        i_cut = zeros(eltype(cutoff), (1, ))
    elseif (cutoff >= 1.05) && (M_tip_tan < 1)
        i_cut = zeros(eltype(cutoff), (1, ))
    else
        i_cut = zeros(eltype(cutoff), (1, ))
    end

    return i_cut
end

function fan(settings::PyObject, data::PyObject, ac::PyObject, n_t::Int64, shielding::Array{Float64, 2}, M_0, c_0, T_0, rho_0, theta, DTt_f_star, mdot_f_star, N_f_star, A_f_star, d_f_star, comp::String)

    # Initialize solution
    T = eltype(DTt_f_star)
    msap_fan = zeros(T, (n_t, settings.N_f))

    for i in range(1, n_t, step=1)

        ### Extract the inputs
        delta_T_fan = DTt_f_star[i] * T_0[i]  # Total temperature rise across fan [R]
        rpm = N_f_star[i] * 60 * c_0[i] / (d_f_star[i] * sqrt(settings.A_e))  # Shaft speed [rpm]
        M_tip_tan = (d_f_star[i] * sqrt(settings.A_e) / 2) * rpm * 2 * π / 60 / c_0[i]  # Tangential (i.e., radius*omega) tip Mach number: ! Doesn"t need correction
        mdot_fan = mdot_f_star[i] * rho_0[i] * c_0[i] * settings.A_e  # Airflow [kg/s]
        bpf = rpm * ac.B_fan / 60. / (1 .- M_0[i] * cos.(theta[i] * π / 180))  # Blade passing frequency, [Hz]
        flow_M = mdot_fan / (rho_0[i] .* A_f_star[i] * settings.A_e * c_0[i])  # Fan face flow Mach number (assumes ambient and fan face static densities are equal): !!!!!!!
        M_tip = (M_tip_tan^2 + flow_M^2)^0.5  # Relative (i.e., helical) tip Mach number: ! Doesn"t need correction

        # Temperature-flow power base term:
        rho_sl = 1.22514
        c_sl = 340.29395
        if settings.fan_BB_method == "kresja"
            tsqem = 40 * log10(delta_T_fan * 1.8) + 10 * log10(2.20462 * mdot_fan / (1 - M_0[i] * cos(theta[i] * π / 180)).^4) - 20 * log10(settings.r_0) - 10*log10(rho_sl^2 * c_sl^4)
        else  # All other methods:
            tsqem = 20 * log10(delta_T_fan * 1.8) + 10 * log10(2.20462 * mdot_fan / (1 - M_0[i] * cos(theta[i] * π / 180)).^4) - 20 * log10(settings.r_0) - 10*log10(rho_sl^2 * c_sl^4)
        end

        # Calculate individual noise components
        if comp == "fan_inlet"
            spl_IB = inlet_broadband(settings, theta[i], M_tip, tsqem, ac.M_d_fan, ac.RSS_fan)
            spl_DB = zeros(T, 1)
            tonlv_I = inlet_tones(settings, theta[i], M_tip, tsqem, ac.M_d_fan, ac.RSS_fan)
            tonlv_X = zeros(T, 1)
        elseif comp == "fan_discharge"
            spl_IB = zeros(T, 1)
            spl_DB = discharge_broadband(settings, theta[i], M_tip, tsqem, ac.M_d_fan, ac.RSS_fan)
            tonlv_I = zeros(T, 1)
            tonlv_X = discharge_tones(settings, theta[i], M_tip, tsqem, ac.M_d_fan, ac.RSS_fan)
        end

        if settings.combination_tones
            dcp = combination_tones(settings,  data.f, theta[i], M_tip, bpf, tsqem)
        else
            dcp = zeros(T, settings.N_f)
        end

        # Calculate if cut-off happens (1) or not (0)
        i_cut = calculate_cutoff(M_tip_tan, ac.B_fan, ac.V_fan)

        # Assign tones_to bands
        dp, dpx = calculate_harmonics(settings, data.f, theta[i], tonlv_I, tonlv_X, i_cut, M_tip, bpf, comp)

        # Final calculations;  cycle through frequencies and assign values:
        if settings.fan_BB_method == "allied_signal"
            # Eqn 2 or Figure 3A:
            # if data.f[j] / bpf < 2:
            # FLOGinlet,exit = -10 * np.log10(np.exp(-0.35 * (np.log(data.f[j] / bpf / 2.0) / np.log(2.2))^2))
            FLOGinlet = 2.445096095 * (log.(data.f / bpf / 2)).^2
            # elif data.f[j] / bpf > 2:
            # FLOGinlet,exit = -10 * np.log10(np.exp(-2.0 * (np.log(data.f[j] / bpf / 2.0) / np.log(2.2))^2))
            FLOGinlet[data.f/bpf .> 2] = (13.97197769 * (log.(data.f / bpf / 2)).^2)[data.f/bpf .> 2]
            FLOGexit = FLOGinlet

        elseif settings.fan_BB_method == "kresja"
            # Eqn 2 or Figure 3A:
            # FLOGinlet = -10 * np.log10(np.exp(-0.5 * (np.log(data.f[j] / bpf / 4) / np.log(2.2))^2))
            # Which may be simplified as:
            FLOGinlet = 3.4929944 * (log.(data.f / bpf / 4)).^2
            # FLOGexit = -10 * np.log10(np.exp(-0.5 * (np.log(data.f[j] / bpf / 2.5) / np.log(2.2))^2))
            # Which may be simplified as:
            FLOGexit = 3.4929944 * (log.(data.f / bpf / 2.5)).^2
        else
            # For the original or the GE large fan methods:
            # Eqn 2 or Figure 3A:
            # FLOGinlet,exit = -10 * np.log10(np.exp(-0.5 * (np.log(data.f[j] / bpf / 2.5) / np.log(2.2))^2))
            # Which may be simplified as:
            FLOGinlet = 3.4929944 * (log.(data.f ./ bpf / 2.5)).^2
            FLOGexit = FLOGinlet
        end

        if comp == "fan_inlet" # or comp == "inlet BB":
            PLEV = 10 .^(0.1 * (spl_IB .- FLOGinlet))
            PLEV = PLEV .+ dp

        # Add discrete tone and broadband components for exhaust noise:
        elseif comp == "fan_discharge" # or comp == "discharge BB":
            PLEV = 10 .^(0.1 * (spl_DB .- FLOGexit))
            PLEV = PLEV .+ dpx
        else
            throw(DomainError("Invalid component specified."))
        end

        # Add inlet combination tones if needed:
        if (M_tip > 1) && (settings.combination_tones == true)
            PLEV = PLEV .+ dcp
        end

        # Add ambient correction
        msap_j = PLEV * settings.p_ref^2
        for i in range(1, settings.N_f, step=1)
            if PLEV[i] == 0.
                msap_j[i] = zeros(T, 1)
            end
        end

        # Multiply for number of engines
        msap_j = msap_j * ac.n_eng

        # Fan liner suppression
        if settings.fan_liner_suppression
            if comp == "fan_inlet"
                data_angles = data.supp_fi_angles
                data_freq = data.supp_fi_freq
                data_supp = data.supp_fi
            elseif comp == "fan_discharge"
                data_angles = data.supp_fd_angles
                data_freq = data.supp_fd_freq
                data_supp = data.supp_fd
            else
                throw(DomainError("Invalid fan component."))
            end
            f_supp = LinearInterpolation((data_freq, data_angles), data_supp)
            supp = f_supp.(data.f, theta[i].*ones(T, settings.N_f, ))
            msap_j = supp .* msap_j
        end

        # Fan inlet shielding
        if (settings.shielding==true) && (comp == "fan_inlet")
            msap_j = msap_j ./ (10 .^(shielding[i, :] ./ 10.))
        end

        # Normalize msap by reference pressure
        msap_fan[i,:] = clamp.(msap_j, 1e-99, Inf)/settings.p_ref^2
    end

    return msap_fan
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
module Reconstruction

export inverse_compress, inverse_digitize, quantize

function inverse_compress(pieces::Array{Float64,}, start=0.0::AbstractFloat)
    """
    Reconstruct time series from its first value `start` and its `pieces`.
    `pieces` must have (at least) two columns, incremenent and window width, resp.
    A window width w means that the piece ranges from s to s+w.
    In particular, a window width of 1 is allowed.

    Parameters
    ----------
    start -
        First element of original time series. Applies vertical shift in
        reconstruction.

    pieces -
        Array with three columns, each row contains increment, length,
        error for the segment. Only the first two columns are required.

    Returns
    -------
    time_series : Reconstructed time series
    """

    r_time_series = [start];
    num_pieces = size(pieces)[1];
    for j in range(1, stop=num_pieces)
        x = range(0, stop=pieces[j,1]) / (pieces[j,1]) * pieces[j,2];
        y = broadcast(+, r_time_series[length(r_time_series)], x);
        r_time_series = vcat(r_time_series, y[2:length(y)]);
    end
    return r_time_series
end

function inverse_digitize(symbols::Array{Char,1}, model)
    """
    Convert symbolic representation back to compressed representation for reconstruction.

    Parameters
    ----------
    symbols -
        Time series in symbolic representation using unicode characters starting
        with character 'a'.

    centers -
        Centers of clusters from clustering algorithm. Each centre corresponds
        to character in string.

    Returns
    -------
    pieces -
        Time series in compressed format.
    """

    pieces = zeros((0, 2));
    for symbol in symbols
        pc = model.centers[model.hashmap[symbol],:];
        pieces = vcat(pieces, pc');
    end
    return pieces;
end

function quantize(pieces::Array{Float64,})
    """
    Realign window lengths with integer grid.

    Parameters
    ----------
    pieces: Time series in compressed representation.

    Returns
    -------
    pieces: Time series in compressed representation with window length adjusted to integer grid.
    """
    num_pieces = size(pieces)[1];
    if num_pieces == 1
        pieces[1,1] = round(pieces[1,1], digits=0);
    else
        for p in range(1, stop=num_pieces-1)
            corr = round(pieces[p,1]) - pieces[p,1];
            pieces[p,1] = round(pieces[p,1] + corr);
            pieces[p+1,1] = pieces[p+1,1] - corr;
            if pieces[p,1] == 0
                pieces[p,1] = 1;
                pieces[p+1,1] = pieces[p+1,1] - 1;
            end
        end

        pieces[num_pieces,1] = round(pieces[num_pieces,1],digits=0);
    end
    return pieces;
end

end

ext: jl

max_stars_repo_name: nla-group/Julia-ABBA

max_issues_repo_name: nla-group/Julia-ABBA

max_forks_repo_name: nla-group/Julia-ABBA


text:
using ChainRulesCore
import ChainRulesCore: rrule
import .Code

function rrule(::typeof(Code.create_array), A::Type{<:AbstractArray}, T, u::Val{j}, d::Val{dims}, elems...) where {dims, j}
  y = Code.create_array(A, T, u, d, elems...)
  function create_array_pullback(Δ)
    dx = Δ
    (NO_FIELDS, DoesNotExist(), DoesNotExist(), DoesNotExist(), DoesNotExist(), dx..., ntuple(_ -> DoesNotExist(), length(elems) - prod(dims) + j)...)
  end
  y, create_array_pullback
end


text:
using BEAST
using CompScienceMeshes

using Test

Γ = meshrectangle(1.0,1.0,1.0)
γ = meshsegment(1.0,1.0,3)
X = raviartthomas(Γ, γ)

x = divergence(X)

@test numfunctions(X) == 2
@test numfunctions(x) == numfunctions(X)
for (_f,_g) in zip(x.fns, X.fns)
    @test length(_f) == length(_g)
    if length(_f) == 1
        c = _f[1].cellid
        _s = chart(Γ, Γ.faces[c])
        @test _f[1].coeff == 1 / volume(_s)
    else
        @test _f[1].coeff + _f[2].coeff == 0
    end
end

edges = skeleton(Γ, 1)
vps = cellpairs(Γ, edges)
rwg = raviartthomas(Γ, vps)
dvg = divergence(rwg)

D = spzeros(Int,numcells(edges), numcells(Γ))
for (_m,_f) in enumerate(dvg.fns)
    for _s in _f
        _e = _s.cellid
        D[_m,_e] = _s.coeff
    end
end

C = connectivity(edges, Γ)
@test size(C) == size(D')

## Test the curl of the linear Lagrange elements
using CompScienceMeshes
using BEAST
using Test

m = meshrectangle(1.0, 1.0, 0.25)
X = duallagrangec0d1(m)
@test numfunctions(X) == numcells(m)

Y = curl(X)
@test numfunctions(X) == numfunctions(Y)
@test isa(Y, BEAST.RTBasis)
@test isa(refspace(Y), BEAST.BEAST.RTRefSpace)

# test the chain property
Z = divergence(Y)
for fn ∈ Z.fns
    total = zero(scalartype(Z))
    for _sh ∈ fn
        total += _sh.coeff
    end
    @test total ≈ 0
end


text:
using Clang
using CImGui.LibCImGui.CImGui_jll

const HELPER_H = joinpath(dirname(CImGui_jll.libcimgui_path), "..", "include", "helper.h") |> normpath

# create a work context
ctx = DefaultContext()

# parse headers
parse_headers!(ctx, [HELPER_H], args=map(x->"-I"*x, find_std_headers()), includes=[LLVM_INCLUDE])

# settings
ctx.libname = "libcimgui_helper"
ctx.options["is_function_strictly_typed"] = false
ctx.options["is_struct_mutable"] = false

# write output
api_file = joinpath(@__DIR__, "libhelper_api.jl")
api_stream = open(api_file, "w")

for trans_unit in ctx.trans_units
    root_cursor = getcursor(trans_unit)
    push!(ctx.cursor_stack, root_cursor)
    header = spelling(root_cursor)
    @info "wrapping header: $header ..."
    # loop over all of the child cursors and wrap them, if appropriate.
    ctx.children = children(root_cursor)
    for (i, child) in enumerate(ctx.children)
        child_name = name(child)
        child_header = filename(child)
        ctx.children_index = i
        # choose which cursor to wrap
        startswith(child_name, "__") && continue  # skip compiler definitions
        child_name in keys(ctx.common_buffer) && continue  # already wrapped
        child_header != header && continue  # skip if cursor filename is not in the headers to be wrapped

        wrap!(ctx, child)
    end
    @info "writing $(api_file)"
    println(api_stream, "# Julia wrapper for header: $(basename(header))")
    println(api_stream, "# Automatically generated using Clang.jl\n")
    print_buffer(api_stream, ctx.api_buffer)
    empty!(ctx.api_buffer)  # clean up api_buffer for the next header
end
close(api_stream)


text:
using DSGE, ModelConstructors
using Nullables, DataFrames, OrderedCollections, Dates

#############################
# To use:
# Just run in the Julia REPL
# include("run_default.jl")
#############################

##############
# Model Setup
##############
# Instantiate the FRBNY DSGE model object
m = Model1002()

# estimate as of 2015-Q3 using the default data vintage from 2015 Nov 27
m <= Setting(:data_vintage, "181115")
m <= Setting(:date_forecast_start, quartertodate("2018-Q4"))

# The following settings ensure that this script runs in
# a short amount of time. To properly estimate and
# forecast, we recommend either using the default settings
# (i.e. comment out the settings below) or
# changing the settings yourself.
m <= Setting(:n_mh_simulations, 100) # Do 100 MH steps during estimation
m <= Setting(:use_population_forecast, false) # Population forecast not available as data to turn off
m <= Setting(:forecast_block_size, 5) # adjust block size to run on small number of estimations

#############
# Estimation
#############
# Reoptimize parameter vector, compute Hessian at mode, and full posterior
# parameter sampling.
#
# Note some columns will have missing data because not all our data
# is publicly available. Warnings will be thrown after calling load_data as a result.
# However, estimate will still run when data is missing.
df = load_data(m, try_disk = false, check_empty_columns = false, summary_statistics = :none)
data = df_to_matrix(m, df)
estimate(m, data)

# produce LaTeX tables of parameter moments
moment_tables(m)

############
# Forecasts
############
# (1) forecast_one produces forecasts of the variables specified in output_vars
# and stores the forecasts in model units (log deviations from steady state).
# (2) compute_meansbands stores the forecast in a MeansBands type, defined in DSGE.jl,
# which transforms the forecast from model units to "human readable" units, i.e.
# the GDP forecast is transformed from annualized quarterly log per-capita growth rates
# to annualized quarterly aggregate percent change.

# The variables we want to forecast. In this case, all of the model observables
output_vars = [:histobs, :forecastobs]

# Modal forecast (point forecast)
forecast_one(m, :mode, :none, output_vars; check_empty_columns = false)
compute_meansbands(m, :mode, :none, output_vars; check_empty_columns = false)

# Full-distribution forecast (point forecast (mean) and uncertainty bands)
# Optionally add 10 processes to run the forecast in parallel (uncomment the 3 lines below).
# Alternatively, you can load the ClusterManagers package and add processes
# using one of the schedulers such as SGE or Slurm.
# addprocs(10)
# @everywhere using DSGE
# m <= Setting(:use_parallel_workers, true)

forecast_one(m, :full, :none, output_vars; check_empty_columns = false)
compute_meansbands(m, :full, :none, output_vars; check_empty_columns = false)

# Comment out the line below if you did not run the forecast in parallel.
# rmprocs(procs())


text:
module SimpleNLUtils

import SimpleNL: Model, variable, Expression, Expression1, Expression2, Gradient, SparseJacobian, ExpressionSum, Variable, Constant, Parameter, Sink, Field, Field1, num_variables, num_constraints, Constraint, inner, index, fill_sparsity!, non_caching_eval, DUMMY
import SparseArrays: SparseMatrixCSC, sparse
import LinearAlgebra: mul!, norm

export get_terms, get_entries, sparsity, set_KKT_error_evaluator!

include("utils.jl")
include("structured.jl")
include("kkt.jl")

end # module


text:
reload("Missings"); reload("WeakRefStrings"); reload("DataStreams"); reload("CSV"); reload("DataStreamsIntegrationTests")
# NamedTuples
FILE = joinpath(DataStreamsIntegrationTests.DSTESTDIR, "randoms_small.csv")
DF = CSV.read(FILE)
DF2 = CSV.read(FILE)
dfsource = DataStreamsIntegrationTests.Tester("NamedTuple", x->x, false, NamedTuple, (DF,), DataStreamsIntegrationTests.scalartransforms, DataStreamsIntegrationTests.vectortransforms, x->x, x->nothing)
dfsink = DataStreamsIntegrationTests.Tester("NamedTuple", x->x, false, NamedTuple, (DF2,), DataStreamsIntegrationTests.scalartransforms, DataStreamsIntegrationTests.vectortransforms, x->x, x->nothing)

# CSV
FILE2 = joinpath(DataStreamsIntegrationTests.DSTESTDIR, "randoms2_small.csv")
csvsource = DataStreamsIntegrationTests.Tester("CSV.Source", CSV.read, true, CSV.Source, (FILE,), DataStreamsIntegrationTests.scalartransforms, DataStreamsIntegrationTests.vectortransforms, x->x, x->nothing)
csvsink = DataStreamsIntegrationTests.Tester("CSV.Sink", CSV.write, true, CSV.Sink, (FILE2,), DataStreamsIntegrationTests.scalartransforms, DataStreamsIntegrationTests.vectortransforms, x->CSV.read(FILE2; use_mmap=false), x->rm(FILE2))

DataStreamsIntegrationTests.teststream([dfsource, csvsource], [dfsink]; rows=99)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkConfigurationDiagnosticResult <: SwaggerModel
    profile::Any # spec type: Union{ Nothing, NetworkConfigurationDiagnosticProfile } # spec name: profile
    networkSecurityGroupResult::Any # spec type: Union{ Nothing, NetworkSecurityGroupResult } # spec name: networkSecurityGroupResult

    function NetworkConfigurationDiagnosticResult(;profile=nothing, networkSecurityGroupResult=nothing)
        o = new()
        validate_property(NetworkConfigurationDiagnosticResult, Symbol("profile"), profile)
        setfield!(o, Symbol("profile"), profile)
        validate_property(NetworkConfigurationDiagnosticResult, Symbol("networkSecurityGroupResult"), networkSecurityGroupResult)
        setfield!(o, Symbol("networkSecurityGroupResult"), networkSecurityGroupResult)
        o
    end
end # type NetworkConfigurationDiagnosticResult

const _property_map_NetworkConfigurationDiagnosticResult = Dict{Symbol,Symbol}(Symbol("profile")=>Symbol("profile"), Symbol("networkSecurityGroupResult")=>Symbol("networkSecurityGroupResult"))
const _property_types_NetworkConfigurationDiagnosticResult = Dict{Symbol,String}(Symbol("profile")=>"NetworkConfigurationDiagnosticProfile", Symbol("networkSecurityGroupResult")=>"NetworkSecurityGroupResult")
Base.propertynames(::Type{ NetworkConfigurationDiagnosticResult }) = collect(keys(_property_map_NetworkConfigurationDiagnosticResult))
Swagger.property_type(::Type{ NetworkConfigurationDiagnosticResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkConfigurationDiagnosticResult[name]))}
Swagger.field_name(::Type{ NetworkConfigurationDiagnosticResult }, property_name::Symbol) =  _property_map_NetworkConfigurationDiagnosticResult[property_name]

function check_required(o::NetworkConfigurationDiagnosticResult)
    true
end

function validate_property(::Type{ NetworkConfigurationDiagnosticResult }, name::Symbol, val)
end


text:
module test020
using Speckles

function run1()
    Speckles.run_020()
end
export run1

end

import .test020

test020.run1()


text:
include(string(pwd(),"/github/DSGE-HMC/log_like.jl"))

data = randn(100,3)
A = [0.5 0.3 0.2;0.4 0.2 -0.1; -0.5 0.8 -0.9]

loglike(A,data)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
## run "../flow_embedding.py" first to generate data
using Plots
using SparseArrays
using LinearAlgebra
using TextParse
using DelimitedFiles
using PyCall
using ColorSchemes

function read_graph(filename)
    src, dst, val = csvread(filename;
        spacedelim=true, colparsers=[Int,Int,Float64], header_exists=false)[1]
    src .+= 1
    dst .+= 1
    n = maximum(src)
    n = max(maximum(dst),n)
    A = sparse(src, dst, val, n, n)
    A = max.(A,A')
    return findnz(triu(A,1))[1:2]
end

src, dst = read_graph("../../dataset/lawlor-spectra-k32.edgelist.gz")
# x, y = csvread("$../../dataset/lawlor-spectra-k32.edgelist";
#         spacedelim=true, colparsers=[Float64,Float64],
#         header_exists=false)[1]
V = readdlm("../../dataset/lawlor-spectra-k32.coords")
V[:,1] *= -10
V[:,2] *= -4

function plotperm(v)
  return invperm(sortperm(v))
end
# function myscatter(x,y,color;kwargs...)
#   p = sortperm(color)
#   scatter(x[p],y[p], marker_z = color[p],
#     label="", alpha=0.1, markerstrokewidth=0,
#     colorbar=false, framestyle=:none; kwargs...)
# end

##
pfile = "embeddings/spectral_3_bfs.p.gz"

pickle = pyimport("pickle")
gzip = pyimport("gzip")
f = gzip.open(pfile)
data= pickle.load(f)
f.close()
function _pycsr_to_sparse(rowptr, colinds, vals, m, n)
  rowinds = zeros(eltype(colinds),0)
  for i in 1:length(rowptr)-1
    for j in rowptr[i]:rowptr[i+1]-1
      push!(rowinds, i)
    end
  end
  return sparse(rowinds, colinds, vals, m, n)
end
function pycsr_to_sparse(pycsr)
  rowptr = pycsr.indptr
  colinds = pycsr.indices
  vals = pycsr.data
  colinds .+= 1
  m, n = data.shape
  # put in a fucntion so we can do type inference
  return _pycsr_to_sparse(rowptr, colinds, vals, m, n )
end

X = pycsr_to_sparse(data)

flipexp = findall(X[517173,:] .> 0) # this is at the opposite end
for f in flipexp
  X[:,f] = 1.0 .-X[:,f]
end

X = Matrix(X) # convert to Matrix
X ./= 10 # convert to 0, 1
subset = vec((sum(X,dims=2) .> 0))
subsetnodes = findall(subset)
subset_c = vec((sum(X,dims=2) .<= 0))
subsetnodes_c = findall(subset_c)
U,sig,Vt = svd(X.-1)

new_coords = zeros(size(X,1),2)
new_coords[subsetnodes,1] = plotperm(-1 * U[subsetnodes,1])
new_coords[subsetnodes,2] = plotperm(-1 * U[subsetnodes,2])
i = argmax(new_coords[subsetnodes,1])
new_coords[subsetnodes_c,1] .= new_coords[subsetnodes[i],1] + 10000
new_coords[subsetnodes_c,2] .= new_coords[subsetnodes[i],2] + 10000

using Clustering

Xs = transpose(new_coords[subsetnodes,:])

R = kmeans(Xs,100; maxiter=200, display=:iter)

clusters = Dict(i=>[] for i = 1:100)

for i = 1:length(subsetnodes)
    push!(clusters[R.assignments[i]],subsetnodes[i]-1)
end

pushfirst!(PyVector(pyimport("sys")."path"), "/homes/liu1740/Research/LocalGraphClustering/")
lgc = pyimport("localgraphclustering")

G = lgc.GraphLocal("../../dataset/lawlor-spectra-k32.edgelist","edgelist")
conds = []
for i = 1:100
    push!(conds,G.compute_conductance(clusters[i]))
end

max_stars_repo_name: dgleich/flowpaper-code

max_issues_repo_name: dgleich/flowpaper-code

max_forks_repo_name: dgleich/flowpaper-code


text:
######################################################################## # Codes for the paper:
# An Anthropocene Map of Genetic Diversity
#
# Andreia Miraldo, Sen Li, Michael K. Borregaard, Alexander Floréz-Rodriguéz, Shyam Gopalakrishnan, Mirneza Risvanovic, Zhiheng Wang, Carsten Rahbek, Katharine A. Marske & David Nogués-Bravo
#
# Submitted to Science, 2016
# Code in this file by Michael K. Borregaard
#########################################################################

#Load the necessary library
using DataFrames
include("GD_summary_functions.jl")

"""
Do a rarefaction where all grid cell / species combinations with > N sequences get resampled to exactly N sequences
**Parameters**
* 'data' : a matrix generated by create_master_matrices()
* 'N'    : the number of sequences to rarefy down to
"""
function rare(data::DataFrame, N::Int)
    subset = by(dats, [:species, :cell]) do dat
        seqs = unique(vcat(dat[:seq1], dat[:seq2]))
        if size(seqs,1) < N
            return(dat[1:0,:])
        end
        randseqs = sample(seqs, N, replace = false)
        dat = dat[findin(dat[:seq1], randseqs),:]
        dat = dat[findin(dat[:seq2], randseqs),:]
        dat
    end
    ret, = calcGD(subset)
    ret
end

rarefy(data, N, reps::Int) = hcat([rare(data, N)[:GDval] for i in 1:reps]...)

# Read in the matrices for mammals and amphibians and concatenate them
am = ready_data("pairwise_equalarea_amphibians_cytb.csv")
ma = ready_data("pairwise_equalarea_mammals_cytb.csv")
dats = vcat(am, ma)

# Calculate GD values for those grid cells and species with less than 10 sequences
less_than_10s = by(dats, [:species, :cell]) do dat
    seqs = unique(vcat(dat[:seq1], dat[:seq2]))
    size(seqs, 1) < 10 ? dat : dat[1:0,:]
end

Supp_F1_data, = calc_cellvalues(less_than_10s)
writetable("Supp_Less10_numbers.csv", Supp_F1_data)

# Calculate GD values for those grid cells and species with less than 10 sequences
less_than_5s = by(dats, [:species, :cell]) do dat
    seqs = unique(vcat(dat[:seq1], dat[:seq2]))
    size(seqs, 1) < 5 ? dat : dat[1:0,:]
end

Supp_F4_data, = calc_cellvalues(less_than_5s)
writetable("Supp_Less5_numbers.csv", Supp_F4_data)

# Calculate GD values for those grid cell and species with more than 10 sequences
at_least_10s = by(dats, [:species, :cell]) do dat
    seqs = unique(vcat(dat[:seq1], dat[:seq2]))
    size(seqs, 1) >= 10 ? dat : dat[1:0,:]
end

Supp_F2_data, = calcGD(at_least_10s)
writetable("Supp_Atleast10_numbers.csv", Supp_F2_data)

# Calculate GD values for those grid cell and species with more than 5 sequences
at_least_5s = by(dats, [:species, :cell]) do dat
    seqs = unique(vcat(dat[:seq1], dat[:seq2]))
    size(seqs, 1) >= 5 ? dat : dat[1:0,:]
end

Supp_F5_data, = calcGD(at_least_5s)
writetable("Supp_AtLeast5_numbers.csv", Supp_F5_data)
# Calculate GD values from a rarefaction

for N in [5,10]
    # Create a matrix with a 100 realisations
    samples = rarefy(dats, N, 100)

    # Calculate mean and standard deviations from these samples
    rare_means = mean(samples, 2)
    rare_sds = std(samples, 2)
    Supp_f3_data = rarefy(dats, N)
    Supp_f3_data = DataFrame(cell = Supp_f3_data[:site], rare_mean = vec(rare_means), rare_sd = vec(rare_sds))

    writetable("Supp_Rarefy_numbers_$N.csv", Supp_f3_data)
    writecsv("All_rarefied_samples_$N.csv", samples)
end

Supp_F9_data, = calcGD_weighted(dats)
writetable("Supp_WeightedGD_numbers.csv", Supp_F9_data)

max_stars_repo_path: Miraldo et al. 2016/Calculate_Values_for_S5_to_S9.jl

max_stars_repo_name: mkborregaard/codes-for-published-papers

max_issues_repo_name: DalavanCloud/codes-for-published-papers

max_forks_repo_name: DalavanCloud/codes-for-published-papers


text:
# Recurrence parameters as defined by Marwan et al. (2007)

"""
    recurrencerate(x, theiler=0)

Calculate the recurrence rate (RR) of a recurrence matrix, ruling out
the points within the Theiler window.
"""
recurrencerate(x::AbstractMatrix) = countnz(x)/prod(size(x))

function recurrencerate(x::AbstractMatrix; theiler::Integer=0)
    theiler < 0 && error("Theiler window length must be greater than 0")
    diags_remove = -(theiler-1):(theiler-1)
    theiler_points = 0
    theiler_nz = 0
    for d in diags_remove
        theiler_points += length(diag(x,d))
        theiler_nz += countnz(diag(x,d))
    end
    (countnz(x)-theiler_nz)/(prod(size(x))-theiler_points)
end

function tau_recurrence(x::AbstractMatrix{Bool})
    n = minimum(size(x))
    [countnz(diag(x,d))/(n-d) for d in (1:n-1)]
end

# Based on diagonal lines

function diagonalhistogram(x::AbstractMatrix{Bool}; theiler::Integer=1, kwargs...)
    theiler < 0 && error("Theiler window length must be greater than 0")
    bins = [0]
    nbins = 1
    current_diag = 0
    m, n = size(x)
    # Iterate over diagonals - excluding LOI and Theiler window
    # If the matrix is symmetric, examine only the upper triangle
    diag_collection = collect(theiler:n-2)
    xsym = issym(x)
    !xsym && prepend!(diag_collection, collect(-(m-2):-max(theiler,1)))
    for d in diag_collection
        increment = (xsym && d > 0) ? 2 : 1
        previous_cell = false
        first_c = max(1, d+1)
        last_c = min(n, m+d)
        for c = first_c:last_c
            # Operate on existing diagonal line
            if previous_cell
                # Extend length with current cell
                (extend = x[c-d,c]) && (current_diag += 1)
                # If arrived to the end of a line
                # add the current length to the corresponding bin
                if (!extend || (c == last_c)) && (current_diag > 0)
                    # Append new positions to the bins if needed
                    if current_diag > nbins
                        append!(bins, zeros(current_diag - nbins))
                        nbins = current_diag
                    end
                    bins[current_diag] += increment
                    current_diag = 0
                end
            end
            previous_cell = x[c-d,c]
        end
    end
    # Add isolated points in first bin
    allpoints = (theiler == 0) ? countnz(x) : countnz(triu(x, theiler)) + countnz(tril(x,-theiler))
    [allpoints - collect(2:nbins+1)'*bins; bins]
end

"""
    determinism(x; lmin=2, theiler=1)

Calculate the determinism (DET) of a recurrence matrix, ruling out
the points within the Theiler window and diagonals shorter than a minimum value.
"""
function determinism(diag_hist::AbstractVector; lmin=2, kwargs...)
    if lmin < 2
        error("lmin must be 2 or higher")
    end
    nbins = length(diag_hist)
    diag_points = collect(1:nbins) .* diag_hist
    (lmin > nbins) ? 0.0 : sum(diag_points[lmin:nbins])/sum(diag_points)
end

function determinism(x::AbstractMatrix; kwargs...)
    determinism(diagonalhistogram(x; kwargs...); kwargs...)
end

"""
    avgdiag(x; lmin=2, theiler=1)

Calculate the average diagonal length (L) in a recurrence matrix, ruling out
the points within the Theiler window and diagonals shorter than a minimum value.
"""
function avgdiag(diag_hist::AbstractVector; lmin=2, kwargs...)
    if lmin < 2
        error("lmin must be 2 or higher")
    end
    nbins = length(diag_hist)
    diag_points = collect(1:nbins) .* diag_hist
    (lmin > nbins) ? 0.0 : sum(diag_points[lmin:nbins])/sum(diag_hist[lmin:nbins])
end

function avgdiag(x::AbstractMatrix; kwargs...)
    avgdiag(diagonalhistogram(x; kwargs...); kwargs...)
end

"""
    maxdiag(x; lmin=2, theiler=1)

Calculate the longest diagonal (Lmax) in a recurrence matrix, ruling out
the points within the Theiler window and diagonals shorter than a minimum value.
"""

maxdiag(diag_hist::AbstractVector) = length(diag_hist)
maxdiag(x::AbstractMatrix; kwargs...) = maxdiag(diagonalhistogram(x; kwargs...))

"""
    divergence(x; lmin=2, theiler=1)

Calculate the divergence of a recurrence matrix
(actually the inverse of `maxdiag`.
"""
divergence(x; kwargs...) = 1/maxdiag(x; kwargs...)

"""
    entropy(x; lmin=2, theiler=1)

Calculate the entropy of diagonal lengths (ENT) of a recurrence matrix, ruling out
the points within the Theiler window and diagonals shorter than a minimum value.
"""
function entropy(diag_hist::AbstractVector; lmin=2, kwargs...)
    if lmin < 2
        error("lmin must be 2 or higher")
    end
    nbins = length(diag_hist)
    if lmin <= nbins
        prob_bins = diag_hist[lmin:nbins] ./ sum(diag_hist[lmin:nbins])
        prob_bins = prob_bins[find(prob_bins)]
        -sum(prob_bins .* log(prob_bins))
    else
        0.0
    end
end

function entropy(x::AbstractMatrix; kwargs...)
    entropy(diagonalhistogram(x; kwargs...); kwargs...)
end

"""
    trend(x; theiler=1, border=10)

Calculate the trend of recurrences in recurrence matrix towards its edges, ruling out
the points within the Theiler window and in the outermost diagonals.
"""
function trend(npoints::AbstractVector; theiler=1, border=10, kwargs...)
    nmax = length(npoints)
    rrk = npoints./collect(nmax:-1:1)
    m = nmax-border
    w = collect(theiler:m)-m/2
    w'*(rrk[theiler:m]-mean(rrk[theiler:m])) / (w'*w)
end

function trend(x::AbstractMatrix; kwargs...)
    trend(tau_recurrence(x); kwargs...)
end

# Number of l-length sequences, based on diagonals
function countsequences(diag_hist::AbstractVector; lmin=2, kwargs...)
    overlap = (1:length(diag_hist))' - lmin + 1
    overlap[overlap .< 0] = 0
    overlap * diag_hist
end

function countsequences(x::AbstractMatrix; kwargs...)
    countsequences(diagonalhistogram(x; kwargs...); kwargs...)
end

# 2. Based on vertical lines

function verticalhistogram(x::AbstractMatrix{Bool})
    bins = [0]
    nbins = 1
    current_vert = 0
    m, n = size(x)
    # Iterate over columns
    for c = 1:n
        previous_cell = false
        for r = 1:m
            if previous_cell
                (extend = x[r,c]) && (current_vert += 1)
                if (!extend || (r == n)) && (current_vert > 0)
                    if current_vert > nbins
                        append!(bins, zeros(current_vert - nbins))
                        nbins = current_vert
                    end
                    bins[current_vert] += 1
                    current_vert = 0
                end
            end
            previous_cell = x[r,c]
        end
    end
    # Add isolated points in first bin
    [countnz(x) - collect(2:nbins+1)'*bins; bins]
end

"""
    laminarity(x; lmin=2)

Calculate the laminarity (LAM) of a recurrence matrix, ruling out
vertical lines shorter than a minimum value.
"""
laminarity(vert_hist::AbstractVector; kwargs...) = determinism(vert_hist; kwargs...)
laminarity(x::AbstractMatrix; kwargs...) = laminarity(verticalhistogram(x); kwargs...)

"""
    trappingtime(x; lmin=2)

Calculate the trapping time (TT) of a recurrence matrix, ruling out
vertical lines shorter than a minimum value.
"""
trappingtime(vert_hist::AbstractVector; kwargs...) = avgdiag(vert_hist; kwargs...)
trappingtime(x::AbstractMatrix; kwargs...) = trappingtime(verticalhistogram(x); kwargs...)

"""
    maxvert(x; lmin=2)

Calculate the longest vertical line (Vmax) of a recurrence matrix, ruling out
lines shorter than a minimum value.
"""
maxvert(vert_hist::AbstractVector) = length(vert_hist)
maxvert(x::AbstractMatrix; kwargs...) = maxvert(verticalhistogram(x; kwargs...))


text:
"""
    Parameters

Module containing sets of parameters.
"""
module Parameters

export AbstractParameterSet
export EarthParameterSet
abstract type AbstractParameterSet end

struct EarthParameterSet <: AbstractParameterSet end

include("Planet.jl")
include("Atmos.jl")
# include("Ocean.jl")
# include("Land.jl")
# include("Ice.jl")

end

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
using Pkg.Artifacts
using Tar, GZip, SHA

artifact_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")
version = "0.0.1"
uri = "https://github.com/barche/timestables-assets/archive/v$version.tar.gz"

tarball = download(uri)
hash = bytes2hex(open(sha256, tarball))

timeshash = create_artifact() do artifact_dir
    gzopen(tarball) do unzipped_tar
        Tar.extract(unzipped_tar, artifact_dir)
    end
end

bind_artifact!(artifact_toml, "timestables-assets", timeshash; download_info = [(uri, hash)], force=true)


text:
"""Generates *n* points at random on a *N*-dimensional ball."""
function on_ball(N::Int, n::Int)
    p = rand(Normal(), N, n)
    r = sqrt.(sum(abs2, p, dims=1))
    return p./r
end

"""Generates *n* points at random in a *N*-dimensional ball."""
function in_ball(N::Int, n::Int)
    X = reshape(rand(Normal(), n*N), N, n)
    Y = rand(Poisson(), n)
    r = sqrt.(Y' .+ sum(abs2, X, dims=1))
    return X./r
end

"""Generates *m* random *N*-dimensional translation vectors with *μ_bounds* bounds"""
function translations(N::Int, n::Int, μ_bounds::Tuple{Float64,Float64})
    o_min, o_max = μ_bounds
    Dτ = Uniform(o_min, o_max)
    o = on_ball(N, n)
    τ = rand(Dτ, n)
    return (o.*τ' .+ fill(1.0,N,n))./2.
end

"""Generates manifold basis and its complement"""
function lm_basis(N::Int, M::Int)
    @assert N > M
    X = rand(Normal(),N,N+1)
    Y = (X.-mean(X, dims=2))'./sqrt(N)
    F = svd(Y)
    F.V[:,1:M], F.V[:,(M+1):end]
end

"""Generates *n* points of a *M*-dimensional linear manifold cluster
situated in a hypercube."""
function generate_cluster(
                  n::Int,                   # number of points in generated manifold
                  μ::Vector{Float64},       # manifold translation vector
                  B::Matrix{Float64},       # manifold basis
                  B′::Matrix{Float64},      # orthogonal complement to manifold basis
                  DΦ::Vector{T},  # bound of manifold points
                  DE::Vector{T}   # bound of a point extent from a manifold
                ) where T<:Distribution
    N, M = size(B)
    @assert size(B′) == (N, N-M) "NULL space dimention must be $(N-M)"
    @assert length(DΦ) == M
    @assert length(DE) == N-M

    mfld = zeros(N,n)
    c = 1
    while c <= n
        ϕ = map(d->rand(d), DΦ)
        ɛ = map(d->rand(d), DE)
        x = μ + B*ϕ + B′*ɛ
        # check if manifold point outside of unit hypercube
        b = extrema(x)
        if b[1] < -1 || b[2] > 1.
            continue
        end
        mfld[:,c] = x
        c+=1
    end
    return mfld
end

"""Generates *n* points of a *M*-dimensional linear manifold cluster
situated in a hypercube."""
function generate_cluster(
                  n::Int,                   # number of points in generated manifold
                  μ::Vector{Float64},       # manifold translation vector
                  B::Matrix{Float64},       # manifold basis
                  B′::Matrix{Float64},      # orthogonal complement to manifold basis
                  Φ::Vector{Float64},       # bound of manifold points
                  E::Vector{Float64}        # bound of a point extent from a manifold
                )
    DΦ = map(ϕ->Uniform(-ϕ, ϕ), Φ)
    DE = map(ɛ->Uniform(-ɛ, ɛ), E)
    return generate_cluster(n, μ, B, B′, DΦ, DE)
end

"""Generates *n* points of a *M*-dimensional linear manifold cluster
situated in a hypercube."""
function generate_cluster(
                  n::Int,                   # number of points in generated manifold
                  μ::Vector{Float64},       # manifold translation vector
                  B::Matrix{Float64},       # manifold basis
                  B′::Matrix{Float64},      # orthogonal complement to manifold basis
                  Φ::Float64,               # bound of manifold points
                  E::Float64,               # bound of a point extent from a manifold
                )
    N, M = size(B)
    DΦ = fill(Uniform(-Φ, Φ), M)
    DE = fill(Uniform(-E, E), N-M)
    return generate_cluster(n, μ, B, B′, DΦ, DE)
end

"""Generates *n* points of a *M*-dimensional linear manifold cluster
situated in a hypercube using normal distribution."""
function generate_cluster(
                  ::Type{Normal},
                  n::Int,                   # number of points in generated manifold
                  μ::Vector{Float64},       # manifold translation vector
                  B::Matrix{Float64},       # full basis
                  M::Int,                   # cluster dimension
                  Φ::Float64,               # bound of manifold points
                  E::Float64,               # bound of a point extent from a manifold
                )
    N = size(B,1)
    DΦ = fill(Normal(0, Φ), M)
    DE = fill(Normal(0, E), N-M)
    return generate_cluster(n, μ, B[:,1:M], B[:, M+1:end], DΦ, DE)
end

"""Generates *n* points of a *M*-dimensional linear manifold cluster
situated in a hypercube using uniform distribution."""
function generate_cluster(
                  ::Type{Uniform},
                  n::Int,                   # number of points in generated manifold
                  μ::Vector{Float64},       # manifold translation vector
                  B::Matrix{Float64},       # full basis
                  M::Int,                   # cluster dimension
                  Φ::Float64,               # bound of manifold points
                  E::Float64,               # bound of a point extent from a manifold
                )
    N = size(B,1)
    DΦ = fill(Uniform(-Φ, Φ), M)
    DE = fill(Uniform(-E, E), N-M)
    return generate_cluster(n, μ, B[:,1:M], B[:, M+1:end], DΦ, DE)
end

""" Generates *m* linear manifold clusters
where *m* is size of parameter *M* """
function generate(n::Int,                     # number of points in generated manifold
                  N::Int,                     # space dimensionality
                  M::Vector{Int},             # generated manifolds dimensions
                  τ::Tuple{Float64,Float64},  # translation vector bounds
                  Φ::Vector{Vector{Float64}}, # bound of manifold points
                  E::Vector{Vector{Float64}}  # bound of a point extent from a manifold
                )
    m = length(M) # number of manifolds to generate

    @assert ( length(E) == m && length(Φ) == m ) "Bounds should correspond to number of generated manifolds"

    μ = translations(N, m, τ)
    bases = Matrix{Float64}[]
    manifolds = Matrix{Float64}[]
    for i in 1:m
        B, B′ = lm_basis(N, M[i])
        mfld = generate_cluster(n, μ[:,i], B, B′, Φ[i], E[i])
        push!(bases, B)
        push!(manifolds, mfld)
    end
    return manifolds, bases, μ
end

"""Generate manifold dimensions."""
function mdims(N, m, dim_cut::Int = 7)
    M = cumsum(ones(Int,m))
    if m >= N
        n = N >= dim_cut ? dim_cut : N
        M[n:end] .= n
    end
    return M
end

""" Generates *m* linear manifold clusters
where *m* is size of parameter *M* """
function generate(n::Int,                     # number of points in generated manifold
                  N::Int,                     # space dimensionality
                  M::Vector{Int},             # number of manifolds
                  τ::Tuple{Float64,Float64},   # translation vector bounds
                  κ::Float64
                )
    m = length(M)

    # generate bounds
    D = Uniform(0.2, 0.99); rand(D, 1000);
    ΔΦ = [sort(rand(D, 2))::Vector{Float64} for i in 1:m]
    ΔE = [sort(rand(D, 2)*diff(mean(ΔΦ))[1]/κ)::Vector{Float64} for i in 1:m]

    Φ = Vector{Float64}[]
    E = Vector{Float64}[]
    for i in M
        DΦ = Uniform(ΔΦ[i][1], ΔΦ[i][2]); rand(DΦ, 1000);
        DE = Uniform(ΔE[i][1], ΔE[i][2]); rand(DE, 1000);
        push!(Φ, rand(DΦ, M[i]))
        push!(E, rand(DE, N-M[i]))
    end
    return generate(n, N, M, τ, Φ, E)
end

""" Generates *m* linear manifold clusters"""
function generate(n::Int,                     # number of points in generated manifold
                  N::Int,                     # space dimensionality
                  m::Int,                     # number of manifolds
                  τ::Pair{Float64,Float64},   # translation vector bounds
                  κ::Float64;
                  dim_cut::Int = 7
                )
    # generate manifold dimensions
    M = mdims(N-1,m,dim_cut)
    return generate(n, N, M, (τ.first,τ.second), κ)
end

function noise(N::Int, n::Int)
    return rand(N, n)
end


text:
# module Lupin

function Base.getproperty(dict::T, key::Symbol) where {T <: Dict{String}}
    key in fieldnames(T) ? getfield(dict, key) : getindex(dict, String(key))
end

function dict_to_namedtuple(dict::Dict{String})
    ks = keys(dict)
    vals = map(collect(ks)) do k
        v = getindex(dict, k)
        val = if v isa Dict{String}
            dict_to_namedtuple(v)
        else
            v
        end
        Symbol(k) => val
    end
    (; vals...)
end

function NamedTuple(dict::Dict{String})
    dict_to_namedtuple(dict)
end

# module Lupin


text:
@testset "996.number-of-squareful-arrays.jl" begin
    @test num_squareful_perms([1,17,8]) == 2
    @test num_squareful_perms([2, 2, 2]) == 1
    @test num_squareful_perms([2, 2, 7]) == 3
end