"""
    run_omniscape("path/to/config.ini")

Run the Omniscape algorithm using the files and options specified in a .ini file

Visit https://circuitscape.github.io/Omniscape.jl/stable for detailed
documentation
"""
function run_omniscape(path::String)
    start_time = time()
    n_threads = nthreads()
    cfg_user = parse_cfg(path)

    check_missing_args(cfg_user) && return

    cfg = init_cfg()
    update_cfg!(cfg, cfg_user)

    ## Parse commonly called integer arguments
    int_arguments = Dict{String, Int64}()

    int_arguments["block_size"] = Int64(round(parse(Float64,
                                                    cfg["block_size"])))

    check_block_size(int_arguments["block_size"]) &&
        (int_arguments["block_size"] = int_arguments["block_size"] + 1)

    int_arguments["block_radius"] = Int64((int_arguments["block_size"] - 1) / 2)
    int_arguments["radius"] = Int64(round(parse(Float64, cfg["radius"])))
    int_arguments["buffer"] = Int64(round(parse(Float64, cfg["buffer"])))
    int_arguments["n_conditions"] = Int64(round(parse(Float64, cfg["n_conditions"])))

    ## Parse other arguments
    compare_to_future = lowercase(cfg["compare_to_future"])
    precision = cfg["precision"] in SINGLE ? Float32 : Float64

    # flags
    calc_flow_potential = cfg["calc_flow_potential"] in TRUELIST
    calc_normalized_current = cfg["calc_normalized_current"] in TRUELIST
    write_raw_currmap = cfg["write_raw_currmap"] in TRUELIST
    parallelize = cfg["parallelize"] in TRUELIST
    correct_artifacts = cfg["correct_artifacts"] in TRUELIST
    source_from_resistance = cfg["source_from_resistance"] in TRUELIST
    conditional = cfg["conditional"] in TRUELIST
    mask_nodata = cfg["mask_nodata"] in TRUELIST
    resistance_file_is_conductance = cfg["resistance_file_is_conductance"] in TRUELIST
    write_as_tif = cfg["write_as_tif"] in TRUELIST
    allow_different_projections = cfg["allow_different_projections"] in TRUELIST
    reclassify = cfg["reclassify_resistance"] in TRUELIST
    write_reclassified_resistance = cfg["write_reclassified_resistance"] in TRUELIST

    if int_arguments["block_size"] == 1
        correct_artifacts = false
    end

    # other
    source_threshold = parse(Float64, cfg["source_threshold"])
    project_name = cfg["project_name"]
    r_cutoff = parse(Float64, cfg["r_cutoff"])
    file_format = write_as_tif ? "tif" : "asc"

    ## Set number of BLAS threads to 1
    BLAS.set_num_threads(1)

    ## Import sources and resistances
    resistance_raster = read_raster("$(cfg["resistance_file"])", precision)

    resistance = resistance_raster[1]
    wkt = resistance_raster[2]
    transform = resistance_raster[3]

    check_resistance_values(resistance) && return

    # Reclassify resistance layer
    if reclassify
        resistance_old = deepcopy(resistance)
        reclass_table = convert.(precision, readdlm("$(cfg["reclass_table"])"))

        for i in 1:(size(reclass_table)[1])
            resistance[coalesce.(resistance_old .== reclass_table[i, 1], false)] .= reclass_table[i, 2]
        end

        resistance_old = nothing # remove from memory
    end

    # Compute source strengths from resistance if needed
    if source_from_resistance
        source_strength = deepcopy(resistance)
        if !resistance_file_is_conductance
            source_strength = Array{Union{precision, Missing}, 2}(1 ./ source_strength)
        end
        source_strength[coalesce.(source_strength .< 1/r_cutoff, true)] .= 0.0 # handles replacing NoData with 0 as well
    else
        sources_raster = read_raster("$(cfg["source_file"])", precision)
        source_strength = sources_raster[1]

        # Check for raster alignment
        check_raster_alignment(resistance_raster, sources_raster,
                               "resistance_file", "sources_file",
                               allow_different_projections) && return

        # get rid of unneeded raster to save memory
        sources_raster = nothing

        # overwrite nodata with 0
        source_strength[ismissing.(source_strength)] .= 0.0

        # Set values < user-specified threshold to 0
        source_strength[source_strength .< source_threshold] .= 0.0
    end

    int_arguments["nrows"] = size(source_strength, 1)
    int_arguments["ncols"] = size(source_strength, 2)

    # Import conditional rasters and other conditional connectivity stuff
    if conditional
        condition1_raster = read_raster("$(cfg["condition1_file"])", precision)
        condition1 = condition1_raster[1]

        # Check for raster alignment
        check_raster_alignment(resistance_raster, condition1_raster,
                               "resistance_file", "condition1_file",
                               allow_different_projections) && return

        # get rid of unneedecheck_rasterd raster to save memory
        condition1_raster = nothing

        if compare_to_future == "1" || compare_to_future == "both"
            condition1_future_raster = read_raster("$(cfg["condition1_future_file"])", precision)
            condition1_future = condition1_future_raster[1]

            # Check for raster alignment
            check_raster_alignment(resistance_raster, condition1_future_raster,
                                   "resistance_file", "condition1_future_file",
                                   allow_different_projections) && return

            # get rid of unneeded raster to save memory
            condition1_future_raster = nothing
        else
            condition1_future = condition1
        end

        if int_arguments["n_conditions"] == 2
            condition2_raster = read_raster("$(cfg["condition2_file"])", precision)
            condition2 = condition2_raster[1]

            # Check for raster alignment
            check_raster_alignment(resistance_raster, condition2_raster,
                                   "resistance_file", "condition2_file",
                                   allow_different_projections) && return

            # get rid of unneeded raster to save memory
            condition2_raster = nothing

            if compare_to_future == "2" || compare_to_future == "both"
                condition2_future_raster = read_raster("$(cfg["condition2_future_file"])", precision)
                condition2_future = condition2_future_raster[1]

                # Check for raster alignment
                check_raster_alignment(resistance_raster, condition2_future_raster,
                                       "resistance_file", "condition2_future_file",
                                       allow_different_projections) && return

                # get rid of unneeded raster to save memory
                condition2_future_raster = nothing
            else
                condition2_future = condition2
            end

        else
            condition2 = Array{Union{Missing, precision}, 2}(undef, 1, 1)
            condition2_future = condition2
        end
    else
        condition1 = Array{Union{Missing, precision}, 2}(undef, 1, 1)
        condition2 = Array{Union{Missing, precision}, 2}(undef, 1, 1)
        condition1_future = condition1
        condition2_future = condition2
    end

    # get rid of unneeded raster to save memory
    resistance_raster = nothing

    comparison1 = cfg["comparison1"]
    comparison2 = cfg["comparison2"]
    condition1_lower = parse(Float64, cfg["condition1_lower"])
    condition2_lower = parse(Float64, cfg["condition2_lower"])
    condition1_upper = parse(Float64, cfg["condition1_upper"])
    condition2_upper = parse(Float64, cfg["condition2_upper"])

    ## Setup Circuitscape configuration
    cs_cfg_dict = init_csdict(cfg)
    cs_cfg = Circuitscape.init_config()
    Circuitscape.update!(cs_cfg, cs_cfg_dict)

    ## Calculate targets
    targets = get_targets(source_strength,
                          int_arguments,
                          precision)

    ## Circuitscape calls in loop over targets
    n_targets = size(targets, 1)

    ## Define parameters for cs
    # Get flags
    o = Circuitscape.OutputFlags(false, false,
                                 false, false,
                                 false, false,
                                 false, false)

    precision_name = precision == Float64 ? "double" : "single"
    ## Add parallel workers
    if parallelize
        println("Starting up Omniscape. Using $(n_threads) workers in parallel. Using $(precision_name) precision...")

        cum_currmap = fill(convert(precision, 0.),
                           int_arguments["nrows"],
                           int_arguments["ncols"],
                           n_threads)

        if calc_flow_potential || calc_normalized_current
            fp_cum_currmap = fill(convert(precision, 0.),
                                  int_arguments["nrows"],
                                  int_arguments["ncols"],
                                  n_threads)
        else
            # Hacky fix -- a later function needs fp_cum_currmap to be an array
            fp_cum_currmap = Array{precision, 3}(undef, 1, 1, 1)
        end
    else
        println("Starting up Omniscape. Running in serial using 1 worker. Using $(precision_name) precision...")
        cum_currmap = fill(convert(precision, 0.),
                          int_arguments["nrows"],
                          int_arguments["ncols"],
                          1)

        if calc_flow_potential || calc_normalized_current
            fp_cum_currmap = fill(convert(precision, 0.),
                                  int_arguments["nrows"],
                                  int_arguments["ncols"],
                                  1)
        else
            fp_cum_currmap = Array{precision, 3}(undef, 1, 1, 1)
        end
    end

    solver = "cg+amg"

    # n_cells = int_arguments["nrows"] * int_arguments["ncols"]
    # if n_cells <= 2000000
    #     solver = "cholmod" # TODO: "cholmod" not available in advanced mode
    # end

    flags = Circuitscape.RasterFlags(true, false, true,
                                     false, false,
                                     false, Symbol("rmvsrc"),
                                     cfg["connect_four_neighbors_only"] in TRUELIST,
                                     false, solver, o)

    if correct_artifacts
        println("Calculating block artifact correction array...")
        correction_array = calc_correction(int_arguments,
                                           cs_cfg,
                                           flags,
                                           o,
                                           conditional,
                                           condition1,
                                           condition1_future,
                                           condition2,
                                           condition2_future,
                                           comparison1,
                                           comparison2,
                                           condition1_lower,
                                           condition1_upper,
                                           condition2_lower,
                                           condition2_upper,
                                           precision)

    else
        correction_array = Array{precision, 2}(undef, 1, 1)
    end

    ## Calculate and accumulate currents on each worker
    println("Solving moving window targets...")

    ## Create progress object
    p = Progress(n_targets; dt = 0.25, barlen = 60)

    if parallelize
        parallel_batch_size = Int64(round(parse(Float64, cfg["parallel_batch_size"])))
        n_batches = Int(ceil(n_targets / parallel_batch_size))

        @threads for i in 0:(n_batches - 1)
            start_ind = parallel_batch_size * i + 1
            end_ind = min(n_targets, start_ind + parallel_batch_size - 1)

            for j in start_ind:end_ind
                solve_target!(j,
                              n_targets,
                              int_arguments,
                              targets,
                              source_strength,
                              resistance,
                              cs_cfg,
                              flags,
                              o,
                              calc_flow_potential || calc_normalized_current,
                              correct_artifacts,
                              conditional,
                              condition1,
                              condition1_future,
                              condition2,
                              condition2_future,
                              comparison1,
                              comparison2,
                              condition1_lower,
                              condition1_upper,
                              condition2_lower,
                              condition2_upper,
                              correction_array,
                              cum_currmap,
                              fp_cum_currmap,
                              precision,
                              resistance_file_is_conductance)
                next!(p)
            end
        end
    else
        for i in 1:n_targets
            solve_target!(i,
                          n_targets,
                          int_arguments,
                          targets,
                          source_strength,
                          resistance,
                          cs_cfg,
                          flags,
                          o,
                          calc_flow_potential || calc_normalized_current,
                          correct_artifacts,
                          conditional,
                          condition1,
                          condition1_future,
                          condition2,
                          condition2_future,
                          comparison1,
                          comparison2,
                          condition1_lower,
                          condition1_upper,
                          condition2_lower,
                          condition2_upper,
                          correction_array,
                          cum_currmap,
                          fp_cum_currmap,
                          precision,
                          resistance_file_is_conductance)
            next!(p)
        end
    end

    ## Set some objects to nothing to free up memory
    source_strength = nothing
    condition1 = nothing
    condition1_future = nothing
    condition2 = nothing
    condition2_future = nothing
    GC.gc()

    ## Collapse 3-dim cum current arrays to 2-dim via sum
    cum_currmap = dropdims(sum(cum_currmap, dims = 3), dims = 3)

    if calc_flow_potential || calc_normalized_current
        fp_cum_currmap = dropdims(sum(fp_cum_currmap, dims = 3), dims = 3)
    end

    ## Normalize by flow potential
    if calc_normalized_current
        normalized_cum_currmap = cum_currmap ./ fp_cum_currmap
        # replace NaNs with 0's
        normalized_cum_currmap[isnan.(normalized_cum_currmap)] .= 0
    end

    ## create new directory if project_name already exists
    dir_suffix = 1
    while isdir(string(project_name, "_$(dir_suffix)"))
        dir_suffix+=1
    end
    isdir(project_name) && (project_name = string(project_name, "_$(dir_suffix)"))

    mkpath(project_name)

    # Copy .ini file to output directory
    cp(path, "$(project_name)/config.ini")

    ## Overwrite no data
    if mask_nodata
        if calc_normalized_current
            normalized_cum_currmap[ismissing.(resistance)] .= -9999
        end
        if calc_flow_potential
            fp_cum_currmap[ismissing.(resistance)] .= -9999
        end
        cum_currmap[ismissing.(resistance)] .= -9999
    end

    # Get rid of resistance (save first if needed)
    if reclassify && write_reclassified_resistance
        resistance[ismissing.(resistance)] .= -9999
        write_raster("$(project_name)/classified_resistance",
                     convert(Array{precision, 2}, resistance),
                     wkt,
                     transform,
                     file_format)
    end
    resistance = nothing
    GC.gc()

    ## Write outputs
    if write_raw_currmap
        write_raster("$(project_name)/cum_currmap",
                     cum_currmap,
                     wkt,
                     transform,
                     file_format)
    end

    if calc_flow_potential
        write_raster("$(project_name)/flow_potential",
                     fp_cum_currmap,
                     wkt,
                     transform,
                     file_format)
    end


    if calc_normalized_current
        write_raster("$(project_name)/normalized_cum_currmap",
                     normalized_cum_currmap,
                     wkt,
                     transform,
                     file_format)
    end

    println("Done!")
    println("Time taken to complete job: $(round(time() - start_time; digits = 4)) seconds")

    println("Outputs written to $(string(pwd(),"/",project_name))")

    ## Return outputs, depending on user options
    if calc_normalized_current && !calc_flow_potential
        return cum_currmap, normalized_cum_currmap
    elseif !calc_normalized_current && calc_flow_potential
        return cum_currmap, fp_cum_currmap
    elseif calc_normalized_current && calc_flow_potential
        return cum_currmap, fp_cum_currmap, normalized_cum_currmap
    else
        return cum_currmap
    end
end

