Created
March 20, 2020 19:25
-
-
Save vlandau/8917ec2f23ebb3a0bc4f9cc894fde2f2 to your computer and use it in GitHub Desktop.
Omniscape.jl memory allocations for functions.jl, with 4903×4405 resistance surface
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
- 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