Merge pull request #167 from Circuitscape/RA/include

Fixes for include pairs bugs

.travis.yml

@@ -1,7 +1,6 @@
 # Documentation: http://docs.travis-ci.com/user/languages/julia/
 language: julia
 julia:
-  - nightly
   - 1.0
 os:
   - linux

appveyor.yml

@@ -1,10 +1,8 @@
 environment:
   matrix:
   - julia_version: 1
-  - julia_version: nightly
 
 platform:
-  - x86 # 32-bit
   - x64 # 64-bit
 
 # # Uncomment the following lines to allow failures on nightly julia

src/consts.jl

@@ -19,6 +19,18 @@ const PAIRS_AAGRID = 4
 const PAIRS_LIST = 5
 const TRUELIST = ["True", "true", "1"]
 
+const FILE_TYPE_NPY = 1
+const FILE_TYPE_AAGRID = 2
+const FILE_TYPE_TXTLIST = 3
+const FILE_TYPE_INCL_PAIRS_AAGRID = 4
+const FILE_TYPE_INCL_PAIRS = 5
+
+const FILE_HDR_GZIP = "\x1f\x8b\x08"
+const FILE_HDR_NPY = "\x93NUMPY"
+const FILE_HDR_AAGRID = "ncols"
+const FILE_HDR_INCL_PAIRS_AAGRID = "min"
+const FILE_HDR_INCL_PAIRS = "mode"
+
 # Constants for logging
 const NONE = ["NONE", "None", "none"]
 const INFO = ["INFO", "info", "Info"]

src/core.jl

@@ -55,9 +55,6 @@ function single_ground_all_pairs(data::GraphData{T,V}, flags, cfg, log = true) w
         amg_solver_path(data, flags, cfg, log)
     else
         csinfo("Solver used: CHOLMOD")
-        if eltype(data.G) == Float32
-            cswarn("CHOLMOD solver mode works only in double precision")
-        end
         bs = parse(Int, cfg["cholmod_batch_size"])
         _cholmod_solver_path(data, flags, cfg, log, bs)
     end
@@ -105,7 +102,8 @@ function amg_solver_path(data::GraphData{T,V}, flags, cfg, log)::Matrix{T} where
 
     get_shortcut_resistances = false
     if is_raster && !write_volt_maps && !write_cur_maps && 
-            !write_cum_cur_map_only && !write_max_cur_maps
+            !write_cum_cur_map_only && !write_max_cur_maps && 
+            isempty(exclude)
         get_shortcut_resistances = true
         csinfo("Triggering resistance calculation shortcut")
         num, d = get_num_pairs_shortcut(cc, points, exclude)
@@ -170,45 +168,43 @@ function amg_solver_path(data::GraphData{T,V}, flags, cfg, log)::Matrix{T} where
                 comp_j = V(comp_j)
                 J = findall(x -> x == pj, points)
 
-                # Forget excluded pairs
-                ex = false
-                for c_i in I, c_j in J
-                    if (c_i, c_j) in exclude
-                        ex = true
-                        break
-                    end
-                end
-                ex && continue
-
                 if pi == pj
                     continue
                 end
 
-                # Initialize currents
-                current = zeros(T, size(matrix, 1))
-                current[comp_i] = -1
-                current[comp_j] = 1
-
-                # Solve system
-                # csinfo("Solving points $pi and $pj")
-                log && csinfo("Solving pair $(d[(pi,pj)]) of $num")
-                t2 = @elapsed v = solve_linear_system(cfg, matrix, current, P)
-                csinfo("Time taken to solve linear system = $t2 seconds")
-                v .= v .- v[comp_i]
-
-                # Calculate resistance
-                r = v[comp_j] - v[comp_i]
+                # Forget excluded pairs
+                ex = false
+                for c_i in I
+                    for c_j in J
+                        if (c_i, c_j) in exclude
+                            continue
+                        end
 
-                # Return resistance value
-                for c_i in I, c_j in J
-                    push!(ret, (c_i, c_j, r))
-                    if get_shortcut_resistances
-                        resistances[c_i, c_j] = r
-                        resistances[c_j, c_i] = r
+                        # Initialize currents
+                        current = zeros(T, size(matrix, 1))
+                        current[comp_i] = -1
+                        current[comp_j] = 1
+
+                        # Solve system
+                        # csinfo("Solving points $pi and $pj")
+                        log && csinfo("Solving pair $(d[(pi,pj)]) of $num")
+                        t2 = @elapsed v = solve_linear_system(cfg, matrix, current, P)
+                        csinfo("Time taken to solve linear system = $t2 seconds")
+                        v .= v .- v[comp_i]
+
+                        # Calculate resistance
+                        r = v[comp_j] - v[comp_i]
+
+                        # Return resistance value
+                        push!(ret, (c_i, c_j, r))
+                        if get_shortcut_resistances
+                            resistances[c_i, c_j] = r
+                            resistances[c_j, c_i] = r
+                        end
+                        output = Output(points, v, (orig_pts[c_i], orig_pts[c_j]),
+                                        (comp_i, comp_j), r, V(c_j), cum)
+                        postprocess(output, component_data, flags, shortcut, cfg)
                     end
-                    output = Output(points, v, (orig_pts[c_i], orig_pts[c_j]),
-                                    (comp_i, comp_j), r, V(c_j), cum)
-                    postprocess(output, component_data, flags, shortcut, cfg)
                 end
             end
 
@@ -306,7 +302,8 @@ function _cholmod_solver_path(data::GraphData{T,V}, flags,
 
     get_shortcut_resistances = false
     if is_raster && !write_volt_maps && !write_cur_maps && 
-            !write_cum_cur_map_only  && !write_max_cur_maps
+            !write_cum_cur_map_only  && !write_max_cur_maps &&
+            isempty(exclude)
         get_shortcut_resistances = true
         csinfo("Triggering resistance calculation shortcut")
         num, d = get_num_pairs_shortcut(cc, points, exclude)
@@ -359,23 +356,18 @@ function _cholmod_solver_path(data::GraphData{T,V}, flags,
                 comp_j = V(something(findfirst(isequal(pj), comp),0))
                 J = findall(x -> x == pj, points)
 
-                # Forget excluded pairs
-                ex = false
-                for c_i in I, c_j in J
-                    if (c_i, c_j) in exclude
-                        ex = true
-                        break
-                    end
-                end
-                ex && continue
-
                 if pi == pj
                     continue
                 end
 
+                # Forget excluded pairs
                 for c_i in I, c_j in J
-                    push!(cholmod_batch, 
-                      CholmodNode((comp_i, comp_j), (V(c_i), V(c_j))))
+                    if (c_i, c_j) in exclude
+                        continue
+                    else
+                        push!(cholmod_batch, 
+                          CholmodNode((comp_i, comp_j), (V(c_i), V(c_j))))
+                    end
                 end
             end
         end

src/io.jl

@@ -121,7 +121,7 @@ function _ascii_grid_read_header(habitat_file, f)
     RasterMeta(ncols, nrows, xllcorner, yllcorner, cellsize, nodata, file_type)
 end
 
-function _guess_file_type(filename, f)
+#=function _guess_file_type(filename, f)
     s = readline(f)
     seek(f, 0)
 
@@ -137,6 +137,26 @@ function _guess_file_type(filename, f)
         throw("Check file format")
     end
 
+end=#
+
+function _guess_file_type(filename, f)
+
+    hdr = readline(f)
+    seek(f, 0)
+
+    if startswith(hdr, FILE_HDR_NPY)
+        filetype = FILE_TYPE_NPY
+    elseif startswith(lowercase(hdr), FILE_HDR_AAGRID)
+        filetype = FILE_TYPE_AAGRID
+    elseif startswith(hdr, FILE_HDR_INCL_PAIRS_AAGRID)
+        filetype = FILE_TYPE_INCL_PAIRS_AAGRID
+    elseif startswith(hdr, FILE_HDR_INCL_PAIRS)
+        filetype = FILE_TYPE_INCL_PAIRS
+    else
+        filetype = FILE_TYPE_TXTLIST
+    end
+
+    return filetype
 end
 
 function read_polymap(T, file::String, habitatmeta;
@@ -173,13 +193,13 @@ function read_point_map(V, file, habitatmeta)
 
     f = endswith(file, ".gz") ? GZip.open(file, "r") : open(file, "r")
     filetype = _guess_file_type(file, f)
-    _points_rc = filetype == TXTLIST ? readdlm(file) :
+    _points_rc = filetype == FILE_TYPE_TXTLIST ? readdlm(file) :
                         read_polymap(V, file, habitatmeta)
 
     i = V[]
     j = V[]
     v = V[]
-    if filetype == TXTLIST
+    if filetype == FILE_TYPE_TXTLIST
         I = _points_rc[:,2]
         J = _points_rc[:,3]
         v = _points_rc[:,1]
@@ -228,7 +248,7 @@ function read_source_and_ground_maps(T, V, source_file, ground_file, habitatmeta
     f = endswith(ground_file, "gz") ? Gzip.open(ground_file, "r") : open(ground_file, "r")
     filetype = _guess_file_type(ground_file, f)
 
-    if filetype == AAGRID
+    if filetype == FILE_TYPE_AAGRID
         ground_map = read_polymap(T, ground_file, habitatmeta; nodata_as = -1)
         ground_map = map(T, ground_map)
     else
@@ -240,7 +260,7 @@ function read_source_and_ground_maps(T, V, source_file, ground_file, habitatmeta
     f = endswith(source_file, "gz") ? Gzip.open(source_file, "r") : open(source_file, "r")
     filetype = _guess_file_type(source_file, f)
 
-    if filetype == AAGRID
+    if filetype == FILE_TYPE_AAGRID
         source_map = read_polymap(T, source_file, habitatmeta)
         source_map = map(T, source_map)
     else
@@ -261,36 +281,43 @@ function read_source_and_ground_maps(T, V, source_file, ground_file, habitatmeta
     source_map, ground_map
 end
 
-function read_included_pairs(V, file)
+function read_included_pairs(V, filename)
 
-    f = endswith(file, "gz") ? Gzip.open(file, "r") : open(file, "r")
-    filetype = _guess_file_type(file, f)
+    f = endswith(filename, "gz") ? Gzip.open(filename, "r") : open(filename, "r")
+    filetype = _guess_file_type(filename, f)
     minval = 0
     maxval = 0
     mode = :undef
 
-    if filetype == PAIRS_AAGRID
-        open(file, "r") do f
+    if filetype == FILE_TYPE_INCL_PAIRS_AAGRID
+        open(filename, "r") do fV
             minval = parse(Float64, split(readline(f))[2])
             maxval = parse(Float64, split(readline(f))[2])
         end
-        included_pairs = readdlm(file, skipstart=2)
+        included_pairs = readdlm(filename, skipstart=2)
         point_ids = V.(included_pairs[:,1])
         deleteat!(point_ids, 1)
         included_pairs = included_pairs[2:end, 2:end]
         map!(x -> x > maxval ? 0 : x, included_pairs, included_pairs)
         idx = findall(x -> x >= minval, included_pairs)
         mode = :include
         bin = map(x -> x >= minval ? V(1) : V(0), included_pairs)
-        IncludeExcludePairs(mode, point_ids, bin)
-    else
-        open(file, "r") do f
+        return IncludeExcludePairs(mode, point_ids, bin)
+    elseif filetype == FILE_TYPE_INCL_PAIRS
+        open(filename, "r") do f
             mode = Symbol(split(readline(f))[2])
         end
-        included_pairs = readdlm(file, skipstart = 1)
+        included_pairs = readdlm(filename, skipstart = 1)
+        if size(included_pairs, 1) == 1
+            pl = zeros(V, 1,2)
+            pl[1,:] = included_pairs
+            included_pairs = pl
+        end
         point_ids = V.(sort!(unique(included_pairs)))
-        if point_ids[1] == 0
-            deleteat!(point_ids, 1)
+        idx = findall(x -> x == 0 , point_ids)
+        if length(idx) > 0
+            deleteat!(point_ids, idx)
+            cswarn("Code to include pairs is activated, some entries did not match with focal node file. Some focal nodes may have been dropped")
         end
 
         mat = zeros(V, size(point_ids, 1), size(point_ids, 1))
@@ -303,8 +330,10 @@ function read_included_pairs(V, file)
                 mat[idx2,idx1] = 1
             end
         end
-        IncludeExcludePairs(mode, point_ids, mat)
+
+        return IncludeExcludePairs(mode, point_ids, mat)
     end
+
 end
 
 function get_network_data(T, V, cfg)::NetworkData{T,V}

src/logging.jl

@@ -3,7 +3,7 @@ const fmt = x -> Dates.format(x, "yyyy-mm-dd HH:MM:SS")
 const ui_interface = Ref{Function}((x,y) -> nothing)
 
 csinfo(msg) = (@info(string(fmt(Dates.now())) * " : " * msg); ui_interface[](msg, :info))
-cswarn(msg) = (@warn(string(fmt(Dates.now())) * " : " * msg); ui_interace[](msg, :warn))
+cswarn(msg) = (@warn(string(fmt(Dates.now())) * " : " * msg); ui_interface[](msg, :warn))
 
 function update_logging!(cfg)
 

src/raster/onetoall.jl

@@ -82,7 +82,7 @@ function onetoall_kernel(data::RasData{T,V}, flags, cfg)::Matrix{T} where {T,V}
         s = copy(z)
         n = points_unique[i]
         if use_included_pairs
-            for j = 1:num_points_to_solve
+            for j = 1:size(point_ids,1)
                 if i != j && included_pairs.include_pairs[i,j] == mode
                     exclude = point_ids[j]
                     map!(x -> x == exclude ? 0 : x, point_map, point_map)

src/raster/pairwise.jl

@@ -237,14 +237,14 @@ function generate_exclude_pairs(points_rc, included_pairs::IncludeExcludePairs{V
     mat = included_pairs.include_pairs
     mode = included_pairs.mode == :include ? 0 : 1    
 
-        prune_points!(points_rc, included_pairs.point_ids)
-        for j = 1:size(mat, 2)
-            for i = 1:size(mat, 1)
-                if mat[i,j] == mode
-                    push!(exclude_pairs_array, (i,j))
-                end
+    prune_points!(points_rc, included_pairs.point_ids)
+    for j = 1:size(mat, 2)
+        for i = 1:size(mat, 1)
+            if mat[i,j] == mode && mat[j,i] == mode
+                push!(exclude_pairs_array, (i,j))
             end
         end
+    end
 
     exclude_pairs_array
 end

src/run.jl

@@ -14,11 +14,12 @@ Inputs:
 function compute(path::String)
     cfg = parse_config(path)
     update_logging!(cfg)
-    #=if cfg["use_64bit_indexing"] in TRUELIST
-        global INT = Int64
-    end=#
     write_config(cfg)
     T = cfg["precision"] in SINGLE ? Float32 : Float64
+    if T == Float32 && cfg["solver"] in CHOLMOD
+        cswarn("CHOLMOD solver mode works only in double precision. Switching precision to double.")
+        T = Float64
+    end
     V = cfg["use_64bit_indexing"] in TRUELIST ? Int64 : Int32
     csinfo("Precision used: $(cfg["precision"])")
     is_parallel = cfg["parallelize"] in TRUELIST