Add reader for my own netCDF based format

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [weakdeps]

src/Data/Data.jl

@@ -1,9 +1,14 @@
 module Data
 using IMDP, SparseArrays
+using NCDatasets
 
 include("bmdp-tool.jl")
 export read_bmdp_tool_file, write_bmdp_tool_file
 
 include("prism.jl")
 export read_prism_file, write_prism_file
+
+include("imdp.jl")
+export read_imdp_jl_file, write_imdp_jl_file
+
 end

src/Data/bmdp-tool.jl

@@ -91,7 +91,7 @@ end
 function write_bmdp_tool_file(path, mdp::IntervalMarkovDecisionProcess, terminal_states)
     prob = transition_prob(mdp)
     l, g = lower(prob), gap(prob)
-    nsrc = num_src(prob)
+    num_columns = num_src(prob)
     sptr = IMDP.stateptr(mdp)
     act = actions(mdp)
 
@@ -109,7 +109,7 @@ function write_bmdp_tool_file(path, mdp::IntervalMarkovDecisionProcess, terminal
         end
 
         s = 1
-        for j in 1:nsrc
+        for j in 1:num_columns
             action = act[j]
 
             if sptr[s + 1] == j
@@ -131,3 +131,38 @@ function write_bmdp_tool_file(path, mdp::IntervalMarkovDecisionProcess, terminal
         end
     end
 end
+
+function write_bmdp_tool_file(path, mdp::IntervalMarkovChain, terminal_states)
+    prob = transition_prob(mdp)
+    l, g = lower(prob), gap(prob)
+    num_columns = num_src(prob)
+
+    number_states = num_states(mdp)
+    number_terminal = length(terminal_states)
+
+    open(path, "w") do io
+        println(io, number_states)
+        println(io, 1)  # number_actions
+        println(io, number_terminal)
+
+        for terminal_state in terminal_states
+            println(io, terminal_state - 1)
+        end
+
+        for j in 1:num_columns
+            src = j - 1
+
+            column_lower = view(l, :, j)
+            I, V = SparseArrays.findnz(column_lower)
+
+            for (i, v) in zip(I, V)
+                dest = i - 1
+                pl = v
+                pu = pl + g[i, j]
+
+                transition = "$src 0 $dest $pl $pu"
+                println(io, transition)
+            end
+        end
+    end
+end

src/Data/imdp.jl

@@ -0,0 +1,57 @@
+
+function read_imdp_jl_file(path)
+    dataset = Dataset(path)
+
+    n = dataset.attrib["num_states"]
+    initial_state = dataset.attrib["initial_state"]
+    model = dataset.attrib["model"]
+
+    @assert model ∈ ["imdp", "imc"]
+    @assert dataset.attrib["rows"] == "to"
+    @assert dataset.attrib["cols"] ∈ ["from", "from/action"]
+    @assert dataset.properties["format"] == "sparse_csc"
+
+    lower_colptr = convert.(Int32, dataset["lower_colptr"][:])
+    lower_rowval = convert.(Int32, dataset["lower_rowval"][:])
+    lower_nzval = dataset["lower_nzval"][:]
+    P̲ = SparseMatrixCSC(Int32(n + 1), Int32(n), lower_colptr, lower_rowval, lower_nzval)
+
+    upper_colptr = convert.(Int32, dataset["upper_colptr"][:])
+    upper_rowval = convert.(Int32, dataset["upper_rowval"][:])
+    upper_nzval = dataset["upper_nzval"][:]
+    P̅ = SparseMatrixCSC(Int32(n + 1), Int32(n), upper_colptr, upper_rowval, upper_nzval)
+
+    prob = MatrixIntervalProbabilities(; lower = P̲, upper = P̅)
+
+    if model == "imdp"
+        mdp_or_mc = read_imdp_jl_mdp(dataset, prob, initial_state)
+    elseif model == "imc"
+        mdp_or_mc = read_imdp_jl_mc(dataset, prob, initial_state)
+    end
+
+    # IMPORTANT! Otherwise the file cannot be reopened until the OS has released the file handle.
+    close(dataset)
+
+    return mdp_or_mc
+end
+
+function read_imdp_jl_mdp(dataset, prob, initial_state)
+    @assert dataset.attrib["cols"] == "from/action"
+
+    stateptr = convert.(Int32, dataset["stateptr"][:])
+    action_vals = dataset["action_vals"][:]
+
+    mdp = IntervalMarkovDecisionProcess(prob, stateptr, action_vals, Int32(initial_state))
+    return mdp
+end
+
+function read_imdp_jl_mc(dataset, prob, initial_state)
+    @assert dataset.attrib["cols"] == "from"
+
+    mc = IntervalMarkovChain(prob, Int32(initial_state))
+    return mc
+end
+
+function write_imdp_jl_file(path, mdp_or_mc)
+    # TODO: implement
+end

src/Data/prism.jl

@@ -1,13 +1,13 @@
 
-function write_prism_file(path_without_file_ending, mdp, terminal_states)
-    write_prism_states_file(path_without_file_ending, mdp)
-    write_prism_labels_file(path_without_file_ending, mdp, terminal_states)
-    write_prism_transitions_file(path_without_file_ending, mdp)
+function write_prism_file(path_without_file_ending, mdp_or_mc, terminal_states)
+    write_prism_states_file(path_without_file_ending, mdp_or_mc)
+    write_prism_labels_file(path_without_file_ending, mdp_or_mc, terminal_states)
+    write_prism_transitions_file(path_without_file_ending, mdp_or_mc)
     return write_prism_props_file(path_without_file_ending)
 end
 
-function write_prism_states_file(path_without_file_ending, mdp)
-    number_states = num_states(mdp)
+function write_prism_states_file(path_without_file_ending, mdp_or_mc)
+    number_states = num_states(mdp_or_mc)
 
     lines = Vector{String}(undef, 1 + number_states)
     lines[1] = "(s)"
@@ -20,8 +20,8 @@ function write_prism_states_file(path_without_file_ending, mdp)
     return write(path_without_file_ending * ".sta", join(lines, "\n"))
 end
 
-function write_prism_labels_file(path_without_file_ending, mdp, terminal_states)
-    istate = initial_state(mdp) - 1
+function write_prism_labels_file(path_without_file_ending, mdp_or_mc, terminal_states)
+    istate = initial_state(mdp_or_mc) - 1
 
     lines = Vector{String}(undef, 2 + length(terminal_states))
     lines[1] = "0=\"init\" 1=\"deadlock\" 2=\"goal\""
@@ -35,13 +35,13 @@ function write_prism_labels_file(path_without_file_ending, mdp, terminal_states)
     return write(path_without_file_ending * ".lab", join(lines, "\n"))
 end
 
-function write_prism_transitions_file(path_without_file_ending, mdp)
+function write_prism_transitions_file(path_without_file_ending, mdp::IntervalMarkovDecisionProcess)
     number_states = num_states(mdp)
 
     prob = transition_prob(mdp)
     l, g = lower(prob), gap(prob)
 
-    transition_lines = num_src(prob)
+    num_columns = num_src(prob)
     num_transitions = nnz(l)
 
     sptr = IMDP.stateptr(mdp)
@@ -53,7 +53,7 @@ function write_prism_transitions_file(path_without_file_ending, mdp)
 
         s = 1
         action_idx = 0
-        for j in 1:transition_lines
+        for j in 1:num_columns
             action = act[j]
 
             if sptr[s + 1] == j
@@ -79,6 +79,36 @@ function write_prism_transitions_file(path_without_file_ending, mdp)
     end
 end
 
+function write_prism_transitions_file(path_without_file_ending, mc::IntervalMarkovChain)
+    number_states = num_states(mc)
+
+    prob = transition_prob(mc)
+    l, g = lower(prob), gap(prob)
+
+    num_columns = num_src(prob)
+    num_transitions = nnz(l)
+
+    open(path_without_file_ending * ".tra", "w") do io
+        println(io, "$number_states 1 $num_transitions")  # number_states number_choices number_transitions
+
+        for j in 1:num_columns
+            src = j - 1
+
+            column_lower = view(l, :, j)
+            I, V = SparseArrays.findnz(column_lower)
+
+            for (i, v) in zip(I, V)
+                dest = i - 1
+                pl = v
+                pu = pl + g[i, j]
+                pl = max(pl, 1e-12)
+
+                println(io, "$src 0 $dest [$pl,$pu] mc")
+            end
+        end
+    end
+end
+
 function write_prism_props_file(path_without_file_ending)
     line = "Pmaxmin=? [ F \"goal\" ]"