Skip to content

Instantly share code, notes, and snippets.

@vlandau
Created March 20, 2020 19:25
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save vlandau/8917ec2f23ebb3a0bc4f9cc894fde2f2 to your computer and use it in GitHub Desktop.
Save vlandau/8917ec2f23ebb3a0bc4f9cc894fde2f2 to your computer and use it in GitHub Desktop.
Omniscape.jl memory allocations for functions.jl, with 4903×4405 resistance surface
- abstract type Data end
-
- function clip(
- A::Array{Float64, 2};
- x::Int64,
- y::Int64,
- distance::Int64
- )
-
0 dim1 = size(A)[1]
0 dim2 = size(A)[2]
-
348257488 dist = [sqrt((j - x)^2 + (i - y)^2) for i = 1:dim1, j = 1:dim2]
-
348259328 clipped = deepcopy(A)
350704560 clipped[dist .> distance] .= -9999
-
0 clipped
- end
-
-
- function get_targets(
- source_array::Array{Float64, 2},
- arguments::Dict{String, Int64}
- )
-
0 block_size = arguments["block_size"]
0 block_radius = arguments["block_radius"]
0 nrows = arguments["nrows"]
0 ncols = arguments["ncols"]
-
0 start = (block_size + 1) / 2
-
0 xs = [start:block_size:ncols;]
0 ys = [start:block_size:nrows;]
-
0 ground_points = zeros(Float64,(length(xs)*length(ys), 2))
-
- let
- c = 1
0 for i = 1:length(xs)
0 for j = 1:length(ys)
0 ground_points[c, 1] = xs[i]
0 ground_points[c, 2] = ys[j]
0 c += 1
- end
- end
- end
-
- # create column ground_points[, 3] to hold source strengths for each target
0 ground_points = cat(ground_points,
- zeros(size(ground_points)[1], 1);
- dims = 2
- )
-
- # populate ground_points[, 3] with sum of sources in block of size
- # block_size centered on each target
0 for i = 1:size(ground_points)[1]
0 xlower = Int64(ground_points[i, 1] - block_radius)
0 xupper = min(Int64(ground_points[i, 1] + block_radius), ncols)
0 ylower = Int64(ground_points[i, 2] - block_radius)
0 yupper = min(Int64(ground_points[i, 2] + block_radius), nrows)
-
0 ground_points[i, 3] = sum(source_array[ylower:yupper, xlower:xupper])
- end
-
- # get rid of ground_points with strength below 0
0 targets = ground_points[ground_points[:, 3] .> 0, 1:3]
0 targets
- end
-
- # x and y defined by targets object. Ultimately the for loop will be done by
- # iterating through rows of targets object
- function get_source(
- source_array::Array{Float64, 2},
- arguments::Dict{String, Int64},
- conditional::Bool,
- condition1_present::Array{Float64, 2},
- condition1_future::Array{Float64, 2},
- condition2_present::Array{Float64, 2},
- condition2_future::Array{Float64, 2},
- comparison1::String,
- comparison2::String,
- condition1_lower::Float64,
- condition1_upper::Float64,
- condition2_lower::Float64,
- condition2_upper::Float64;
- x::Int64,
- y::Int64,
- strength::Float64
- )
-
0 block_radius = arguments["block_radius"]
0 radius = arguments["radius"]
0 buffer = arguments["buffer"]
0 nrows = arguments["nrows"]
0 ncols = arguments["ncols"]
-
0 source_subset = clip(source_array,
- x = x,
- y = y,
- distance = radius)
-
175244400 source_subset[source_subset .== -9999] .= 0.0
- # Set any sources inside target to NoData
0 xlower = x - block_radius
0 xupper = min(x + block_radius, ncols)
0 ylower = y - block_radius
0 yupper = min(y + block_radius, nrows)
-
128 source_subset[ylower:yupper, xlower:xupper] .= 0
-
- # Extract subset for faster solve times
0 xlower_buffered = max(x - radius - buffer, 1)
0 xupper_buffered = min(x + radius + buffer, ncols)
0 ylower_buffered = max(y - radius - buffer, 1)
0 yupper_buffered = min(y + radius + buffer, nrows)
-
1455456 source_subset = source_subset[ylower_buffered:yupper_buffered,
- xlower_buffered:xupper_buffered]
-
- # allocate total current equal to target "strength", divide among sources
- # according to their source strengths
684240 source_sum = sum(source_subset[source_subset .> 0])
2673536 source_subset[source_subset .> 0] .=
- (source_subset[source_subset .> 0] * strength) / source_sum
-
0 if conditional
0 source_target_match!(source_subset,
- arguments["n_conditions"],
- condition1_present,
- condition1_future,
- condition2_present,
- condition2_future,
- comparison1,
- comparison2,
- condition1_lower,
- condition1_upper,
- condition2_lower,
- condition2_upper,
- ylower,
- yupper,
- xlower,
- xupper,
- ylower_buffered,
- yupper_buffered,
- xlower_buffered,
- xupper_buffered
- )
- end
0 source_subset
- end
-
- function source_target_match!(source_subset::Array{Float64,2},
- n_conditions::Int64,
- condition1_present::Array{Float64,2},
- condition1_future::Array{Float64,2},
- condition2_present::Array{Float64,2},
- condition2_future::Array{Float64,2},
- comparison1::String,
- comparison2::String,
- condition1_lower::Float64,
- condition1_upper::Float64,
- condition2_lower::Float64,
- condition2_upper::Float64,
- ylower::Int64,
- yupper::Int64,
- xlower::Int64,
- xupper::Int64,
- ylower_buffered::Int64,
- yupper_buffered::Int64,
- xlower_buffered::Int64,
- xupper_buffered::Int64
- )
0 con1_present_subset = condition1_present[ylower_buffered:yupper_buffered,
- xlower_buffered:xupper_buffered]
0 if comparison1 == "within"
0 value1 = median(condition1_future[ylower:yupper, xlower:xupper])
0 source_subset[((con1_present_subset .- value1) .> condition1_upper) .|
- ((con1_present_subset .- value1) .< condition1_lower)] .= 0
0 elseif comparison1 == "equals"
0 value1 = mode(condition1_future[ylower:yupper, xlower:xupper])
0 source_subset[con1_present_subset .!= value1] .= 0
- end
-
0 if n_conditions == 2
0 con2_present_subset = condition2_present[ylower_buffered:yupper_buffered,
- xlower_buffered:xupper_buffered]
0 if comparison2 == "within"
0 value2 = median(condition2_future[ylower:yupper, xlower:xupper])
0 source_subset[((con2_present_subset .- value2) .> condition2_upper) .|
- ((con2_present_subset .- value2) .< condition2_lower)] .= 0
0 elseif comparison2 == "equals"
0 value2 = mode(condition2_future[ylower:yupper, xlower:xupper])
0 source_subset[con2_present_subset .!= value2] .= 0
- end
- end
- end
-
- function get_ground(
- arguments::Dict{String, Int64};
- x::Int64,
- y::Int64
- )
-
0 radius = arguments["radius"]
0 buffer = arguments["buffer"]
0 nrows = arguments["nrows"]
0 ncols = arguments["ncols"]
-
0 xlower_buffered = Int64(max(x - radius - buffer, 1))
0 xupper_buffered = Int64(min(x + radius + buffer, ncols))
0 ylower_buffered = Int64(max(y - radius - buffer, 1))
0 yupper_buffered = Int64(min(y + radius + buffer, nrows))
-
172781840 ground = fill(0.0,
- nrows,
- ncols)
0 ground[y, x] = Inf
-
557520 output = ground[ylower_buffered:yupper_buffered,
- xlower_buffered:xupper_buffered]
0 output
- end
-
- function get_resistance(
- raw_resistance::Array{Float64, 2},
- arguments::Dict{String, Int64};
- x::Int64,
- y::Int64
- )
-
0 radius = arguments["radius"]
0 buffer = arguments["buffer"]
0 nrows = arguments["nrows"]
0 ncols = arguments["ncols"]
-
0 xlower_buffered = Int64(max(x - radius - buffer, 1))
0 xupper_buffered = Int64(min(x + radius + buffer, ncols))
0 ylower_buffered = Int64(max(y - radius - buffer, 1))
0 yupper_buffered = Int64(min(y + radius + buffer, nrows))
-
0 resistance_clipped = clip(raw_resistance,
- x = x,
- y = y,
- distance = radius + buffer)
-
1115040 resistance = 1 ./ resistance_clipped[ylower_buffered:yupper_buffered,
- xlower_buffered:xupper_buffered]
- end
-
-
- function calculate_current(
- resistance::Array{Float64, 2},
- source::Array{Float64, 2},
- ground::Array{Float64, 2},
- flags::Circuitscape.RasterFlags,
- cs_cfg::Dict{String, String}
- )
-
0 T = Float64
- V = Int64
-
- # get raster data
- cellmap = resistance
240 polymap = Matrix{V}(undef, 0, 0)
- source_map = source
- ground_map = ground
816 points_rc = (V[], V[], V[])
240 strengths = Matrix{T}(undef, 0, 0)
-
576 included_pairs = Circuitscape.IncludeExcludePairs(:undef,
- V[],
- Matrix{V}(undef,0,0))
-
0 hbmeta = Circuitscape.RasterMeta(size(cellmap)[2],
- size(cellmap)[1],
- 0.,
- 0.,
- 1.,
- -9999.,
- 0)
-
384 rasterdata = Circuitscape.RasData(cellmap,
- polymap,
- source_map,
- ground_map,
- points_rc,
- strengths,
- included_pairs,
- hbmeta)
-
- # Generate advanced data
0 data = Circuitscape.compute_advanced_data(rasterdata, flags)
-
0 G = data.G
0 nodemap = data.nodemap
0 polymap = data.polymap
0 hbmeta = data.hbmeta
0 sources = data.sources
0 grounds = data.grounds
0 finitegrounds = data.finite_grounds
0 cc = data.cc
- src = data.src
0 check_node = data.check_node
- source_map = data.source_map # Need it for one to all mode
- cellmap = data.cellmap
-
- # Flags
- is_raster = flags.is_raster
- is_alltoone = flags.is_alltoone
- is_onetoall = flags.is_onetoall
- write_v_maps = flags.outputflags.write_volt_maps
- write_c_maps = flags.outputflags.write_cur_maps
- write_cum_cur_map_only = flags.outputflags.write_cum_cur_map_only
-
2353392 volt = zeros(eltype(G), size(nodemap))
8391248 ind = findall(x -> x != 0,nodemap)
240 f_local = Vector{eltype(G)}()
- solver_called = false
240 voltages = Vector{eltype(G)}()
2353392 outvolt = Circuitscape.alloc_map(hbmeta)
2353392 outcurr = Circuitscape.alloc_map(hbmeta)
-
0 for c in cc
0 if check_node != -1 && !(check_node in c)
- continue
- end
-
- # a_local = laplacian(G[c, c])
0 a_local = G[c,c]
1712352 s_local = sources[c]
1712352 g_local = grounds[c]
-
0 if sum(s_local) == 0 || sum(g_local) == 0
- continue
- end
-
192 if finitegrounds != [-9999.]
0 f_local = finitegrounds[c]
- else
- f_local = finitegrounds
- end
-
0 voltages = Circuitscape.multiple_solver(cs_cfg,
- a_local,
- s_local,
- g_local,
- f_local)
-
0 local_nodemap = Circuitscape.construct_local_node_map(nodemap,
- c,
- polymap)
-
- solver_called = true
-
0 Circuitscape.accum_currents!(outcurr,
- voltages,
- cs_cfg,
- a_local,
- voltages,
- f_local,
- local_nodemap,
- hbmeta)
- end
-
0 outcurr
- end
-
- function solve_target!(
- i::Int64,
- n_targets::Int64,
- int_arguments::Dict{String, Int64},
- targets::Array{Float64, 2},
- sources_raw::Array{Float64, 2},
- resistance_raw::Array{Float64, 2},
- cs_cfg::Dict{String, String},
- o::Circuitscape.OutputFlags,
- calc_flow_potential::Bool,
- correct_artifacts::Bool,
- conditional::Bool,
- condition1_present::Array{Float64, 2},
- condition1_future::Array{Float64, 2},
- condition2_present::Array{Float64, 2},
- condition2_future::Array{Float64, 2},
- comparison1::String,
- comparison2::String,
- condition1_lower::Float64,
- condition1_upper::Float64,
- condition2_lower::Float64,
- condition2_upper::Float64,
- correction_array::Array{Float64, 2},
- cum_currmap::Array{Float64, 3},
- fp_cum_currmap::Array{Float64, 3}
- )
-
- ## get source
400 @info "Solving target $(i) of $(n_targets)"
0 x_coord = Int64(targets[i, 1])
0 y_coord = Int64(targets[i, 2])
0 source = get_source(sources_raw,
- int_arguments,
- conditional,
- condition1_present,
- condition1_future,
- condition2_present,
- condition2_future,
- comparison1,
- comparison2,
- condition1_lower,
- condition1_upper,
- condition2_lower,
- condition2_upper;
- x = x_coord,
- y = y_coord,
- strength = float(targets[i, 3]))
-
- ## get ground
0 ground = get_ground(int_arguments,
- x = x_coord,
- y = y_coord)
-
- ## get resistance
0 resistance = get_resistance(resistance_raw,
- int_arguments,
- x = x_coord,
- y = y_coord)
-
0 grid_size = size(source)
- n_cells = prod(grid_size)
-
- solver = "cg+amg"
-
- # if n_cells <= 2000000
- # solver = "cholmod" # FIXME: "cholmod" not available in advanced mode
- # end
-
48 flags = Circuitscape.RasterFlags(true, false, true,
- false, false,
- false, Symbol("rmvsrc"),
- false, false, solver, o)
-
- ## Run circuitscape
0 curr = calculate_current(resistance, source, ground, flags, cs_cfg)
-
- ## If normalize = True, calculate null map and normalize
0 if calc_flow_potential == true
0 @info "Calculating flow potential for target $(i) of $(n_targets)"
0 null_resistance = fill(1., grid_size)
-
0 flow_potential = calculate_current(null_resistance,
- source,
- ground,
- flags,
- cs_cfg)
- end
-
0 if correct_artifacts
898304 correction_array2 = deepcopy(correction_array)
- lowerxcut = 1
0 upperxcut = size(correction_array, 2)
- lowerycut = 1
0 upperycut = size(correction_array, 1)
-
0 if x_coord > int_arguments["ncols"] - (int_arguments["radius"] + int_arguments["buffer"])
0 upperxcut = upperxcut -
- (upperxcut - grid_size[2])
- end
-
0 if x_coord < int_arguments["radius"] + int_arguments["buffer"] + 1
0 lowerxcut = upperxcut - grid_size[2] + 1
- end
-
0 if y_coord > int_arguments["nrows"] - (int_arguments["radius"] + int_arguments["buffer"])
0 upperycut = upperycut -
- (upperycut - grid_size[1])
- end
-
0 if y_coord < int_arguments["radius"] + int_arguments["buffer"] + 1
0 lowerycut = upperycut - grid_size[1] + 1
- end
-
557520 correction_array2 = correction_array[lowerycut:upperycut,
- lowerxcut:upperxcut]
-
557520 curr = curr .* correction_array2
-
0 if calc_flow_potential
0 flow_potential = flow_potential .* correction_array2
- end
- end
-
- ## Accumulate values
0 xlower = max(x_coord - int_arguments["radius"] - int_arguments["buffer"],
- 1)
0 xupper = min(x_coord + int_arguments["radius"] + int_arguments["buffer"],
- int_arguments["ncols"])
0 ylower = max(y_coord - int_arguments["radius"] - int_arguments["buffer"],
- 1)
0 yupper = min(y_coord + int_arguments["radius"] + int_arguments["buffer"],
- int_arguments["nrows"])
-
557520 cum_currmap[ylower:yupper, xlower:xupper, threadid()] .=
- cum_currmap[ylower:yupper, xlower:xupper, threadid()] .+ curr
-
0 if calc_flow_potential == true
0 fp_cum_currmap[ylower:yupper, xlower:xupper, threadid()] .=
- fp_cum_currmap[ylower:yupper, xlower:xupper, threadid()] .+ flow_potential
- end
- end
-
- function calc_correction(
- arguments::Dict{String, Int64},
- cs_cfg::Dict{String, String},
- o,
- conditional::Bool,
- condition1_present::Array{Float64, 2},
- condition1_future::Array{Float64, 2},
- condition2_present::Array{Float64, 2},
- condition2_future::Array{Float64, 2},
- comparison1::String,
- comparison2::String,
- condition1_lower::Float64,
- condition1_upper::Float64,
- condition2_lower::Float64,
- condition2_upper::Float64
- )
-
- # This may not apply seamlessly in the case (if I add the option) that source strengths
- # are not adjusted by target weight, but stay the same according to their
- # original values. Something to keep in mind...
-
16 solver = "cg+amg"
-
- # if (arguments["radius"]+1)^2 <= 2000000
- # solver = "cholmod" # FIXME: "cholmod" not available in advanced mode
- # end
-
48 flags = Circuitscape.RasterFlags(true, false, true,
- false, false,
- false, Symbol("keepall"),
- false, false, solver, o)
-
897936 temp_source = fill(1.,
- arguments["radius"] * 2 + arguments["buffer"] * 2 + 1,
- arguments["radius"] * 2 + arguments["buffer"] * 2 + 1)
-
0 temp_source_clip = clip(temp_source,
- x = arguments["radius"] + arguments["buffer"] + 1,
- y = arguments["radius"] + arguments["buffer"] + 1,
- distance = arguments["radius"])
-
898304 source_null = deepcopy(temp_source_clip)
719536 n_sources = sum(source_null[source_null .!= -9999])
-
215760 source_null[source_null .== -9999] .= 0.0
0 source_null[arguments["radius"] + arguments["buffer"] + 1,
- arguments["radius"] + arguments["buffer"] + 1] = 0.0
719632 source_null[source_null .!= 0.0] .= 1 / (n_sources - 1)
-
0 source_block = get_source(temp_source,
- arguments,
- conditional,
- condition1_present,
- condition1_future,
- condition2_present,
- condition2_future,
- comparison1,
- comparison2,
- condition1_lower,
- condition1_upper,
- condition2_lower,
- condition2_upper,
- x = (arguments["radius"] + arguments["buffer"] + 1),
- y = (arguments["radius"] + arguments["buffer"] + 1),
- strength = float(arguments["block_size"] ^ 2))
-
0 resistance = clip(temp_source,
- x = arguments["radius"] + arguments["buffer"] + 1,
- y = arguments["radius"] + arguments["buffer"] + 1,
- distance = arguments["radius"] + arguments["buffer"])
-
897936 ground = fill(0.0,
- size(source_null))
-
0 ground[arguments["radius"] + arguments["buffer"] + 1,
- arguments["radius"] + arguments["buffer"] + 1] = Inf
-
0 block_null_current = calculate_current(resistance,
- source_block,
- ground,
- flags,
- cs_cfg)
-
0 null_current = calculate_current(resistance,
- source_null,
- ground,
- flags,
- cs_cfg)
1377936 null_current_total = fill(0.,
- arguments["radius"] * 2 + arguments["buffer"] * 2 + arguments["block_size"],
- arguments["radius"] * 2 + arguments["buffer"] * 2 + arguments["block_size"])
-
0 for i in 1:arguments["block_size"]
0 for j in 1:arguments["block_size"]
11782716192 null_current_total[i:(i + arguments["radius"] * 2 + arguments["buffer"] * 2),
- j:(j + arguments["radius"] * 2 + arguments["buffer"] * 2)] += null_current
- end
- end
-
897936 null_current_total = null_current_total[(arguments["block_radius"] + 1):(size(null_current, 1) + arguments["block_radius"]),
- (arguments["block_radius"] + 1):(size(null_current, 2) + arguments["block_radius"])]
-
215760 null_current_total[block_null_current .== 0.] .= 0
-
215760 null_current_total[null_current_total .== 0.0] .= 1.0
215760 block_null_current[block_null_current .== 0.0] .= 1.0
-
897936 correction = null_current_total ./ block_null_current
-
0 correction
- end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment