text:
function _polygon_error_handler(request_body_dictionary::Dict{String, Any})::Tuple

    # initialize -
    error_response_dictionary = Dict{String,Any}()

    # what are my error keys?
    error_keys = keys(request_body_dictionary)
    for key ∈ error_keys
        error_response_dictionary[key] = request_body_dictionary[key]
    end

    # return -
    return (error_response_dictionary, nothing)
end

function _process_previous_close_polygon_call_response(body::String)

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error?
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        ticker=String[],
        volume=Float64[],
        volume_weighted_average_price=Float64[],
        open=Float64[],
        close=Float64[],
        high=Float64[],
        low=Float64[],
        timestamp=DateTime[],
        number_of_transactions=Int[]
    )

    # fill in the header dictionary -
    header_keys = [
        "ticker", "queryCount", "adjusted", "status", "request_id", "count"
    ]

    # check - do we have a count (if not resturn zero)
    get!(request_body_dictionary, "count", 0)

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # check - do we have a resultsCount field?
    results_count = get!(request_body_dictionary, "resultsCount", 0)
    if (results_count == 0) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    # populate the results DataFrame -
    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        # build a results tuple -
        result_tuple = (

            ticker = result_dictionary["T"],
            volume = result_dictionary["v"],
            volume_weighted_average_price = result_dictionary["vw"],
            open = result_dictionary["o"],
            close = result_dictionary["c"],
            high = result_dictionary["h"],
            low = result_dictionary["l"],
            timestamp = unix2datetime(result_dictionary["t"]*(1/1000)),
            number_of_transactions = result_dictionary["n"]
        )

        # push that tuple into the df -
        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_aggregates_polygon_call_response(body::String)

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error?
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        volume=Float64[],
        volume_weighted_average_price=Float64[],
        open=Float64[],
        close=Float64[],
        high=Float64[],
        low=Float64[],
        timestamp=DateTime[],
        number_of_transactions=Int[]
    )

    # fill in the header dictionary -
    header_keys = [
        "ticker", "queryCount", "adjusted", "status", "request_id", "count"
    ]

    # check - do we have a count (if not resturn zero)
    get!(request_body_dictionary, "count", 0)
    get!(request_body_dictionary, "ticker", "N/A")

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # check - do we have a resultsCount field?
    results_count = get!(request_body_dictionary, "resultsCount", 0)
    if (results_count == 0) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    # populate the results DataFrame -
    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        # set some defaults in case missing fields -
        get!(result_dictionary, "vw", 0.0)
        get!(result_dictionary, "n", 0)

        # build a results tuple -
        result_tuple = (

            volume = result_dictionary["v"],
            volume_weighted_average_price = result_dictionary["vw"],
            open = result_dictionary["o"],
            close = result_dictionary["c"],
            high = result_dictionary["h"],
            low = result_dictionary["l"],
            timestamp = unix2datetime(result_dictionary["t"]*(1/1000)),
            number_of_transactions = result_dictionary["n"]
        )

        # push that tuple into the df -
        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_daily_open_close_call_response(body::String)

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error?
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR" || status_flag == "NOT_FOUND")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        ticker=String[],
        volume=Float64[],
        open=Float64[],
        close=Float64[],
        high=Float64[],
        low=Float64[],
        from=Date[],
        afterHours=Float64[],
        preMarket=Float64[]
    )

    header_keys = [
        "status"
    ]
    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # build a results tuple -
    result_tuple = (

        volume = request_body_dictionary["volume"],
        open = request_body_dictionary["open"],
        close = request_body_dictionary["close"],
        high = request_body_dictionary["high"],
        low = request_body_dictionary["low"],
        from = Date(request_body_dictionary["from"]),
        ticker = request_body_dictionary["symbol"],
        afterHours = request_body_dictionary["afterHours"],
        preMarket = request_body_dictionary["preMarket"]
    )

    # push that tuple into the df -
    push!(df, result_tuple)

    # return -
    return (header_dictionary, df)
end

function _process_options_reference_call_response(body::String)

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error?
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        cfi=String[],
        contract_type=String[],
        exercise_style=String[],
        expiration_date=Date[],
        primary_exchange=String[],
        shares_per_contract=Int64[],
        strike_price=Float64[],
        ticker = String[],
        underlying_ticker = String[]
    )

    # fill in the header dictionary -
    header_keys = [
        "status", "request_id", "count", "next_url"
    ]
    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # populate the results DataFrame -
    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        # build a results tuple -
        result_tuple = (

            cfi = result_dictionary["cfi"],
            contract_type = result_dictionary["contract_type"],
            exercise_style = result_dictionary["exercise_style"],
            expiration_date = Date(result_dictionary["expiration_date"]),
            primary_exchange = result_dictionary["primary_exchange"],
            shares_per_contract = result_dictionary["shares_per_contract"],
            strike_price = result_dictionary["strike_price"],
            ticker = result_dictionary["ticker"],
            underlying_ticker = result_dictionary["underlying_ticker"]
        )

        # push that tuple into the df -
        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_ticker_news_call_response(body::String)

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error?
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame();

    # fill in the header dictionary -
    header_keys = [
        "status", "request_id", "count", "next_url"
    ]
    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # return -
    return (header_dictionary, df)
end

function _process_ticker_details_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error? can be due to stick or date
    status_flag = request_body_dictionary["status"]
    if (status_flag == "NOT_FOUND" || status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

            ticker = String[],
            name = String[],
            market = String[],
            locale = String[],
            primary_exchange = String[],
            type = String[],
            active = Bool[],
            currency_name = String[],
            cik = String[],
            composite_figi = String[],
            share_class_figi = String[],
            market_cap = Float64[],
            phone_number = String[],
            address1 = String[],
            city = String[],
            state = String[],
            postal_code = String[],
            description = String[],
            sic_code = String[],
            sic_description = String[],
            ticker_root = String[],
            homepage_url = String[],
            total_employees = Float64[],
            list_date = String[],
            logo_url = String[],
            icon_url = String[],
            share_class_shares_outstanding = Float64[],
            weighted_shares_outstanding = Float64[],
            delisted_utc = String[]
        )

    # fill in the header dictionary -
    header_keys = [
                "status", "request_id"
        ];

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # if no results we return nothing
    if (request_body_dictionary["results"] == Any[]) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    results_array = request_body_dictionary["results"]

    #get if no values
    get!(results_array, "primary_exchange", "N/A")
    get!(results_array, "type", "N/A")
    get!(results_array, "cik", "N/A")
    get!(results_array, "composite_figi", "N/A")
    get!(results_array, "share_class_figi", "N/A")
    get!(results_array, "market_cap", NaN)
    get!(results_array, "phone_number", "N/A")

    #get if no values - address dictionary
    get!(results_array, "address", Dict("address1" => "N/A", "city" => "N/A", "state" => "N/A", "postal_code" => "N/A"))
    get!(results_array["address"], "address1", "N/A")
    get!(results_array["address"], "city", "N/A")
    get!(results_array["address"], "state", "N/A")
    get!(results_array["address"], "postal_code", "N/A")

    get!(results_array, "description", "N/A")
    get!(results_array, "sic_code", "N/A")
    get!(results_array, "sic_description", "N/A")
    get!(results_array, "ticker_root", "N/A")
    get!(results_array, "homepage_url", "N/A")
    get!(results_array, "total_employees", NaN)
    get!(results_array, "list_date", "N/A")

    #get if no values - branding dictionary
    get!(results_array, "branding", Dict("logo_url" => "N/A",  "icon_url" => "N/A"))
    get!(results_array["branding"], "logo_url", "N/A")
    get!(results_array["branding"], "icon_url", "N/A")

    get!(results_array, "share_class_shares_outstanding", NaN)
    get!(results_array, "share_class_shares_outstanding", NaN)
    get!(results_array, "delisted_utc", "N/A")

    result_tuple = (

                ticker = results_array["ticker"],
                name = results_array["name"],
                market = results_array["market"],
                locale = results_array["locale"],
                primary_exchange = results_array["primary_exchange"],
                type = results_array["type"],
                active = results_array["active"],
                currency_name = results_array["currency_name"],
                cik = results_array["cik"],
                composite_figi = results_array["composite_figi"],
                share_class_figi = results_array["share_class_figi"],
                market_cap = results_array["market_cap"],
                phone_number = results_array["phone_number"],
                address1 = results_array["address"]["address1"],
                city = results_array["address"]["city"],
                state = results_array["address"]["state"],
                postal_code = results_array["address"]["postal_code"],
                description = results_array["description"],
                sic_code = results_array["sic_code"],
                sic_description = results_array["sic_description"],
                ticker_root = results_array["ticker_root"],
                homepage_url = results_array["homepage_url"],
                total_employees = results_array["total_employees"],
                list_date = results_array["list_date"],
                logo_url = results_array["branding"]["logo_url"],
                icon_url = results_array["branding"]["icon_url"],
                share_class_shares_outstanding = results_array["share_class_shares_outstanding"],
                weighted_shares_outstanding = results_array["weighted_shares_outstanding"],
                delisted_utc = results_array["delisted_utc"]
            )
    push!(df, result_tuple)

    # return -
    return (header_dictionary, df)
end

function _process_market_holidays_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

# initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        exchange = String[],
        name = String[],
        date = Date[],
        status = String[],
        open = String[],
        close = String[]
    );

    # populate the results DataFrame -
    results_array = request_body_dictionary
    for result_dictionary ∈ results_array

        get!(result_dictionary, "open", "N/A")
        get!(result_dictionary, "close", "N/A")

        # build a results tuple -
        result_tuple = (

            exchange = result_dictionary["exchange"],
            name = result_dictionary["name"],
            date = Date(result_dictionary["date"]),
            status = result_dictionary["status"],
            open = result_dictionary["open"],
            close = result_dictionary["close"]
        )

        # push that tuple into the df -
        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_exchanges_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error?
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        id = Int[],
        type = String[],
        asset_class = String[],
        locale = String[],
        name = String[],
        acronym = String[],
        mic = String[],
        operating_mic = String[],
        participant_id = String[],
        url = String[]

    );

    # fill in the header dictionary -
    header_keys = [
            "status", "request_id", "count"
    ];

    # check - do we have a count (if not resturn zero)
    results_count = get!(request_body_dictionary, "count", 0)

    # if no result we return nothing
    if (results_count == 0) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # populate the results DataFrame -
    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array
        # set some defaults in case missing fields -
        get!(result_dictionary, "acronym", "N/A")
        get!(result_dictionary, "mic", "N/A")
        get!(result_dictionary, "operating_mic", "N/A")
        get!(result_dictionary, "participant_id", "N/A")
        get!(result_dictionary, "url", "N/A")

        result_tuple = (

                    id = result_dictionary["id"],
                    type = result_dictionary["type"],
                    asset_class = result_dictionary["asset_class"],
                    locale = result_dictionary["locale"],
                    name = result_dictionary["name"],
                    acronym = result_dictionary["acronym"],
                    mic = result_dictionary["mic"],
                    operating_mic = result_dictionary["operating_mic"],
                    participant_id = result_dictionary["participant_id"],
                    url = result_dictionary["url"]
                )

        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_stock_splits_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        execution_date = Date[],
        split_from = Int[],
        split_to = Int[],
        ticker = String[]

    );

    # fill in the header dictionary -
    header_keys = [
            "status", "request_id"
    ];

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # if no result we return nothing
    if (request_body_dictionary["results"] == Any[]) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    # populate the results DataFrame -
    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        result_tuple = (

                    execution_date = Date(result_dictionary["execution_date"]),
                    split_from = result_dictionary["split_from"],
                    split_to = result_dictionary["split_to"],
                    ticker = result_dictionary["ticker"]
                )

        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_market_status_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    #initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

    serverTime = String[],
    market = String[],
    earlyHours = Bool[],
    afterHours = Bool[],
    nyse  = String[],
    nasdaq  = String[],
    otc = String[],
    fx  = String[],
    crypto  = String[]
    );

    # populate the results DataFrame -
    results_array = request_body_dictionary

    # build a results tuple -
    result_tuple = (

        serverTime = results_array["serverTime"],
        market = results_array["market"],
        earlyHours = results_array["earlyHours"],
        afterHours = results_array["afterHours"],
        nyse = results_array["exchanges"]["nyse"],
        nasdaq = results_array["exchanges"]["nasdaq"],
        otc = results_array["exchanges"]["otc"],
        fx = results_array["currencies"]["fx"],
        crypto = results_array["currencies"]["crypto"]
        )

    # push that tuple into the df -
    push!(df, result_tuple)

    # return -
    return (header_dictionary, df)
end

function _process_dividends_call_response(body::String)

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)

    # before we do anything - check: do we have an error? can be due to stick or date
    status_flag = request_body_dictionary["status"]
    if (status_flag == "NOT_FOUND" || status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    #initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

            cash_amount = Float64[],
            declaration_date = Date[],
            dividend_type = String[],
            ex_dividend_date = Date[],
            frequency = Int[],
            pay_date = Date[],
            record_date = Date[],
            ticker =String[]
        )

    # fill in the header dictionary -
    header_keys = [
                "status", "request_id","next_url"
        ];

    #fill in next_url if no value
    get!(request_body_dictionary,"next_url","N/A")

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # if no results we return nothing
    if (request_body_dictionary["results"] == Any[]) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        result_tuple = (

                    cash_amount = result_dictionary["cash_amount"],
                    declaration_date = Date(result_dictionary["declaration_date"]),
                    dividend_type = result_dictionary["dividend_type"],
                    ex_dividend_date = Date(result_dictionary["ex_dividend_date"]),
                    frequency = result_dictionary["frequency"],
                    pay_date = Date(result_dictionary["pay_date"]),
                    record_date = Date(result_dictionary["record_date"]),
                    ticker = result_dictionary["ticker"],

                )

        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_tickers_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)
    # before we do anything - check: do we have an error? can be due to stick or date
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        ticker = String[],
        name = String[],
        market = String[],
        locale = String[],
        active = Bool[],
        primary_exchange = String[],
        type = String[],
        currency_name = String[],
        cik = String[],
        composite_figi = String[],
        share_class_figi = String[],
        last_updated_utc = String[],
        delisted_utc = String[]

        )

    # fill in the header dictionary -
    header_keys = [
                "status", "request_id", "count", "next_url"
        ];

    #fill in next_url if no value
    get!(request_body_dictionary,"next_url","N/A")
    get!(request_body_dictionary,"count",0)

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # if no results we return nothing
    if (request_body_dictionary["results"] == Any[]) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        #get if no values
        get!(result_dictionary, "primary_exchange", "N/A")
        get!(result_dictionary, "type", "N/A")
        get!(result_dictionary, "currency_name", "N/A")
        get!(result_dictionary, "cik", "N/A")
        get!(result_dictionary, "composite_figi", "N/A")
        get!(result_dictionary, "share_class_figi", "N/A")
        get!(result_dictionary, "last_updated_utc", "N/A")
        get!(result_dictionary, "delisted_utc", "N/A")

        result_tuple = (

                    ticker = result_dictionary["ticker"],
                    name = result_dictionary["name"],
                    market = result_dictionary["market"],
                    locale = result_dictionary["locale"],
                    active = result_dictionary["active"],
                    primary_exchange = result_dictionary["primary_exchange"],
                    type = result_dictionary["type"],
                    currency_name = result_dictionary["currency_name"],
                    cik = result_dictionary["cik"],
                    composite_figi = result_dictionary["composite_figi"],
                    share_class_figi = result_dictionary["share_class_figi"],
                    last_updated_utc = result_dictionary["last_updated_utc"],
                    delisted_utc = result_dictionary["delisted_utc"]
                )

        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_conditions_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)
    # before we do anything - check: do we have an error? can be due to stick or date
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        name = String[],#
        asset_class = String[],#
        data_types =  Array{Array{String,1},1}(),#
        id = Int[],#
        abbreviation = String[],
        description = String[],
        exchange = String[],
        legacy = String[],
        sip_mapping_CTA = String[],
        sip_mapping_OPRA = String[],
        sip_mapping_UTP = String[],
        type = String[],#
        consolidated_updates_high_low = String[],
        consolidated_updates_open_close = String[],
        consolidated_updates_volume = String[],
        market_center_updates_high_low = String[],
        market_center_updates_open_close = String[],
        market_center_updates_volume = String[]

        )

    # fill in the header dictionary -
    header_keys = [
                "status", "request_id", "count", "next_url"
        ];

    #fill in next_url if no value
    get!(request_body_dictionary,"next_url","N/A")
    get!(request_body_dictionary,"count",0)

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # if no results we return nothing
    if (request_body_dictionary["results"] == Any[]) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        #get if no values
        get!(result_dictionary, "abbreviation", "N/A")
        get!(result_dictionary, "description", "N/A")
        get!(result_dictionary, "exchange", "N/A")
        get!(result_dictionary, "legacy", "N/A")
        get!(result_dictionary["sip_mapping"], "CTA", "N/A")
        get!(result_dictionary["sip_mapping"], "OPRA", "N/A")
        get!(result_dictionary["sip_mapping"], "UTP", "N/A")
        get!(result_dictionary, "update_rules", Dict("consolidated"  => Dict("updates_high_low"   => "N/A",
                                                                             "updates_open_close" => "N/A",
                                                                             "updates_volume"     => "N/A"),
                                                     "market_center" => Dict("updates_high_low"   => "N/A",
                                                                             "updates_open_close" => "N/A",
                                                                             "updates_volume"     => "N/A")))
        result_tuple = (

                    name = result_dictionary["name"],
                    asset_class = result_dictionary["asset_class"],
                    data_types = result_dictionary["data_types"],
                    id = result_dictionary["id"],
                    abbreviation = result_dictionary["abbreviation"],
                    description = result_dictionary["description"],
                    exchange = string(result_dictionary["exchange"]),
                    legacy = string(result_dictionary["legacy"]),
                    sip_mapping_CTA = result_dictionary["sip_mapping"]["CTA"],
                    sip_mapping_OPRA = result_dictionary["sip_mapping"]["OPRA"],
                    sip_mapping_UTP = result_dictionary["sip_mapping"]["UTP"],
                    type = result_dictionary["type"],
                    consolidated_updates_high_low = string(result_dictionary["update_rules"]["consolidated"]["updates_high_low"]),
                    consolidated_updates_open_close = string(result_dictionary["update_rules"]["consolidated"]["updates_open_close"]),
                    consolidated_updates_volume = string(result_dictionary["update_rules"]["consolidated"]["updates_volume"]),
                    market_center_updates_high_low = string(result_dictionary["update_rules"]["market_center"]["updates_high_low"]),
                    market_center_updates_open_close = string(result_dictionary["update_rules"]["market_center"]["updates_open_close"]),
                    market_center_updates_volume = string(result_dictionary["update_rules"]["market_center"]["updates_volume"])

                )

        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end

function _process_stock_financials_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)
    # before we do anything - check: do we have an error? can be due to stick or date
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String, Any}()
    df = DataFrame(

        financials = String[],
        content = String[],
        label = String[],
        order = Int[],
        unit = String[],
        value = Float64[],
        formula = String[],
        xpath = String[]
        )

    # fill in the header dictionary -
    header_keys = [
                "status", "request_id", "count", "next_url"
        ];

    #addressing missing value in header dictionary -
    get!(request_body_dictionary,"next_url","N/A")
    get!(request_body_dictionary,"count",0)

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # if no results we return nothing
    if (request_body_dictionary["results"] == Any[]) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    #pull out financials dictionary
    results_array = request_body_dictionary["results"]
    financials_dictionary = Dict{String, Any}()

    for result_dictionary ∈ results_array

        #addressing missing values
        get!(result_dictionary, "start_date", "N/A")
        get!(result_dictionary, "end_date", "N/A")

        #update header dictionary
        header_dictionary["company_name"] = result_dictionary["company_name"]
        header_dictionary["cik"] = result_dictionary["cik"]
        header_dictionary["fiscal_period"] = result_dictionary["fiscal_period"]
        header_dictionary["fiscal_year"] = result_dictionary["fiscal_year"]
        header_dictionary["source_filing_file_url"] = result_dictionary["source_filing_file_url"]
        header_dictionary["source_filing_url"] = result_dictionary["source_filing_url"]
        header_dictionary["start_date"] = result_dictionary["start_date"]
        header_dictionary["end_date"] = result_dictionary["end_date"]

        #pull out financials dictionary
        financials_dictionary =  result_dictionary["financials"]

    end

    #filling dataframe with desired information
    for i in keys(financials_dictionary)
        for j in keys(financials_dictionary[i])

            #addressing missing values
            get!(financials_dictionary[i][j], "formula", "N/A")
            get!(financials_dictionary[i][j], "xpath", "N/A")
            result_tuple = (

                financials = i,
                content = j,
                label = financials_dictionary[i][j]["label"],
                order = financials_dictionary[i][j]["order"],
                unit = financials_dictionary[i][j]["unit"],
                value = financials_dictionary[i][j]["value"],
                formula = financials_dictionary[i][j]["formula"],
                xpath = financials_dictionary[i][j]["xpath"]
            )
            push!(df, result_tuple)
        end
    end

    # return -
    return (header_dictionary, df)
end

function _process_ticker_types_call_response(body::String) #ycpan

    # convert to JSON -
    request_body_dictionary = JSON.parse(body)
    # before we do anything - check: do we have an error? can be due to stick or date
    status_flag = request_body_dictionary["status"]
    if (status_flag == "ERROR")
        return _polygon_error_handler(request_body_dictionary)
    end

    # initialize -
    header_dictionary = Dict{String,Any}()
    df = DataFrame(

        code = String[],
        description = String[],
        asset_class = String[],
        locale = String[]
        )

    # fill in the header dictionary -
    header_keys = [
                "status", "request_id", "count"
        ];

    #fill in next_url if no value
    get!(request_body_dictionary,"count",0)

    for key ∈ header_keys
        header_dictionary[key] = request_body_dictionary[key]
    end

    # if no results we return nothing
    if (request_body_dictionary["results"] == Any[]) # we have no results ...
        # return the header and nothing -
        return (header_dictionary, nothing)
    end

    results_array = request_body_dictionary["results"]
    for result_dictionary ∈ results_array

        result_tuple = (

                    code = result_dictionary["code"],
                    description = result_dictionary["description"],
                    asset_class = result_dictionary["asset_class"],
                    locale = result_dictionary["locale"]
                )

        push!(df, result_tuple)
    end

    # return -
    return (header_dictionary, df)
end


text:
# %% permutation entropy
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

using Random, Combinatorics, Statistics
Random.seed!(356)
o = 3
N = 8
x = rand(N)

fig = figure(figsize = (figx, figx/3.5))
ax1 = subplot2grid((1, 3), (0, 0))
ax1.plot(1:N, x, marker = "o", color = "C0", mfc = "C2", ms = 15, zorder = 99)
ax1.set_title("timeseries", size = 28)
ax1.spines["bottom"].set_position("center")
ax1.spines["right"].set_color("none")
ax1.spines["top"].set_color("none")
ax1.set_xticks([])
ax1.set_yticks([])
ax1.set_xlabel("\$t\$")
ax1.set_ylabel("\$x\$", rotation = 0)
ax1.xaxis.set_label_coords(1.0, 0.5)
ax1.yaxis.set_label_coords(-0.05, 0.9)
ax1.set_ylim(-0.1, 1.1)
# ax1.axis("off")

# add pattern indications
for (s, n) in zip(("2", "3", "3"), (2, 4, 6))
    xs = [n, n, n+o-1, n+o-1]
    ym = minimum(x[xs[1]:xs[end]]) - 0.1
    yd = ym - 0.02
    ys = [ym, yd, yd, ym] #.+ 0.05
    ax1.plot(xs, ys, color = "k", lw = 1)
    ax1.text(mean(xs), minimum(ys)-0.02, "#=$s", ha="center", va = "top")
end

# ax1.annotate("4", (2, -0.2), xytext = (4, -0.2), ha="center")

ax2 = subplot2grid((1, 3), (0, 1), colspan = 2)
p = permutations([0, 0.5, 1], o) |> collect
counts = [0, 1, 3, 1, 0, 1]
for (i, a) in enumerate(p)
    ax2.plot((1:o) .+ i*o .+ 1, a, marker = "o",
    color = "C0", mfc = "C2", ms = 10)
    ax2.text(o÷2 + i*o + 1, 1.2, "$i")
    ax2.text(o÷2 + i*o + 1, -0.5, "$(counts[i])")
end
ax2.text(o+1, 1.2, "#", ha = "right")
ax2.text(o+1, 0.5, "pattern", ha = "right")
ax2.text(o+1, -0.5, "count", ha = "right")
ax2.set_ylim(-1, 1.5)
ax2.set_xlim(0, ax2.get_xlim()[2])
ax2.set_title("relative amplitude permutations, \$d=$o\$", size = 28)
ax2.axis("off")

fig.tight_layout()
fig.subplots_adjust(top = 0.9, bottom = 0.02, right = 0.95, left = 0.05, wspace=0.1, hspace = 0.1)

wsave(plotsdir("6", "permentropy"), fig)

ext: jl

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# Advent of Code 2020 - Day 22

test_str = "Player 1:
9
2
6
3
1

Player 2:
5
8
4
7
10"

test_list = split(test_str, '\n')

list = readlines("day22_input.txt")

function play_combat(list)
    num_cards = div(length(list)+1, 2) - 2
    decks = fill(0, 2, 2*num_cards)
    decks[1, 1:num_cards] = parse.(Int, list[2:num_cards+1])
    decks[2, 1:num_cards] = parse.(Int, list[num_cards+4:end])
    N = num_cards
    while 0 < N < 2*num_cards
        cards = decks[:, 1]
        decks[:, 1:2*num_cards-1] = decks[:, 2:end]
        winner = argmax(cards)
        slot = winner == 1 ? N : 2*num_cards-N
        N += winner == 1 ? 1 : -1
        decks[winner, slot:slot+1] = [cards[winner] cards[3-winner]]
    end

    winner = N == 0 ? 2 : 1
    points = sum([decks[winner, j]*(2*num_cards-j+1) for j in 1:2*num_cards])

    return points
end

@show play_combat(test_list) == 306
@show play_combat(list) == 32495

function play_recursive_combat(list)
    num_cards = div(length(list)+1, 2) - 2
    decks = fill(0, 2, 2*num_cards)
    decks[1, 1:num_cards] = parse.(Int, list[2:num_cards+1])
    decks[2, 1:num_cards] = parse.(Int, list[num_cards+4:end])

    function playgame(decks)
        num_cards = div(size(decks)[2], 2)
        history1 = Set()
        history2 = Set()
        gamerepeated = false
        while minimum(sum(decks, dims=2)) > 0 && !gamerepeated
            if  decks[1,:] ∈ history1 || decks[2,:] ∈ history2
                gamerepeated = true
                winner = 1
            else
                push!(history1, decks[1,:])
                push!(history2, decks[2,:])
                cards = decks[:, 1]
                decks[:, 1:2*num_cards-1] = decks[:, 2:end]
                if sum(sum(decks.>0, dims=2) .≥ cards) == 2
                    new_decks = copy(decks)
                    new_decks[1,cards[1]+1:end] .= 0
                    new_decks[2,cards[2]+1:end] .= 0
                    winner = playgame(new_decks)
                else
                    winner = argmax(cards)
                end
                slot = findfirst(==(0), decks[winner,:])
                decks[winner, slot:slot+1] = [cards[winner] cards[3-winner]]
            end
        end
        return winner
    end

    winner = playgame(decks)

    points = sum([decks[winner, j]*(2*num_cards-j+1) for j in 1:2*num_cards])
    return points
end

@show play_recursive_combat(test_list) == 291
@show play_recursive_combat(list) == 32665

nothing


text:
using RCall
using DataFrames
using Turing

bayesplot =  rimport("bayesplot")

# https://github.com/JuliaInterop/RCall.jl/issues/333
function mcmc_recover_hist(posterior, true_parameters; kwargs...)
    #TableTraits.isiterabletable(posterior) || error("Argument `posterior` is not an iterable table.")
    df = DataFrame(posterior)
    return mcmc_recover_hist(df, true_parameters; kwargs...)
end
function mcmc_recover_hist(df::DataFrame, true_parameters; kwargs...)
    cols = Symbol.(names(df))
    true_params = Float64[]
    for c in cols
        val = get(true_parameters, c) do
            error("Did not find key $c in `true_parameters`.")
        end
        push!(true_params, val)
    end
    return bayesplot.mcmc_recover_hist(df, true_params; kwargs...)
end
function mcmc_recover_hist(posterior::MCMCChains.Chains, true_parameters; kwargs...)
    #TableTraits.isiterabletable(posterior) || error("Argument `posterior` is not an iterable table.")
    df = DataFrame(MCMCChains.get_sections(posterior, :parameters))
    df = df[:, Not([:iteration, :chain])]
    return mcmc_recover_hist(df, true_parameters; kwargs...)
end

function _ppc_conv(y, yrep::Vector{<:NamedTuple})
    @assert length(first(yrep)) == 1
    @assert length(y) == length(yrep[1][1])
    ymat = transpose(reduce(hcat, map(first, yrep)))
    (y, ymat)
end

function ppc_dens_overlay(y, yrep; kwargs...)
    bayesplot.ppc_dens_overlay(_ppc_conv(y, yrep)...; kwargs...)
end

function ppc_rootogram(y, yrep; kwargs...)
    bayesplot.ppc_rootogram(_ppc_conv(y, yrep)...; kwargs...)
end

function ppc_stat(y, yrep; kwargs...)
    bayesplot.ppc_stat(_ppc_conv(y, yrep)...; kwargs...)
end

function ppc_intervals(y, yrep; kwargs...)
    bayesplot.ppc_intervals(_ppc_conv(y, yrep)...; kwargs...)
end

function ppc_stat_grouped(y, yrep; kwargs...)
    bayesplot.ppc_stat_grouped(_ppc_conv(y, yrep)...; kwargs...)
end

function mcmc_hist(x; kwargs...)
    bayesplot.mcmc_hist(x; kwargs...)
end

function mcmc_trace(x; kwargs...)
    bayesplot.mcmc_trace(x; kwargs...)
end

function mcmc_scatter(x; kwargs...)
    bayesplot.mcmc_scatter(x; kwargs...)
end

function mcmc_parcoord(x; kwargs...)
    bayesplot.mcmc_parcoord(x; kwargs...)
end

function mcmc_areas(x; kwargs...)
    bayesplot.mcmc_areas(x; kwargs...)
end

function bayesplot_grid(plots...; kwargs...)
    bayesplot.bayesplot_grid(plots...; kwargs...)
end


text:
using Random
using Distributions
using MultiAgentPathFinding
using MultiAgentAllocationTransit
using BenchmarkTools
using Statistics
using JSON
using Logging
using JLD2
using NearestNeighbors
using StaticArrays
global_logger(SimpleLogger(stderr, Logging.Error))

rng = MersenneTwister(2345)

# Script arguments relating to transit files etc.
const city_params_file = ARGS[1]
const stop_coords_file = ARGS[2]
const trips_file = ARGS[3]
const drone_params_file = ARGS[4]
const bb_params_file = ARGS[5]
const out_file_pref = ARGS[6]
const city_travel_time_estimates = ARGS[7]
const N_DEPOTS = parse(Int64, ARGS[8])
const N_AGENTS = parse(Int64, ARGS[9])
const N_TRIALS = parse(Int64, ARGS[10])

const TRANSIT_CAP_RANGE = (3, 5)
const ECBS_WEIGHT = 1.1

## NOTE: EXAMPLE CALL
# julia scripts/benchmark_replanning.jl data/sfmta/sf_params.toml data/sfmta/stop_to_coords.json data/sfmta/trips.json
# data/drone_params.toml data/sfmta/sf_bb_params.toml data/sf_replan_20trials data/sfmta/sf_halton_tt_estimates.jld2 5 10 20

function main()

    # Load the various parameters
    city_params = parse_city_params(city_params_file)
    bb_params = parse_city_params(bb_params_file)
    drone_params = parse_drone_params(drone_params_file)

    lat_dist = Uniform(bb_params.lat_start, bb_params.lat_end)
    lon_dist = Uniform(bb_params.lon_start, bb_params.lon_end)

    # Transit Graph Preprocessing
    tg = load_transit_graph_latlong(stop_coords_file, trips_file, TRANSIT_CAP_RANGE, rng)
    tg, stop_idx_to_trips = transit_graph_preprocessing(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean, drone_params)

    # Load Halton stuff
    @load city_travel_time_estimates halton_nn_tree city_halton_points travel_time_estimates

    replan_results = Dict("cap_range"=>TRANSIT_CAP_RANGE, "weight"=>ECBS_WEIGHT,
                        "depots"=>N_DEPOTS, "agents"=>N_AGENTS)

    rp_one_times = Float64[]
    rp_one_msps = Float64[]
    rp_all_times = Float64[]
    rp_all_msps = Float64[]

    # No need to ignore first trial now!
    TRIAL = 1
    while TRIAL <= N_TRIALS+1

        N_SITES = 2*N_AGENTS

        println("Trial $(TRIAL)")

        # Generate depots and sites
        depots = [LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist))) for i = 1:N_DEPOTS]
        sites = [LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist))) for i = 1:N_SITES]
        depot_sites = vcat(depots, sites)

        # Off transit preprocessing
        otg = OffTransitGraph(depots = depots, sites = sites)
        aug_trips_fws_dists = augmented_trip_meta_graph_fws_dists(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                                                              length(depots), length(sites),
                                                              vcat(depots, sites),
                                                              drone_params)
        state_graph, depot_sites_to_vtx, trip_to_vtx_range = setup_state_graph(tg, otg)

        # Set up the shadow env to do task allocation with
        env = MAPFTransitEnv(off_transit_graph = otg, transit_graph = tg, state_graph = state_graph,
                             agent_states = AgentState[], depot_sites_to_vtx = depot_sites_to_vtx, trip_to_vtx_range = trip_to_vtx_range,
                             stop_idx_to_trips = stop_idx_to_trips, aug_trips_fws_dists = aug_trips_fws_dists,
                             drone_params = drone_params, dist_fn = MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                             curr_site_points = [], threshold_global_conflicts = 10)

        # run the task allocation, obtain the agent tasks and true number of agents
        cost_fn(i, j) = allocation_cost_wrapper_estimate(env, ECBS_WEIGHT, N_DEPOTS, N_SITES,
                                                        halton_nn_tree, travel_time_estimates, i, j)
        agent_tours = task_allocation(N_DEPOTS, N_SITES, N_AGENTS,
                                      depot_sites, cost_fn)

        println("Task Allocation Done!")

        agent_tasks = get_agent_task_set(agent_tours, N_DEPOTS, N_SITES)
        true_n_agents = length(agent_tasks)

        # Now set n_agent dependent env params
        env.agent_states = [AgentState(task=agt_task) for agt_task in agent_tasks]
        env.curr_site_points = zeros(Int64, true_n_agents)
        env.plan_ref_times = zeros(true_n_agents)

        # Setup ECBS Solver; Run search
        initial_states = Vector{MAPFTransitVertexState}(undef, true_n_agents)
        for i = 1:true_n_agents
            initial_states[i] = env.state_graph.vertices[env.depot_sites_to_vtx[string("d-", agent_tasks[i].origin)]]
        end

        solver = ECBSSolver{MAPFTransitVertexState,MAPFTransitAction,Float64,Makespan,MAPFTransitConflict,MAPFTransitConstraints,MAPFTransitEnv}(env = env, weight = ECBS_WEIGHT)

        solution = search!(solver, initial_states)

        println("Initial solution obtained!")

        # Compute env diagnostics and throw away if invalid
        set_solution_diagnostics!(env, solution)
        n_valid_firstpaths = length(env.valid_path_dists)

        if n_valid_firstpaths > 0.67 * true_n_agents

            # Copy env and solution for replanning
            env_copy = deepcopy(env)
            solution_copy = deepcopy(solution)

            # Actually run the replanning
            did_replan_indiv, el_time_indiv = replan_individual!(env, solution, N_DEPOTS, N_SITES, agent_tours, ECBS_WEIGHT)
            did_replan_collec, el_time_collec, new_soln_collec = replan_collective!(env_copy, solution_copy, N_DEPOTS, N_SITES, agent_tours, ECBS_WEIGHT)

            # Enter only if both valid
            if did_replan_indiv && did_replan_collec

                println("Replanning done successfully!")

                # Now look at times and makespans
                msp_one = maximum([s.cost for s in solution])
                msp_all = maximum([s.cost for s in new_soln_collec])

                push!(rp_one_times, el_time_indiv)
                push!(rp_one_msps, msp_one)

                push!(rp_all_times, el_time_collec)
                push!(rp_all_msps, msp_all)

                # Increment trial
                TRIAL = TRIAL + 1

            else

                println("No replan! Indiv - $(did_replan_indiv); Collec - $(did_replan_collec)")

            end
        end
    end

    replan_results["results"] = Dict("replan_one_times" => rp_one_times,
                                     "replan_one_makespans" => rp_one_msps,
                                     "replan_all_times" => rp_all_times,
                                     "replan_all_makespans" => rp_all_msps)

    out_file = string(out_file_pref, "_", N_DEPOTS, "deps_", N_AGENTS, "_agts.json")
    open(out_file, "w") do f
        JSON.print(f, replan_results, 2)
    end
end

main()


text:
using JGUI
using Base.Test

w = window()
f = vbox(w); push!(w, f)

b = button(f, "click me")
push!(f, b); raise(w)

b[:icon] = :ok                  # set icon
@test b[:value] == "click me"   # getValue

b[:value] = "New value"         # setValue
@test b[:value] == "New value"

## Callbacks
ctr = [1]
connect(b, "clicked", () -> ctr[1] = ctr[1] + 1)

## invoke
notify(b.model, "clicked")
@test ctr[1] ==  2


text:
abstract type AbstractMesh end

struct Cheb <: AbstractMesh end
struct FiniteDifferenceOrder2  <: AbstractMesh end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Ocean Color Classifications
#
# Here we apply simple & standard classifiers to `ocean color data` over a wide region.
#
# _Credit: sample data file was provided by T. Jackson from [Plymouth Marine Lab](https://www.pml.ac.uk)_
#
# <img src="../figs/cbiomes-01.png" alt="Drawing" style="height: 100px;"/>

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Activate packages for later use
#
# It is assumed that listed packages have aleary been installed using `julia`'s package manager (documentation available [here](https://docs.julialang.org/en/)).
# -

using OceanColorData, Plots, NCDatasets

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Observed wavebands
#
# Currently, the `OC-CCI` [satellite data set](https://esa-oceancolour-cci.org) provides remotely sensed reflectance at 6 wavelengths (`wv_cci` in `nm`)
# -

wv_cci=[412, 443, 490, 510, 555, 670];

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Optical classification using reflectances
#
# `Fuzzy logic` classifiers defined in [Moore et al 2009](https://doi.org/10.1016/j.rse.2009.07.016) and [Jackson et al 2017](http://dx.doi.org/10.1016/j.rse.2017.03.036) can be used to assign optical class memberships from an `Rrs` vector. While Moore et al define `n=8` classes using an in-situ database, Jackson et al instead define `n=14` classes using a satellite database. The latter benefits from better data coverage across all of the ecological provinces of the global ocean and is used in `OC-CCI`.
#
# In both cases the classifier is encoded in a mean reflectance spectra (`M[i][1:6]`) and a covariance matrix (`S[i][1:6,1:6]`) provided for each optical class (`i` in `1:n`). Class memberships are then derived by computing the squared Mahalanobis distance to each `M[i]` and passing the result to cumulative chi-squared distribution function (Equations 11 and 12 in [Moore et al 2011](https://doi.org/10.1109/36.942555)).

# + {"slideshow": {"slide_type": "subslide"}}
(M,Sinv)=Jackson2017();

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Read In Data Sample
#
# _Credit: sample data file was provided by T. Jackson from [Plymouth Marine Lab](https://www.pml.ac.uk)_

# +
dir0="../samples/";
fil=dir0*"ESACCI-OC-RRS-sample-fv4.0.nc"
ds = Dataset(fil)

Rrs_412=ds["Rrs_412"]; Rrs_443=ds["Rrs_443"]; Rrs_490=ds["Rrs_490"];
Rrs_510=ds["Rrs_510"]; Rrs_555=ds["Rrs_555"]; Rrs_670=ds["Rrs_670"];
lon=ds["lon"]; lat=ds["lat"];

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# Plot one of the wave bands as an example.

# + {"slideshow": {"slide_type": "-"}}
c=transpose(reverse(ds["Rrs_490"][:,:,1],dims=2))
heatmap(lon,reverse(lat),c,title="Rrs at 490nm")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Compute memberships
#
# Apply Classification to our 2D data sample and plot out a map.

# +
#Find points that have a full set of input
tmp=fill(false,size(Rrs_412))
for ii in eachindex(Rrs_412)
    !ismissing(Rrs_412[ii]) ? tmp[ii]=!ismissing(Rrs_443[ii].*Rrs_490[ii].*Rrs_510[ii].*Rrs_555[ii].*Rrs_670[ii]) : nothing
end
ii=findall(tmp);

#Apply classifier
mbrshp=fill(NaN,(1440,960,14))
for jj=1:length(ii);
    kk=ii[jj]
    Rrs_tmp=[Rrs_412[kk] Rrs_443[kk] Rrs_490[kk] Rrs_510[kk] Rrs_555[kk] Rrs_670[kk]]
    mbrshp[kk[1],kk[2],:]=FuzzyClassification(M,Sinv,vec(Rrs_tmp))
end
# + {"slideshow": {"slide_type": "subslide"}}
c=transpose(reverse(mbrshp[:,:,10],dims=2))
heatmap(lon,reverse(lat),c,title="Membership for class 10")

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
2############################################################################
#
#  Quadrotor.jl
#
#  Purpose:    Top level quadrotor file
#
#  Author:     Edward Daughtery
#
#  Date:       05 August 2016
#
#
############################################################################
# Mass Properties
#    This function provides the mass properties for a simple quadrotor
function Mass_Props()
    mass = 0.5;  # Kg

    II = [
        2.32  0.0    0.0
        0.0   2.32   0.0
        0.0   0.0    4.0
    ]*10.0^-3; ## kg-m^2

    cg_vec_body = [0;0;0.0];

    massprops = MassProps(mass,II,cg_vec_body);

    return massprops;
end

############################################################################
# Initialize Quadrotor Structure
function QuadInit( time::Float64, dt::Float64, position_init::WGS_Pos, vel_init::Velocity_States, euler_init::Euler_Angles, rate_init::Angle_Rates )

    # Mass Properties
    massprop = Mass_Props();

    # Airframe Initialization
    quad_airframe = Airframe_Initialize( time, time, dt, position, velocity, euler, angle_rates, massprop );

    # Navigation Message Intialization
    zvec = [0.0;0.0;0.0];
    nav_state = NavStates( time, position, velocity, rate_init, zvec, zvec, euler, quad_airframe.quat_bi );

    # Motor Initialization
    wn = 220.0;
    zeta = 0.707;
    rpm_pos_lim = 100000.0*2.0*pi;
    rpm_neg_lim = 0.0;
    rpm_rate_pos_lim = 250.0;
    rpm_rate_neg_lim = -250.0;

    motor1 = Initialize_Act(dt, wn, zeta,
                            rpm_pos_lim,
                            rpm_neg_lim,
                            rpm_rate_pos_lim,
                            rpm_rate_neg_lim
                            );

    motor2 = Initialize_Act(dt, wn, zeta,
                            rpm_pos_lim,
                            rpm_neg_lim,
                            rpm_rate_pos_lim,
                            rpm_rate_neg_lim
                            );

    motor3 = Initialize_Act(dt, wn, zeta,
                            rpm_pos_lim,
                            rpm_neg_lim,
                            rpm_rate_pos_lim,
                            rpm_rate_neg_lim
                            );

    motor4 = Initialize_Act(dt, wn, zeta,
                            rpm_pos_lim,
                            rpm_neg_lim,
                            rpm_rate_pos_lim,
                            rpm_rate_neg_lim
                            );

    # Motor Array
    motor_array = Motors(motor1, motor2, motor3, motor4);

    # Initial Motor Commands
    zz = 20.0;
    motorCmd = MotorCmd(zz, zz, zz, zz);

    ## Software Initialization
    swMsg = Initialize_SW(time, dt, nav_state);

    # Creating Quadrotor Structure
    quad_rotor = Quadrotor( time, dt, quad_airframe, motor_array, motorCmd, swMsg);

    return quad_rotor;
end

############################################################################
# Quadrotor Update
function QuadUpdate(quad::Quadrotor)

    quad_rotor = quad;

    ## Extracting Motor Commands
    motorcmd = quad_rotor.motorCmd;

    ## Taking Measurements
    time = quad_rotor.time = quad_rotor.dt;
    quad_rotor.time = time;
    position = quad_rotor.airframe.position;
    velocity = quad_rotor.airframe.velocity;
    rate = quad_rotor.airframe.angle_rates;
    rate_meas = quad_rotor.airframe.angle_rates.wb_bi_body;
    accel_meas = quad_rotor.airframe.accel_meas;
    euler = quad_rotor.airframe.euler;
    quat_bi = quad_rotor.airframe.quat_bi;
    nav_state = NavStates( time, position, velocity, rate, rate_meas, accel_meas, euler, quat_bi );

    # Parameter Update
    motor_array = quad.motors;
    motor1 = motor_array.motor1;
    motor2 = motor_array.motor2;
    motor3 = motor_array.motor3;
    motor4 = motor_array.motor4;

    # Software Update
    swMsg = Update_SW(motorcmd, nav_state, quad_rotor.software);
    quad_rotor.software = swMsg;
    quad_rotor.motorCmd = quad_rotor.software.motorcmd;

    # Actuator Update
    motor_cmds = swMsg.motorcmd;
    motor1 = Update_Act( motor1, motor_cmds.motor1_cmd );
    motor2 = Update_Act( motor2, motor_cmds.motor2_cmd );
    motor3 = Update_Act( motor3, motor_cmds.motor3_cmd );
    motor4 = Update_Act( motor4, motor_cmds.motor4_cmd );
    motor_array.motor1 = motor1;
    motor_array.motor2 = motor2;
    motor_array.motor3 = motor3;
    motor_array.motor4 = motor4;
    quad_rotor.motors = motor_array;

    # Quadrotor Forces and Moments
    QuadForceArray = Quad_Forces( motor_array );
    ftot_rotor_body = QuadForceArray.f1body +
        QuadForceArray.f2body +
        QuadForceArray.f3body +
        QuadForceArray.f4body;
    mtot_rotor_body = QuadForceArray.m1body +
        QuadForceArray.m2body +
        QuadForceArray.m3body +
        QuadForceArray.m4body;

    # Airframe Update
    quad_rotor.airframe = Airframe_Update(quad_rotor.airframe, ftot_rotor_body, mtot_rotor_body );

    return quad_rotor;
end


text:
# Basic experimentation with SparseCholeskyCovarSampler indicated
# that it is hard to not run into exceptions such as non-definite matrices etc
# when the dimensionality increases. Since high dimensionality is the motivation
# for using a sparse variant I have not experimented further with this...

# This is not a general solution but it works for specifically my use case and
# when the operator is *
function Base.broadcast(op, v::Array{Tv,2}, A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
  I, J = findn(A)
  V = zeros(nnz(A))
  vn, vm = size(v)
  if vn >= 1 && vm == 1
    for l in 1:nnz(A)
      row = I[l]
      V[l] = op(v[row], A.nzval[l])
    end
  elseif vn == 1 && vm >= 1
    for l in 1:nnz(A)
      col = J[l]
      V[l] = op(v[col], A.nzval[l])
    end
  else
    throw(ArgumentError("invalid dimensions"))
  end
  sparse(I, J, V)
end

mutable struct SparseCholeskyCovarSampler <: CovarianceMatrixSampler
  C::SparseMatrixCSC{Float64,Int}
  sqrtC::SparseMatrixCSC{Float64,Int}

  SparseCholeskyCovarSampler(n) = begin
    new(speye(n,n), speye(n,n))
  end
end

function update_covariance_matrix!(cms::SparseCholeskyCovarSampler, delta, a)
  C = a * cms.C + (1 - a) * delta
  cms.C = C # triu(C) + triu(C,1)' # Ensure C is symmetric. Should not be needed, investigate...
end

function decompose!(cms::SparseCholeskyCovarSampler)
  try
    #println("droptol!(C)")
    #cms.C = Base.droptol!(cms.C, 1e-4)
    #println("nnz(C) = $(nnz(cms.C)), min = $(minimum(cms.C.nzval)), max = $(maximum(cms.C.nzval))")
    #println("cholfact(C)")
    t = cholfact(cms.C)
    #println("sparse(t)")
    s = sparse(t)
    #println("s'")
    sqrtC = s'
    #println("droptol!(sqrtC)")
    #cms.sqrtC = Base.droptol!(sqrtC, 1e-4)
    #show(cms.sqrtC)
  catch error
    # We don't update if there is some problem
    println("ERROR: ", error)
  end
end

function multivariate_normal_sample(cms::SparseCholeskyCovarSampler, n, m)
  #if n*m > 100
    # We select a random density in [0.10, 0.60]
  #  density = 0.10 + 0.60 * rand()
  #else
  #  density = 0.90
  #end
  cms.sqrtC * randn(n, m)
end


text:
function cstrt_messine_ngon(n)
if mod(log2(n),1)==0 && n>=4
    # initialization
    d = pi/4 - asin(cos(pi/n)/sqrt(2));
    x = zeros(Int(n/2));
    for i = 1:Int(n/2)
        x[i] = pi/n + (-1)^i*d;
    end
    # construction
    a = zeros(n-1); b = zeros(n-1);
    a[Int(n/2-1)] = sin(x[1]); a[Int(n/2+1)] = -a[Int(n/2-1)];
    b[Int(n/2-1)] = cos(x[1]); b[Int(n/2+1)] = b[Int(n/2-1)];
    for i = 2:Int(n/2-1)
        a[Int(mod(i*(n-1)/2-mod(i,2)/2,n))] = a[Int(mod((i-1)*(n-1)/2-mod(i-1,2)/2,n))] - (-1)^i*sin(sum(x[1:i]));
        a[Int(n-mod(i*(n-1)/2-mod(i,2)/2,n))] = -a[Int(mod(i*(n-1)/2-mod(i,2)/2,n))];
        b[Int(mod(i*(n-1)/2-mod(i,2)/2,n))] = b[Int(mod((i-1)*(n-1)/2-mod(i-1,2)/2,n))] - (-1)^i*cos(sum(x[1:i]));
        b[Int(n-mod(i*(n-1)/2-mod(i,2)/2,n))] = b[Int(mod(i*(n-1)/2-mod(i,2)/2,n))];
    end
    a[Int(n/2)] = 0; b[Int(n/2)] = 1;
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
using OpenSAFT, Test

@testset "SAFT methods single phase" begin
    test_system = PCSAFT(["ethanol"])
    @test OpenSAFT.volume(test_system, 1E5, 298) ≈ 5.9069454773905654e-5 rtol = 0.01 #returns incorrect value
    @test OpenSAFT.sat_pure(test_system, 298)[1] ≈ 7904.2246894075415 rtol = 0.01
    @test OpenSAFT.enthalpy_vap(test_system, 298) ≈ 41720.97771100548 rtol = 0.01
    @test_broken OpenSAFT.crit_pure(test_system)[1] ≈ 1 rtol = 1E10 #T_scale not defined
    @test OpenSAFT.pressure(test_system, 1, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.entropy(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.chemical_potential(test_system, 1E5, 298)[1] ≈ 1 rtol = 1E10
    @test OpenSAFT.internal_energy(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.enthalpy(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.gibbs_free_energy(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.helmholtz_free_energy(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.isochoric_heat_capacity(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.isobaric_heat_capacity(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.isothermal_compressibility(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.isentropic_compressibility(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.speed_of_sound(test_system, 1E5, 298) ≈ 1 rtol = 1E10 #requires that the model has Mr
    @test OpenSAFT.isobaric_expansivity(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.joule_thomson_coefficient(test_system, 1E5, 298) ≈ 1 rtol = 1E10
    @test OpenSAFT.second_virial_coeff(test_system, 298) ≈ 1 rtol = 1E10
end


text:
## ---------------- msh ---------------
function importmsh(filename::String)
    fid = open(filename,"r")
    readline(fid)
    line = readline(fid)
    v_,f_,d_ = split(line," ")
    version = parse(Float64,v_)
    filetype = parse(Int,f_)
    datasize = parse(Int,d_)
    readline(fid)
    if version == 4.1
        elements,nodes = import_msh_4(fid)
    elseif version == 2.2
        elements,nodes = import_msh_2(fid)
    else
        println("Version does not match!")
    end
    return elements, nodes
end

function import_msh_4(fid::IO) end

function import_msh_2(fid::IO)
    etype = Dict(1=>:Seg2,2=>:Tri3,3=>:Quad,15=>:Poi1)
    nodes = Dict{Symbol,Vector{Float64}}()
    elements = Dict{String,Any}()
    physicalnames = Dict{Int,String}()
    for line in eachline(fid)
        if line == "\$PhysicalNames"
            numPhysicalNames = parse(Int,readline(fid))
            physicalnames=>Dict{Int,String}()
            for i in 1:numPhysicalNames
                line = readline(fid)
                d_,p_,n_ = split(line," ")
                dimension = parse(Int,d_)
                physicalTag = parse(Int,p_)
                name = strip(n_,'\"')
                physicalnames[physicalTag] = name
            end
            readline(fid)
        elseif line == "\$Nodes"
            line = readline(fid)
            nₚ = parse(Int,line)
            x = zeros(nₚ)
            y = zeros(nₚ)
            z = zeros(nₚ)
            for i in 1:nₚ
                line = readline(fid)
                t_,x_,y_,z_ = split(line," ")
                tag = parse(Int,t_)
                x[i] = parse(Float64,x_)
                y[i] = parse(Float64,y_)
                z[i] = parse(Float64,z_)
            end
            nodes[:x] = x
            nodes[:y] = y
            nodes[:z] = z
            readline(fid)
        elseif line == "\$Elements"
            line = readline(fid)
            nₑ = parse(Int,line)
            for i in 1:nₑ
                line = readline(fid)
                elmN_,elmT_,numT_,phyT_,elmE_,l_... = split(line," ")
                elmNumber = parse(Int,elmN_)
                elmType = parse(Int,elmT_)
                numTag = parse(Int,numT_)
                phyTag = parse(Int,phyT_)
                elmEntary = parse(Int,elmE_)
                nodeList = parse.(Int,l_)
                name = physicalnames[phyTag]
                type = etype[elmType]
                haskey(elements,name) ? push!(elements[name],Element{type}(nodes,nodeList...)) : elements[name]=Element{type}[Element{type}(nodes,nodeList...)]
            end
            return elements, nodes
        end
    end
end

function importmsh(filename::String,config::Dict{Any,Any})
    elms, nodes = importmsh(filename)
    elements = Dict{String,Any}()
    for (name,cfg) in config
        Type = eval(Meta.parse(cfg["𝓒"]["type"]))
        if haskey(cfg,"𝓖")
            QType = Meta.parse(cfg["𝓖"]["type"])
            if haskey(cfg["𝓖"],"tag")
                elms_𝓖 = Type<:ReproducingKernel{SNode} ? ReproducingKernel{SNode,:Linear1D,:□,:CubicSpline}(elms[cfg["𝓖"]["tag"]]) : elms[cfg["𝓖"]["tag"]]
                elms_𝓒 = elms[cfg["𝓒"]["tag"]]∩elms_𝓖
                set𝓖!(elms_𝓖,QType)
                elems = Type(elms_𝓒,elms_𝓖)
                elements[name] = elems
            else
                elems = Type(elms[cfg["𝓒"]["tag"]])
                set𝓖!(elems,QType)
                elements[name] = elems
            end
            if haskey(cfg["𝓖"],"𝝭")
                ss = cfg["𝓖"]["𝝭"]
                ss = [Meta.parse(s) for s in ss]
                set_storage_𝝭!(elements[name],ss...)
            end
        else
            elements[name] = Type(elms[cfg["𝓒"]["tag"]])
        end
        if haskey(cfg,"𝗠")
            ss = Meta.parse.(cfg["𝗠"])
            set_memory_𝗠!(elements[name],ss...)
        end
        if haskey(cfg,"𝝭")
            ss = Meta.parse.(cfg["𝝭"])
            set_memory_𝝭!(elements[name],ss...)
        end
    end
    return elements,nodes
end


text:
struct CFG
  graph::Vector{Vector{Int}}
end

function CFG(ir::IR)
  graph = Vector{Int}[]
  for b in blocks(ir)
    push!(graph, Int[])
    for c in successors(b)
      push!(graph[end], c.id)
    end
  end
  return CFG(graph)
end

Base.:(==)(a::CFG, b::CFG) = a.graph == b.graph

Base.length(c::CFG) = length(c.graph)
Base.getindex(c::CFG, b::Integer) = c.graph[b]

function Base.transpose(cfg::CFG)
  cfg′ = CFG([Int[] for _ = 1:length(cfg)])
  for i = 1:length(cfg), j in cfg[i]
    push!(cfg′[j], i)
  end
  return cfg′
end

Base.adjoint(cfg::CFG) = transpose(cfg)

function definitions(b::Block)
  defs = [Variable(i) for i = 1:length(b.ir.defs) if b.ir.defs[i][1] == b.id]
  append!(defs, arguments(b))
end

function usages(b::Block)
  uses = Set{Variable}()
  prewalk(b) do x
    x isa Variable && push!(uses, x)
    return x
  end
  return uses
end

function dominators(cfg; entry = 1)
  preds = cfg'
  blocks = [1:length(cfg.graph);]
  doms = Dict(b => Set(blocks) for b in blocks)
  while !isempty(blocks)
    b = popfirst!(blocks)
    # We currently special case the first block here,
    # since Julia sometimes creates blocks with no predecessors,
    # which otherwise throw off the analysis.
    ds = isempty(preds[b]) ? Set([b, entry]) :
      push!(intersect([doms[c] for c in preds[b]]...), b)
    if ds != doms[b]
      doms[b] = ds
      for c in cfg[b]
        c in blocks || push!(blocks, c)
      end
    end
  end
  return doms
end

function domtree(cfg; entry = 1)
  doms = dominators(cfg, entry = entry)
  doms = Dict(b => filter(c -> b != c && b in doms[c], 1:length(cfg)) for b in 1:length(cfg))
  children(b) = filter(c -> !(c in union(map(c -> doms[c], doms[b])...)), doms[b])
  tree(b) = Pair{Int,Any}(b,tree.(children(b)))
  tree(entry)
end

function idoms(cfg; entry = 1)
  ds = zeros(Int, length(cfg))
  _idoms((a, bs)) = foreach(((b, cs),) -> (ds[b] = a; _idoms(b=>cs)), bs)
  _idoms(domtree(cfg, entry = entry))
  return ds
end

function domorder(ir, start = 1; full = false)
  tree = domtree(CFG(ir), entry = start)
  flatten((b,cs)) = vcat(b, flatten.(cs)...)
  tree = flatten(tree)
  if full
    for b = 1:length(ir.blocks)
      b in tree || push!(tree, b)
    end
  end
  return tree
end

domorder!(ir::IR, start = 1) = permute!(ir, domorder(ir, start, full = true))

function verify(ir::IR)
  @assert issorted(domorder(ir)) "Blocks are not in domtree order."
  doms = dominators(ir)
  # TODO check definitions within a block
  for (b, ds) in doms
    defs = union(definitions.(ds)...)
    for x in setdiff(usages(b), defs)
      error("Variable $x in block $(b.id) is not defined.")
    end
  end
  return
end

function renumber(ir)
  p = Pipe(ir)
  for (v, st) in p
    ex = st.expr
    if isbits(ex) # Trivial expressions can be inlined
      delete!(p, v)
      substitute!(p, v, substitute(p, ex))
    end
  end
  return finish(p)
end

function merge_returns!(ir)
  bs = [b for b in blocks(ir) if isreturn(b)]
  length(bs) == 1 && bs[1] == blocks(ir)[end] && return ir
  block!(ir)
  ret = argument!(blocks(ir)[end])
  return!(ir, ret)
  for b in bs
    branches(b)[end] = branch(length(ir.blocks), arguments(branches(b)[end])[1])
  end
  return ir
end

function expand!(ir::IR)
  worklist = blocks(ir)
  spats = Dict(b => Dict() for b in blocks(ir))
  while !isempty(worklist)
    b = pop!(worklist)
    b.id == 1 && continue
    defs = definitions(b)
    uses = usages(b)
    for v in setdiff(uses, defs)
      haskey(spats[b], v) && continue
      spats[b][v] = argument!(b, v)
      for c in predecessors(b)
        c in worklist || push!(worklist, c)
      end
    end
    ir.blocks[b.id] = prewalk(x -> get(spats[b], x, x), ir.blocks[b.id])
  end
  return ir
end

function prune!(ir::IR)
  worklist = blocks(ir)
  while !isempty(worklist)
    b = popfirst!(worklist)
    isempty(arguments(b)) && continue
    brs = filter(br -> br.block == b.id, [br for a in blocks(ir) for br in branches(a)])
    isempty(brs) && continue
    inputs = [setdiff(in, (a,)) for (a, in) in zip(arguments(b), zip(arguments.(brs)...))]
    del = findall(x -> length(x) == 1, inputs)
    rename = Dict(zip(arguments(b)[del], first.(inputs[del])))
    if !isempty(rename)
      prewalk!(x -> get(rename, x, x), ir)
      deletearg!(b, del)
      for c in successors(b)
        c in worklist || push!(worklist, c)
      end
    end
  end
  return ir
end

function ssa!(ir::IR)
  current = 1
  defs = Dict(b => Dict{Slot,Any}() for b in 1:length(ir.blocks))
  todo = Dict(b => Dict{Int,Vector{Slot}}() for b in 1:length(ir.blocks))
  function reaching(b, v)
    haskey(defs[b.id], v) && return defs[b.id][v]
    b.id == 1 && return undef
    x = defs[b.id][v] = argument!(b, insert = false)
    for pred in predecessors(b)
      if pred.id < current
        for br in branches(pred, b)
          push!(br.args, reaching(pred, v))
        end
      else
        push!(get!(todo[pred.id], b.id, Slot[]), v)
      end
    end
    return x
  end
  for b in blocks(ir)
    current = b.id
    rename(ex) = prewalk(x -> x isa Slot ? reaching(b, x) : x, ex)
    for (v, st) in b
      ex = st.expr
      if isexpr(ex, :(=)) && ex.args[1] isa Slot
        defs[b.id][ex.args[1]] = rename(ex.args[2])
        delete!(ir, v)
      else
        ir[v] = rename(ex)
      end
    end
    for i = 1:length(BasicBlock(b).branches)
      BasicBlock(b).branches[i] = rename(BasicBlock(b).branches[i])
    end
    for (succ, ss) in todo[b.id], br in branches(b, succ)
      append!(br.args, [reaching(b, v) for v in ss])
    end
  end
  return ir
end

function reachable_blocks(cfg::CFG)
  bs = Int[]
  reaches(b) = b in bs || (push!(bs, b); reaches.(cfg[b]))
  reaches(1)
  return bs
end

function trimblocks!(ir::IR)
  for b in sort(setdiff(1:length(blocks(ir)), reachable_blocks(CFG(ir))), rev = true)
    deleteblock!(ir, b)
  end
  return ir
end

function inlineable!(ir)
  pushfirst!(ir, Expr(:meta, :inline))
  return ir
end

function log!(ir, msg)
  pushfirst!(ir, xcall(Core, :println, msg))
  return ir
end

totype(T::Type) = T

if isdefined(Core.Compiler, :PartialStruct)
  totype(T::Core.Compiler.PartialStruct) = T.typ
end

function pis!(ir::IR)
  for (v, st) in ir
    ex = st.expr
    ex isa PiNode || continue
    ir[v] = xcall(Core, :typeassert, ex.val, totype(ex.typ))
  end
  return ir
end


text:
function ser(x::Array{Float64,N}, thresh::Float64=0.01) where N
  @dprint "computing signal enhancement ratios"
  @assert thresh > 0.0
  dims = size(x)
  nt = dims[1]
  n = prod(dims[2:end])
  SER = zeros(n)
  x = reshape(x, (nt, n))
  thresh *= maximum(x[:])
  for k in 1:n
    SER[k] = sum(x[1:3,k]) > thresh ? sum(x[end-3:end,k]) / sum(x[1:3,k]) : 0.0
  end
  reshape(SER, dims[2:end])
end

function r1eff(S::Array{Float64,M}, R10::Array{Float64,N}, TR::Float64, flip::Float64) where {M,N}
  @dprint "converting DCE signal to effective R1"
  @assert 0.0 < flip "flip angle must be positive"
  @assert 0.0 < TR && TR < 1.0 "TR must be in units of ms"
  dims = size(S)
  nt = dims[1]
  n = prod(dims[2:end])
  S = reshape(S, (nt, n))
  S0 = mean(S[1:2,:], dims=1)
  A = copy(S)
  R1 = similar(S)
  for k = 1:n
    E0 = exp(-R10[k] * TR)
    A[:,k] = A[:,k] / S0[k] # normalize by pre-contrast signal
    for t in 1:nt
      E = (1.0 - A[t,k] + A[t,k]*E0 - E0*cos(flip)) /
        (1.0 - A[t,k]*cos(flip) + A[t,k]*E0.*cos(flip) - E0*cos(flip))
      R1[t,k] = E > 0.0 ? (-1.0 / TR) * log(E) : 0.0
    end
  end
  reshape(R1, dims)
end

function tissueconc(R1::Array{Float64,M}, R10::Array{Float64,N}, r1::Float64) where {M,N}
  @dprint "converting effective R1 to tracer tissue concentration Ct"
  @assert r1 > 0.0
  dims = size(R1)
  nt = dims[1]
  xidxs = find(R10)
  n = prod(dims[2:end])
  R1 = reshape(R1, (nt, n))
  Ct = similar(R1)
  for x in xidxs, t in 1:nt
    Ct[t,x] = R1[t,x] > 0.0 ? (R1[t,x] - R10[x]) / r1 : 0.0
  end
  R1 = reshape(R1, dims)
  reshape(Ct, dims)
end

function fitr1(x, flip_angles::Vector{Float64}, TR::Float64,
               resid_thresh::Float64=0.01)
  @dprint "fitting R1 relaxation rate to multi-flip data"
  sizein = size(x)
  n = prod(sizein[2:end])
  nangles = sizein[1]
  @assert nangles == length(flip_angles)
  x = reshape(x, (nangles, n))
  p0 = [maximum(x), 1.0]
  model(x,p) = spgreqn(x, p, TR)
  idxs = find(mean(x, dims=1) .> 0.1*maximum(x))

  params, resid = nlsfit(model, x, idxs, flip_angles, p0)
  S0 = reshape(params[1,:], sizein[2:end])
  R10 = reshape(params[2,:], sizein[2:end])
  (R10, S0, resid)
end

function fitdce(Ct::Array{Float64,M}, mask::BitArray{N}, t::Vector{Float64},
                   Cp::Vector{Float64}; models=[2], residthresh=1.0, Ktmax=5.0) where {M,N}
  @dprint "fitting DCE data"
  sizein = size(Ct)
  n = prod(sizein[2:end])
  nt = sizein[1]
  @assert nt == length(t)
  Ct = reshape(Ct, (nt, n))
  idxs = find(mask)

  nidxs = length(idxs)
  nmodels = length(models)
  @assert nmodels > 0 "at least one model must be specified"
  resid = Inf*ones(n)
  params = zeros(3, n)
  modelmap = zeros(UInt8, n)
  if 5 in models
    @dprint "attempting linearized Extended Tofts-Kety model"
    runtimeLL = @elapsed p, r = fitETM(t, Ct, Cp)
    @dprint "Fitted $nt x $nidxs points in $runtimeLL seconds"
    for k in idxs
      if r[k] <= resid[k]
        resid[k] = r[k]
        params[1,k] = p[1,k]
        params[2,k] = p[1,k] / p[2,k] # ve = Ktrans/kep (i.e. p[1]/p[2])
        params[3,k] = p[3,k]
        modelmap[k] = 5
      end
    end
  end
  if 4 in models
    @dprint "attempting linearized Standard Tofts-Kety model"
    runtimeLL = @elapsed p, r = fitTM(t, Ct, Cp)
    @dprint "Fitted $nt x $nidxs points in $runtimeLL seconds"
    for k in idxs
      if r[k] <= resid[k]
        resid[k] = r[k]
        params[1,k] = p[1,k]
        params[2,k] = p[1,k] / p[2,k] # ve = Ktrans/kep (i.e. p[1]/p[2])
        params[3,k] = 0.0
        modelmap[k] = 4
      end
    end
  end
  if 3 in models
    @dprint "attempting Extended Tofts-Kety model"
    p0 = [0.01, 0.01, 0.01]
    f3(x,p) = extendedtoftskety(x, p, Cp)
    p, r, dof = nlsfit(f3, Ct, idxs, t, p0)
    p[2,idxs] = p[1,idxs] ./ p[2,idxs]
    r = dropdims(sum(abs2, r, dims=1), dims=1) / dof
    for k in idxs
      if r[k] <= resid[k]
        resid[k] = r[k]
        params[:,k] = p[:,k]
        modelmap[k] = 3
      end
    end
  end
  if 2 in models
    @dprint "attempting Standard Tofts-Kety model"
    p0 = [0.01, 0.01]
    f2(x,p) = toftskety(x, p, Cp)
    p, r, dof = nlsfit(f2, Ct, idxs, t, p0)
    r = dropdims(sum(abs2, r, dims=1), dims=1) / dof
    for k in idxs
      if r[k] <= resid[k]
        resid[k] = r[k]
        params[1,k] = p[1,k]
        params[2,k] = p[1,k] / p[2,k]  # ve = Ktrans / kep
        params[3,k] = 0.0
        modelmap[k] = 2
      end
    end
  end
  if 1 in models
    @dprint "attempting plasma-only model"
    p0 = [0.01]
    f1(x,p) = onecompartment(x, p, Cp)
    p, r, dof = nlsfit(f1, Ct, idxs, t, p0)
    r = squeeze(sum(abs2, r, 1), 1) / dof
    for k in idxs
      if r[k] <= resid[k]
        resid[k] = r[k]
        params[3,k] = p[1,k]
        params[1:2,k] = 0.0
        modelmap[k] = 1
      end
    end
  end
  for k in idxs
    if resid[k] <= residthresh
      params[1,k] = clamp.(params[1,k], eps(), Ktmax)
      params[2:end,k] =  clamp.(params[2:end,k], eps(), 1.0)
    else
      params[:,k] .= 0.0
      mask[k] = false
    end
  end
  params = reshape(params, [size(params,1), sizein[2:end]...]...)
  resid = reshape(resid, sizein[2:end])
  modelmap = reshape(modelmap, sizein[2:end])
  (params, resid, modelmap)
end

function fitdata(opts::Dict)
  @dprint "running models"

  # option order of precedence
  # 1. function arguments
  # 2. values in "datafile"
  # 3. defaults
  def = defaults()
  infile = haskey(opts,"datafile") ? opts["datafile"] : def["datafile"]
  mat = matread(infile)
  opts = merge(def, merge(mat, opts))
  validate(opts)

  # load DCE and R1 data
  relaxivity = opts["relaxivity"]
  TR = opts["TR"]
  models = opts["models"]
  outfile = opts["outfile"]
  Cp = vec(opts["Cp"])  # colon makes sure that we get a vector

  # parallel workers are better than multithreaded BLAS for this problem
  # run julia with the '-p <n>' flag to start with n workers
  BLAS.set_num_threads(1)
  # startworkers(opts["workers"])
  # require("DCEMRI.jl")

  if haskey(opts, "R10") && haskey(opts, "S0")
    @dprint "found existing R1 map"
    R10 = opts["R10"]
    S0 = opts["S0"]
  else
    @dprint "found multi-flip data"
    t1data = opts["T1data"]
    t1flip = vec(opts["T1flip"]) * pi / 180.0
    @assert length(t1flip) == size(t1data,1) # must have one flip angle per T1 image
    R10, S0 = fitr1(t1data, t1flip, TR)
  end

  dcedata = opts["DCEdata"]
  dceflip = opts["DCEflip"] * pi / 180.0
  t = vec(opts["t"]) / 60.0  # convert time to min

  dims = size(dcedata)[2:end]
  @assert dims == size(R10) "R10 map and DCE images have different dimensions"

  # MAIN: run postprocessing steps
  SER = ser(dcedata)
  if haskey(opts, "mask") # if mask specified, use it
      mask = BitArray(opts["mask"])
  else   # use SER threshold
      mask = SER .> opts["SERcutoff"]
  end
  R1 = r1eff(dcedata, R10, TR, dceflip)
  Ct = tissueconc(R1, R10, relaxivity)
  @. Ct[!isfinite(Ct)] = 0
  params, resid, modelmap = fitdce(Ct, mask, t, Cp, models=models)

  results = Dict()
  results["t"] = t
  results["Cp"] = Cp
  results["R10"] = R10
  results["S0"] = S0
  results["SER"] = SER
  results["mask"] = Array(mask)
  results["models"] = models
  results["modelmap"] = modelmap
  results["R1"] = R1
  results["Ct"] = Ct
  params = reshape(params, (size(params,1), prod(dims)))
  results["Kt"] = reshape(params[1,:], dims)
  results["ve"] = reshape(params[2,:], dims)
  results["vp"] = reshape(params[3,:], dims)
  results["resid"] = resid
  if opts["save"]
    @dprint "saving results to $outfile"
    matwrite(outfile, results)
  end
  results
end
fitdata(filename::AbstractString) = fitdata(datafile=filename) # point to MAT file
fitdata(; kwargs...) = fitdata(kwargs2dict(kwargs))


text:
"""
Generalised reserve loss to track carbon.

Loss is distributed between general and C and N reserves by the fraction θE
"""
@inline reserve_loss!(o, loss) = nothing #reserve_loss!(o.params, o, loss)
@inline reserve_loss!(::HasCNE, o, loss) = begin
    ee = loss * θ # fraction of loss from E reserve
    ecn = loss - ee # fraction on loss from C and N reserves
    ec = ecn/y_E_EC(o)
    en = ecn/y_E_EN(o)
    o.J1[:C,:los] += ec + en + ee
    o.J1[:N,:los] += (ec, en, ee) ⋅ (n_N_EC(o), n_N_EN(o), n_N_E(o))
    nothing
end
@inline reserve_loss!(::HasCN, o, loss) = begin
    ec = loss/y_E_EC(o)
    en = loss/y_E_EN(o)
    o.J1[:C,:los] += ec + en
    o.J1[:N,:los] += (ec, en) ⋅ (n_N_EC(o), n_N_EN(o))
    nothing
end

@inline conversion_loss!(o, loss, dest_n_N) = nothing #conversion_loss!(o.params, o, loss, dest_n_N)
@inline conversion_loss!(::HasCNE, o, loss, dest_n_N) = begin
    ecn = loss * (1 - θE(o)) # fraction on loss from C and N reserves
    ec = ecn/y_E_EC(o)
    en = ecn/y_E_EN(o)
    o.J1[:C,:los] += ec + loss * (θE(o) - 1) # + en
    o.J1[:N,:los] += (ec, en, loss * (θ - dest_n_N/n_N_E(o))) ⋅ (n_N_EC(o), n_N_EN(o), n_N_E(o))
end
@inline conversion_loss!(::HasCN, o, loss, dest_n_N) = begin
    ec = loss/y_E_EC(o)
    en = loss/y_E_EN(o)
    o.J1[:C,:los] += ec + loss # + en
    o.J1[:N,:los] += (ec, en) ⋅ (n_N_EC(o), n_N_EN(o))
end

stoich_merge_losses(Jc1, Jn1, Jc2, Jn2, JEcn, n_c, n_n, n_Ecn) = begin
    lossa = Jc1 - Jc2 + Jn1 - Jn2 - JEcn
    lossb = (Jc1 - Jc2, Jn1 - Jn2, -JEcn) ⋅ (n_c, n_n, n_Ecn)
    lossa, lossb
end


text:
using SudokuSolver
using Test

@testset "SudokuSolver.jl" begin
    s6ans = [3 7 9 8 1 6 2 5 4; 8 6 4 2 7 5 1 9 3; 5 2 1 9 4 3 6 8 7; 6 5 7 4 3 1 9 2 8;
    4 9 8 5 6 2 7 3 1; 1 3 2 7 8 9 5 4 6; 2 1 3 6 9 4 8 7 5; 7 4 5 1 2 8 3 6 9; 9 8 6 3 5 7 4 1 2]
    @test solvesudoku(SudokuSolver.PuzzleExamples[6]) == s6ans
    @test solvesudoku(SudokuSolver.PuzzleExamples[6], solver=SudokuSolver.naivesolver!) == s6ans
    @test solvesudoku(SudokuSolver.PuzzleExamples[6], solver=SudokuSolver.reordersolver!) == s6ans
    @test solvesudoku(SudokuSolver.PuzzleExamples[6], solver=SudokuSolver.prioritysolver!) == s6ans
    invalidpuzzle = "000011054800000000000000000650400000000002730000000000210000800700000300000350000"
    @test solvesudoku(invalidpuzzle) === nothing
    @test solvesudoku(invalidpuzzle, solver=SudokuSolver.naivesolver!, check=false) == SudokuSolver.loadpuzzle(invalidpuzzle)
    @test solvesudoku(invalidpuzzle, solver=SudokuSolver.reordersolver!, check=false) == SudokuSolver.loadpuzzle(invalidpuzzle)
    @test solvesudoku(invalidpuzzle, solver=SudokuSolver.prioritysolver!, check=false) |> size == (9, 9)
end


text:
__precompile__()

module FEBOL

using Distributions: Normal, cdf, MvNormal, pdf
using StatsBase: sample
import StatsBase: entropy
using ParticleFilters

export
    SearchDomain,
    theta!

export
    Vehicle,
    new_pose

export
    Sensor,
    observe,
    update!,
    centroid,
    covariance,
    entropy,
    reset!,

    RangeOnly,
    BearingOnly,
    DirOmni,
    FOV,
    FOVb,
    FOVn,
    FOV3

export
    AbstractFilter,
    DF,
    EKF,
    EIF,
    UKF,

    PF,
    particle,
    particles,
    predict,
    ParticleCollection,
    sample_state,
    sample_states

export
    Policy,
    action,

    RandomPolicy,
    SitPolicy,
    GreedyPolicy,
    GreedyP,
    CirclePolicy,
    SpinPolicy,
    GaussianMPC,
    GreedyPolicy2,

    make_action_list

export
    SimUnit,

    TerminationCondition,
    is_complete,
    StepThreshold,
    MaxNormThreshold,

    CostModel,
    get_cost,
    ConstantCost,
    MoveCost,
    MoveCost2,
    MoveAndRotateCost,

    step!,
    simulate,
    parsim

export
    Model,
    MotionModel,
    StationaryMotion,
    RandomMotion,
    ConstantMotion,
    move_target!,
    move_target

export act!
export makenorm
export print_belief
export LSInitializer, NaiveInitializer
export my_pdf
export true_bearing

const Pose = NTuple{3, Float64}
const Action = NTuple{3, Float64}
const LocTuple = NTuple{4, Float64}
const TargetTuple = NTuple{4, Float64}
const ObsBin = Int64

export TargetTuple

export Pose

# math functions
include("math.jl")

# search domain
include("motion_model.jl")
include("search_domain.jl")

# Sensing models
include("sensors/sensors.jl")

include("vehicle.jl")

# Filters
include("filters/filters.jl")

# information theory
include("information_theory/mutual_information.jl")
include("information_theory/fisher.jl")
include("information_theory/cache.jl")
include("information_theory/kl.jl")

# policies
include("policies/policies.jl")

# Simulation and Plotting
include("simulations/termination.jl")
include("simulations/costs.jl")
include("simulations/simunit.jl")
include("simulations/simulate.jl")
include("simulations/parallel.jl")

end # module


text:
using Base: show_supertypes
#push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

# add https://github.com/CliMA/ClimaCore.jl/#main
# add https://github.com/CliMA/ClimaAtmos.jl/#main

# import required modules
import ClimaCore.Geometry, LinearAlgebra, UnPack
import ClimaCore: Fields, Domains, Topologies, Meshes, DataLayouts, Operators, Geometry, Spaces

using ClimaAtmos

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

using Logging: global_logger
using TerminalLoggers: TerminalLogger

using RecursiveArrayTools

using OrdinaryDiffEq, Test, Random

global_logger(TerminalLogger())

const CI = !isnothing(get(ENV, "CI", nothing))

# general parameters
const FT = Float64

include("dummy_surface_fluxes.jl") # placeholder for SurfaceFluxes.jl

########
# Set up parameters
########

parameters = (
    # timestepping parameters
    Δt_min = 0.02, # minimum model timestep [s]
    timerange = (0.0, 100.0),  # start time and end time [s]
    odesolver = SSPRK33(), # timestepping method from DifferentialEquations.jl (used in both models here)
    nsteps_atm = 1, # no. time steps of atm before coupling
    nsteps_lnd = 1, # no. time steps of lnd before coupling
    saveat = 0.2, # interval at which to save diagnostics [s]

    # atmos domain
    zmin_atm = FT(0.0),     # height of atmos stack bottom
    zmax_atm = FT(200.0), # height of atmos domain top
    n = 30,                 # number of vertical levels

    # atmos physics (eventually will be importing/overwriting CliMAParameters defaults?)
    T_surf = 300.0, # K
    T_min_ref = 230.0, #K
    MSLP = 1e5, # mean sea level pressure [Pa]
    grav = 9.8, # gravitational constant [m /s^2]
    R_d = 287.058, # R dry (gas constant / mol mass dry air)  [J / K / kg]
    C_p = 287.058 * 1.4 / (1.4 - 1), # heat capacity at constant pressure [J / K / kg]
    C_v = 287.058 / (1.4 - 1), # heat capacity at constant volume [J / K / kg]
    R_m = 87.058, # moist R, assumed to be dry [J / K / kg]
    f = 7.29e-5,#5e-5,# Coriolis parameters [1/s]
    ν = 0.1, #0.01,# viscosity, diffusivity
    Cd = 0.01 / (2e2 / 30.0), #drag coeff
    Ch = 0.01 / (2e2 / 30.0), #thermal transfer coefficient
    ug = 1.0,
    vg = 0.0,
    d = sqrt(2.0 * 0.01 / 5e-5), #?

    # soil slab
    h_s = 100,     # depth of the modelled soil model [m]
    c_s = 800,     # specific heat for land (soil)  [J / K / kg]
    κ_s = 0.0,     # soil conductivity [W / m / K] (set to 0 for energy conservation checks)
    T_h = 280,     # temperature of soil at depth h [K]
    ρ_s = 1500,    # density for land (soil) [kg / m^3]

    # radiation
    τ = 0.9,     # atmospheric transmissivity
    α = 0.5,     # surface albedo
    σ = 5.67e-8, # Steffan Boltzmann constant [kg / s^3 / K^4]
    g_a = 0.06,    # aerodynamic conductance for heat transfer [kg / m^2 / s]
    ϵ = 0.98,    # broadband emissivity / absorptivity of the surface
    F_a = 0.0,     # downward LW flux from the atmosphere [W / m^2]
    F_sol = 1361,   # incoming solar TOA radiation [W / m^2]
    τ_d = 10,     # idealized daily cycle period [s]

    # surface fluxes
    λ = FT(0.01),#FT(1e-5)    # coupling transfer coefficient (to be replaced by the bulk formula)

    # Q: ML???
)

########
# Set up atmos domain
########

domain_atm = Domains.IntervalDomain(parameters.zmin_atm, parameters.zmax_atm, x3boundary = (:bottom, :top)) # struct
mesh_atm = Meshes.IntervalMesh(domain_atm, nelems = parameters.n) # struct, allocates face boundaries to 5,6
center_space_atm = Spaces.CenterFiniteDifferenceSpace(mesh_atm) # collection of the above, discretises space into FD and provides coords
face_space_atm = Spaces.FaceFiniteDifferenceSpace(center_space_atm)

########
# Set up inital conditions
########
function init_centers(zc, parameters)
    UnPack.@unpack T_surf, T_min_ref, grav, C_p, MSLP, R_d = parameters

    # temperature
    Γ = grav / C_p
    T = max(T_surf - Γ * zc, T_min_ref)

    # pressure
    p = MSLP * (T / T_surf)^(grav / (R_d * Γ))
    if T == T_min_ref
        z_top = (T_surf - T_min_ref) / Γ
        H_min = R_d * T_min_ref / grav
        p *= exp(-(zc - z_top) / H_min)
    end

    # potential temperature
    θ = T_surf

    # density
    ρ = p / (R_d * θ * (p / MSLP)^(R_d / C_p))

    # velocties
    u = 1.0
    v = 0.0

    return (ρ = ρ, u = u, v = v, ρθ = ρ * θ)
end

function init_faces(zf, parameters)
    return (; w = 0.0 .* zf)
end

# atmos IC state
z_centers = Fields.coordinate_field(center_space_atm)
z_faces = Fields.coordinate_field(face_space_atm)
Yc = init_centers.(z_centers, Ref(parameters))
Yf = init_faces.(z_faces, Ref(parameters))
T_atm_0 = (; Yc = Yc, Yf = Yf)

# land IC state
T_lnd_0 = [260.0] # initiates lnd progostic var

ics = (; atm = T_atm_0, lnd = T_lnd_0)

########
# Set up rhs! and BCs
########

# define model equations:
include("atmos_rhs_fully_coupled.jl")

function ∑tendencies_lnd!(dT_sfc, T_sfc, (parameters, F_sfc), t)
    """
    Slab ocean:
    ∂_t T_sfc = F_sfc + G
    """
    p = parameters
    G = 0.0 # place holder for soil dynamics

    @. dT_sfc = (F_sfc * p.C_p + G) / (p.h_s * p.ρ_s * p.c_s)
end

########
# Set up time steppers
########

# specify timestepping info
stepping = (;
    Δt_min = parameters.Δt_min,
    timerange = parameters.timerange,
    Δt_cpl = max(parameters.nsteps_atm, parameters.nsteps_lnd) * parameters.Δt_min, # period of coupling cycle [s]
    odesolver = parameters.odesolver,
    nsteps_atm = parameters.nsteps_atm,
    nsteps_lnd = parameters.nsteps_lnd,
)

# coupler comm functions which export / import / transform variables
coupler_get(x) = x
coupler_put(x) = x

# Solve the ODE operator
function coupler_solve!(stepping, ics, parameters)
    t = 0.0
    Δt_min = stepping.Δt_min
    Δt_cpl = stepping.Δt_cpl
    t_start = stepping.timerange[1]
    t_end = stepping.timerange[2]
    nsteps_atm = stepping.nsteps_atm
    nsteps_lnd = stepping.nsteps_lnd

    # init coupler fields
    coupler_F_sfc = coupler_put([0.0])
    coupler_T_lnd = coupler_put(copy(ics.lnd))

    # atmos copies of coupler variables
    atm_T_lnd = coupler_get(copy(coupler_T_lnd))
    atm_F_sfc = coupler_get(copy(coupler_F_sfc))

    # SETUP ATMOS
    # put all prognostic variable arrays into a vector and ensure that solve can partition them
    T_atm = ics.atm
    Y_atm = ArrayPartition((T_atm_0.Yc, T_atm_0.Yf, atm_F_sfc))
    prob_atm = ODEProblem(∑tendencies_atm!, Y_atm, (t_start, t_end), (parameters, atm_T_lnd))
    integ_atm = init(prob_atm, stepping.odesolver, dt = Δt_cpl / nsteps_atm, saveat = parameters.saveat)

    # land copies of coupler variables
    T_lnd = ics.lnd
    lnd_F_sfc = copy(coupler_F_sfc)

    # SETUP LAND
    prob_lnd = ODEProblem(∑tendencies_lnd!, T_lnd, (t_start, t_end), (parameters, lnd_F_sfc))
    integ_lnd = init(prob_lnd, stepping.odesolver, dt = Δt_cpl / nsteps_lnd, saveat = parameters.saveat)

    # coupler stepping
    for t in (t_start:Δt_cpl:t_end)

        ## Atmos
        # pre_atmos
        # coupler_get_atmosphere!(atmos_simulation, coupler)
        #   1) get the land state (e.g., temperature) from the coupler
        #   2) reset the flux accumulator to 0

        integ_atm.p[2] .= coupler_get(coupler_T_lnd) # get land temperature and set on atmosphere (Tland is prognostic)
        integ_atm.u.x[3] .= [0.0] # surface flux to be accumulated

        # run atmos
        # NOTE: use (t - integ_atm.t) here instead of Δt_cpl to avoid accumulating roundoff error in our timestepping.
        step!(integ_atm, t - integ_atm.t, true)

        # post_atmos
        # coupler_put_atmosphere!(atmos_simulation, coupler)
        #  1) compute the time-averaged flux to/from the land
        #  2) regrid
        #  3) store that in the coupler

        # negate sign
        coupler_F_sfc .= -coupler_put(integ_atm.u.x[3]) / Δt_cpl

        ## Land

        # pre_land
        # coupler_get_land!(land_simulation, coupler)
        #  1) get thet time-averaged flux from the coupler and save it in the `lnd_F_sfc` parameter of `integ_lnd`
        lnd_F_sfc .= coupler_get(coupler_F_sfc)

        # run land
        step!(integ_lnd, t - integ_lnd.t, true)

        # post land
        # coupler_put_land!(land_simulation, coupler)
        #  1) store required land state (e.g., temperature) in the coupler
        coupler_T_lnd .= coupler_put(integ_lnd.u) # update T_sfc
    end

    return integ_atm, integ_lnd
end

########
# Run simulation
########

integ_atm, integ_lnd = coupler_solve!(stepping, ics, parameters)
sol_atm, sol_lnd = integ_atm.sol, integ_lnd.sol

########
# Visualisation
########

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

dirname = "heat"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

# animation
# anim = Plots.@animate for u in sol_atm.u
#     Plots.plot(u.x[1], xlim=(220,280))
# end
# Plots.mp4(anim, joinpath(path, "heat.mp4"), fps = 10)

# atmos vertical profile at t=0 and t=end
t0_ρθ = parent(sol_atm.u[1].x[1])[:, 4]
tend_ρθ = parent(sol_atm.u[end].x[1])[:, 4]
t0_u = parent(sol_atm.u[1].x[1])[:, 2]
tend_u = parent(sol_atm.u[end].x[1])[:, 2]
t0_v = parent(sol_atm.u[1].x[1])[:, 3]
tend_v = parent(sol_atm.u[end].x[1])[:, 3]
z_centers = collect(1:1:length(tend_u))#parent(Fields.coordinate_field(center_space_atm))[:,1]
Plots.png(Plots.plot([t0_ρθ tend_ρθ], z_centers, labels = ["t=0" "t=end"]), joinpath(path, "T_atm_height.png"))
Plots.png(Plots.plot([t0_u tend_u], z_centers, labels = ["t=0" "t=end"]), joinpath(path, "u_atm_height.png"))
Plots.png(Plots.plot([t0_v tend_v], z_centers, labels = ["t=0" "t=end"]), joinpath(path, "v_atm_height.png"))

# time evolution
atm_sum_u_t = [sum(parent(u.x[1])[:, 4]) for u in sol_atm.u] ./ parameters.n .* parameters.zmax_atm * parameters.C_p # J / m2
lnd_sfc_u_t = [u[1] for u in sol_lnd.u] * parameters.h_s * parameters.ρ_s * parameters.c_s # J / m2

v1 = lnd_sfc_u_t .- lnd_sfc_u_t[1]
v2 = atm_sum_u_t .- atm_sum_u_t[1]
Plots.png(
    Plots.plot(sol_lnd.t, [v1 v2 v1 + v2], labels = ["lnd" "atm" "tot"]),
    joinpath(path, "energy_both_surface_time.png"),
)

# relative error
using Statistics
total = atm_sum_u_t + lnd_sfc_u_t
rel_error = (total .- total[1]) / mean(total)
Plots.png(Plots.plot(sol_lnd.t, rel_error, labels = ["tot"]), joinpath(path, "rel_error_surface_time.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig("output/$(dirname)/heat_end.png", "Heat End Simulation")

# TODO here
# - integrate ClimaCoupler.jl
# - revamp for 2D
# - integrate SurfaceFluxes.jl
# - revamp for 2D

# Questions / Comments
#

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Test, Hijack

@test true

@testset "runtests" begin
    @test true
end

# `include` within nested control-flow constructs
begin
    let cond=true
        let
            if cond
                for i=1:1
                    while i==1
                        try
                            for my_test in ["included.jl"]
                                let my_test2 = my_test
                                    include(my_test2)
                                end
                            end
                            error(":)")
                        catch
                            include("simple.jl")
                        finally
                            include("simple2.jl")
                        end
                        i+=1
                    end
                end
            end
            if false
                include("not-a-file.jl")
            else
                @test true
            end
        end
    end
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
@testset "Schottky null" begin
    τ = [0.40243+0.68413im -0.18138+0.21894im 0.24323-0.13416im 0.00403+0.05085im -0.31818+0.14383im; -0.18137+0.21894im 0.27914+1.01836im -0.09799+0.46222im -0.06566+0.60959im -0.14647+0.37006im; 0.24323-0.13416im -0.09799+0.46222im 0.16663+0.68136im -0.28606+0.02038im 0.18558-0.15061im; 0.00403+0.05085im -0.06566+0.60959im -0.28606+0.02038im 0.04136+1.40560im 0.19025+0.82885im; -0.31818+0.14384im -0.14647+0.37006im 0.18558-0.15061im 0.19025+0.82885im 0.74873+1.01168im ];
    R = RiemannMatrix(τ);
    z = zeros(5);
    @test min_even_theta_char(R) == [[1,0,0,1,0],[1,0,1,1,0]]
    @test hessian(z,R) ≈ hessian(z,τ) atol=1e-3
    @test hessian(z,R, [[1,0,0,1,0],[1,0,1,1,0]])[1,1] ≈ -2.79665+5.29764im atol=1e-3
    @test schottky_null(R) ≈ schottky_null(τ) atol=1e-3
end