diff --git a/docs/_config.yml b/docs/_config.yml
new file mode 100644
index 0000000..fff4ab9
--- /dev/null
+++ b/docs/_config.yml
@@ -0,0 +1 @@
+theme: jekyll-theme-minimal
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..0394e55
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,38 @@
+using Documenter
+#include("../src/QSWalk.jl")
+using QSWalk
+
+# same for contributing and license
+cp(normpath(@__FILE__, "../../LICENSE"), normpath(@__FILE__, "../src/license.md"); force=true)
+
+
+
+makedocs(
+    modules     = [QSWalk],
+    format      = :html,
+    sitename    = "QSWalk",
+    clean       = true,
+    doctest     = true,
+    checkdocs   = :exports,
+    assets 	= ["assets/logo.ico"],
+    pages       = Any[
+		"Home"				=> "index.md",
+        "GKSL master equation" => "gksl.md",
+    	"Demoralization"   => "demoralization.md",
+        "Contributing"		=> "contributing.md",
+        "Citing"	       	=> "citing.md",
+		"Licence"			=> "license.md",
+    ]
+)
+
+deploydocs(
+    deps        = nothing,
+    make        = nothing,
+    repo        = "github.com/QuantumWalks/QSWalk.jl",
+    target      = "build",
+    julia       = "0.6",
+    osname      = "linux"
+)
+
+
+#rm(normpath(@__FILE__, "src/license.md"))
diff --git a/docs/src/assets/logo.ico b/docs/src/assets/logo.ico
new file mode 100644
index 0000000..533f0ed
Binary files /dev/null and b/docs/src/assets/logo.ico differ
diff --git a/docs/src/assets/logo.png b/docs/src/assets/logo.png
new file mode 100644
index 0000000..75936f1
Binary files /dev/null and b/docs/src/assets/logo.png differ
diff --git a/docs/src/citing.md b/docs/src/citing.md
new file mode 100644
index 0000000..f5aa9ce
--- /dev/null
+++ b/docs/src/citing.md
@@ -0,0 +1,19 @@
+## Usage and citing
+
+
+In case of citing, please use the following BibTeX form:
+
+```tex
+@article{glos2018qswalk,
+  title={QSWalk. jl: Julia package for quantum stochastic walks analysis},
+  author={Glos, Adam and Miszczak, Jaros{\l}aw Adam and Ostaszewski, Mateusz},
+  journal={arXiv preprint arXiv:1801.01294},
+  year={2018}
+}
+```
+
+Our package was already used in papers concerning quantum attacks
+* Adam Glos, Jarosław Adam Miszczak, and Mateusz Ostaszewski. 'Limiting properties of stochastic quantum walks on directed graphs.' Journal of Physics A: Mathematical and Theoretical 51.3 (2017): 035304.
+* Krzysztof Domino, Adam Glos, and Mateusz Ostaszewski. 'Superdiffusive quantum stochastic walk definable on arbitrary directed graph'. Quantum Information & Computation, 17(11-12), 973-986.
+
+In case You have used our package for your research, we will be grateful for any information about the paper or the package. With your consent we will provide a link to the paper here.
diff --git a/docs/src/contributing.md b/docs/src/contributing.md
new file mode 100644
index 0000000..e452f5c
--- /dev/null
+++ b/docs/src/contributing.md
@@ -0,0 +1,22 @@
+## Bugs
+
+In the case, you noticed some bugs, please start with an issue with a minimal
+working example of not-working code. If *Exception* is thrown, please provide an
+exception message as well.  If no *Exception* is thrown, but the result is
+wrong, please provide in the issue message correct answer.
+
+In the case you make a pull request, please add a not-working example as a test.
+
+## Improvements
+
+If you can provide a code, which works faster than already existing, please
+check its efficiency for various input data. 
+
+We welcome any ideas concerning the readability and logic of the code as well.
+
+## Development guidelines
+* Post an issue.
+* Wait until the discussion ends.
+* Create assertions on argument types and other requirements.
+* Include necessary references.
+* Write tests.
diff --git a/docs/src/demoralization.md b/docs/src/demoralization.md
new file mode 100644
index 0000000..bb490e9
--- /dev/null
+++ b/docs/src/demoralization.md
@@ -0,0 +1,39 @@
+```@meta
+DocTestSetup = quote
+   using QSWalk
+end
+```
+
+## Nonmoralizing model
+
+Global interaction quantum stochastic walk suffers for creating additional connections. For
+removing such effect, nonmoralizing quantum stochastic walk was introduced, see [here](http://www.rintonpress.com/journals/qiconline.html#v17n1112).
+Such model is constructed in several steps. First, the dimensionality is increased, hence to each
+vertex multidimensional subspaces is attached. Then, Hamiltonian and Lindblad
+operator is increased, furthermore additional Hamiltonian called "local Hamiltonian".
+is introduced.
+
+Below we present additional functionalities typical for nonmoralizing quantum
+stochastic walk. By default the the operator is generalized as in the original [paper](http://www.rintonpress.com/journals/qiconline.html#v17n1112).
+
+```@index
+Order = [:type, :function]
+Modules = [QSWalk]
+Pages   = ["demoralization.md"]
+```
+
+## Full docs
+
+```@docs
+Vertex
+VertexSet
+nm_measurement
+nm_loc_ham
+nm_init
+default_nm_loc_ham
+make_vertex_set
+vlist
+subspace
+fourier_matrix
+vertexsetsize
+```
diff --git a/docs/src/gksl.md b/docs/src/gksl.md
new file mode 100644
index 0000000..6da6adb
--- /dev/null
+++ b/docs/src/gksl.md
@@ -0,0 +1,32 @@
+```@meta
+DocTestSetup = quote
+   using QSWalk
+end
+```
+
+## GKSL master equation
+
+GKSL master equation is general continuous-time open quantum evolution. The master equation base on to form of operators: the Hamiltonian, which describes the evolution of closed system, and the Lindbladian operators which describes the evolution of open system. Basic facts in context of quantum stochastic walks can found [here](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.81.022323). For local interaction, where each Lindblad operator is a matrix with single nonzero element, we created a `local_lind` function which splits matrix into mentioned operators.
+
+Below we present a documentation of basic function used for simulating GKSL master equation.
+
+```@index
+Order = [:type, :function]
+Modules = [QSWalk]
+Pages   = ["gksl.md"]
+```
+
+## Full docs
+
+```@docs
+ket
+bra
+ketbra
+proj
+res
+unres
+evolve_generator(::AbstractMatrix{<:Number}, ::AbstractVector{<:AbstractMatrix{<:Number}}, ::AbstractMatrix{<:Number}, ::Real)
+evolve_operator
+evolve(::AbstractMatrix{<:Number}, ::AbstractMatrix{<:Number})
+local_lind
+```
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..6f6fc64
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,9 @@
+## Welcome to QSWalk.jl
+
+*QSWalk.jl* is a package providing implementation of open continuous-time quantum walk evolution based on GKSL master equation. In particular we have implemented function useful for analysing local interaction, global interaction and nonmoralizing global interaction stochastic quantum walk models.
+
+In GitHub we present an [examples](https://github.com/QuantumWalks/QSWalk.jl/) in .ipynb and .jl extension. The present most of the functionalities of the package. Furthermore, for package explanation we refer the user to our [paper](https://arxiv.org/abs/1801.01294) describing the package, and
+for the basic papers describing and introducing quantum stochastic models:
+* [Quantum stochastic walks: A generalization of classical random walks and quantum walks](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.81.022323) by Whitfield, Rodríguez-Rosario, and Aspuru-Guzik,
+* [Superdiffusive quantum stochastic walk definable on arbitrary directed graph](http://www.rintonpress.com/journals/qiconline.html#v17n1112) by Domino, Glos, and Ostaszewski,
+* [Properties of quantum stochastic walks from the asymptotic scaling exponent](http://www.rintonpress.com/journals/qiconline.html#v18n34) by Domino, Glos, Ostaszewski, Pawela, and Sadowski.
diff --git a/docs/src/license.md b/docs/src/license.md
new file mode 100644
index 0000000..1849e3d
--- /dev/null
+++ b/docs/src/license.md
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2017-2018 Adam Glos, Jarosław Adam Miszczak, Mateusz Ostaszewski
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/examples/ex01_path_propagation.jl b/examples/ex01_path_propagation.jl
index 49822a2..1d5eac9 100644
--- a/examples/ex01_path_propagation.jl
+++ b/examples/ex01_path_propagation.jl
@@ -8,7 +8,6 @@
 using QSWalk
 using PyPlot # for plot
 using LightGraphs # for PathGraph
-using LinearAlgebra
 
 
 # ## Evolution on path graph for fixed time, global vs local interaction
@@ -38,7 +37,7 @@ rho_local = evolve(op_local, proj(midpoint, dim), timepoint);
 
 # To plot the result of the cannonical measurement, we take a diagonal of the states. Note that both *rhoglobal* and *rholocal* are complex matrices, hence we need to take the real part only. Argument *positions* is generated assuming that *midpoint* (the initial position) is at 0. Note that we have very heavy tails in the global interaction case, which confirms fast propagation in this model.
 
-positions = (collect(1:dim) .- midpoint)
+positions = (collect(1:dim) - midpoint)
 plot(positions, real.(diag(rho_local)), "k")
 plot(positions, real.(diag(rho_global)), "b")
 xlabel("position")
@@ -75,7 +74,7 @@ local_states = evolve(op_local, proj(midpoint, dim), timepoints);
 
 secmoment_global = Float64[]
 secmoment_local = Float64[]
-positions = (collect(1:dim) .- midpoint)
+positions = (collect(1:dim) - midpoint)
 
 for (rho_global, rho_local) = zip(global_states, local_states)
   push!(secmoment_global, sum(positions.^2 .* diag(rho_global)))
@@ -90,4 +89,3 @@ plot(timepoints, secmoment_global, "b")
 xlabel("t")
 ylabel(L"\mu_2")
 show()
-
diff --git a/src/demoralization.jl b/src/demoralization.jl
index 8da1404..bcd260e 100644
--- a/src/demoralization.jl
+++ b/src/demoralization.jl
@@ -22,18 +22,18 @@ upper diagonal (equal to `1im`) and lower diagonal (equal to `-1im`).
 
 ```jldoctest
 julia> QSWalk.default_nm_loc_ham(4)
-4×4 SparseMatrixCSC{Complex{Float64}, Int64} with 6 stored entries:
-  [2, 1] = 0.0-1.0im
-  [1, 2] = 0.0+1.0im
-  [3, 2] = 0.0-1.0im
-  [2, 3] = 0.0+1.0im
-  [4, 3] = 0.0-1.0im
-  [3, 4] = 0.0+1.0im
+4×4 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 6 stored entries:
+  [2, 1]  =  0.0-1.0im
+  [1, 2]  =  0.0+1.0im
+  [3, 2]  =  0.0-1.0im
+  [2, 3]  =  0.0+1.0im
+  [4, 3]  =  0.0-1.0im
+  [3, 4]  =  0.0+1.0im
 ```
 """
 function default_nm_loc_ham(size::Int)
   @argumentcheck size>0 "Size of default local Hamiltonian needs to be positive"
-  if size ==  1
+  if size == 1
     return spzeros(ComplexF64, 1, 1)
   else
     spdiagm(1=>im*ones(size-1), -1=>-im*ones(size-1))
@@ -43,15 +43,12 @@ end
 """
 
     nm_loc_ham(vertexset[, hamiltoniansByDegree])
-    nm_loc_ham(vertexset, hamiltoniansByVertex)
 
 Return Hamiltonian acting locally on each vertex from `vertexset` linear
-subspace. In the first form, `hamiltoniansByDegree` is a dictionary `Dict{Int,
+subspace. `hamiltoniansByDegree` is a dictionary `Dict{Int,
 SparseDenseMatrix}`, which, for a given dimension of vertex linear subspace,
 yields a hermitian operator. Only matrices for existing dimensions needs to be
-specified. In the second form, `hamiltoniansByVertex` is a dictionary
-`Dict{Vertex, SparseDenseMatrix}`, which, for a given vertex, yields a hermitian
-operator of the size equal to the dimension of the vertex subspace.
+specified.
 
 *Note:* Value of `vertexset` should be generated by `make_vertex_set` in order
 to match demoralization procedure. Numerical analysis suggests, that
@@ -61,50 +58,36 @@ hamiltonians should be complex valued.
 
 ```jldoctest
 julia> vset = VertexSet([[1, 2], [3, 4]])
-QSWalk.VertexSet(QSWalk.Vertex[QSWalk.Vertex([1, 2]), QSWalk.Vertex([3, 4])])
+VertexSet(Vertex[Vertex([1, 2]), Vertex([3, 4])])
 
-julia> full(nm_loc_ham(vset))
-4×4 Array{Complex{Float64}, 2}:
+julia> Matrix(nm_loc_ham(vset))
+4×4 Array{Complex{Float64},2}:
  0.0+0.0im  0.0+1.0im  0.0+0.0im  0.0+0.0im
  0.0-1.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+1.0im
  0.0+0.0im  0.0+0.0im  0.0-1.0im  0.0+0.0im
 
-julia> A, B = rand(2, 2), rand(2, 2)
-([0.358914 0.183322; 0.379927 0.671986], [0.26643 0.969279; 0.313752 0.636789])
-
-julia> v1, v2 = vlist(vset)
-2-element Array{QSWalk.Vertex, 1}:
- QSWalk.Vertex([1, 2])
- QSWalk.Vertex([3, 4])
-
-julia> nm_loc_ham(vset, Dict(v1  => A, v2  => B))
-4×4 SparseMatrixCSC{Complex{Float64}, Int64} with 8 stored entries:
-  [1, 1] = 0.358914+0.0im
-  [2, 1] = 0.379927+0.0im
-  [1, 2] = 0.183322+0.0im
-  [2, 2] = 0.671986+0.0im
-  [3, 3] = 0.26643+0.0im
-  [4, 3] = 0.313752+0.0im
-  [3, 4] = 0.969279+0.0im
-  [4, 4] = 0.636789+0.0im
-
-julia> nm_loc_ham(VertexSet([[1, 2], [3, 4]]), Dict(2  => A))
-4×4 SparseMatrixCSC{Complex{Float64}, Int64} with 8 stored entries:
-  [1, 1] = 0.358914+0.0im
-  [2, 1] = 0.379927+0.0im
-  [1, 2] = 0.183322+0.0im
-  [2, 2] = 0.671986+0.0im
-  [3, 3] = 0.358914+0.0im
-  [4, 3] = 0.379927+0.0im
-  [3, 4] = 0.183322+0.0im
-  [4, 4] = 0.671986+0.0im
+julia> A = [1 1im; -1im 1]
+2×2 Array{Complex{Int64},2}:
+ 1+0im  0+1im
+ 0-1im  1+0im
+
+julia> nm_loc_ham(vset, Dict{Int,Matrix{ComplexF64}}(2  => A))
+4×4 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 8 stored entries:
+  [1, 1]  =  1.0+0.0im
+  [2, 1]  =  0.0-1.0im
+  [1, 2]  =  0.0+1.0im
+  [2, 2]  =  1.0+0.0im
+  [3, 3]  =  1.0+0.0im
+  [4, 3]  =  0.0-1.0im
+  [3, 4]  =  0.0+1.0im
+  [4, 4]  =  1.0+0.0im
 ```
 """
 function nm_loc_ham(vset::VertexSet,
                     hamiltonians::Dict{Int, <:AbstractMatrix{<:Number}}
                         = Dict(l =>default_nm_loc_ham(l) for l = length.(vlist(vset))))
-  @argumentcheck all([size(h, 1) ==  size(h, 2) for h = values(hamiltonians)])  "Hamiltonians must consists of square matrices"
+  @argumentcheck all([size(h, 1) == size(h, 2) for h = values(hamiltonians)])  "Hamiltonians must consists of square matrices"
   verticeslengths = length.(vlist(vset))
   @assert all([l in keys(hamiltonians) for l = verticeslengths]) "Missing degree in the Dictionary: $verticeslengths needed"
 
@@ -112,9 +95,53 @@ function nm_loc_ham(vset::VertexSet,
   nm_loc_ham(vset, hamiltonianlist)
 end
 
+"""
+
+    nm_loc_ham(vertexset[, hamiltoniansByVertex])
+
+Return Hamiltonian acting locally on each vertex from `vertexset` linear
+subspace. `hamiltoniansByVertex` is a dictionary
+`Dict{Vertex, SparseDenseMatrix}`, which, for a given vertex, yields a hermitian
+operator of the size equal to the dimension of the vertex subspace.
+
+*Note:* Value of `vertexset` should be generated by `make_vertex_set` in order
+to match demoralization procedure. Numerical analysis suggests, that
+hamiltonians should be complex valued.
+
+# Examples
+
+```jldoctest
+julia> vset = VertexSet([[1, 2], [3, 4]])
+VertexSet(Vertex[Vertex([1, 2]), Vertex([3, 4])])
+
+julia> Matrix(nm_loc_ham(vset))
+4×4 Array{Complex{Float64},2}:
+ 0.0+0.0im  0.0+1.0im  0.0+0.0im  0.0+0.0im
+ 0.0-1.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
+ 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+1.0im
+ 0.0+0.0im  0.0+0.0im  0.0-1.0im  0.0+0.0im
+
+julia> A, B = [1 1im; -1im 1], [0 1; 1 0]
+(Complex{Int64}[1+0im 0+1im; 0-1im 1+0im], [0 1; 1 0])
+
+julia> v1, v2 = vlist(vset)
+2-element Array{Vertex,1}:
+ Vertex([1, 2])
+ Vertex([3, 4])
+
+julia> nm_loc_ham(vset, Dict{Vertex,Matrix{Number}}(v1  => A, v2  => B))
+4×4 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 6 stored entries:
+  [1, 1]  =  1.0+0.0im
+  [2, 1]  =  0.0-1.0im
+  [1, 2]  =  0.0+1.0im
+  [2, 2]  =  1.0+0.0im
+  [4, 3]  =  1.0+0.0im
+  [3, 4]  =  1.0+0.0im
+```
+"""
 function nm_loc_ham(vset::VertexSet,
                     hamiltonians::Dict{Vertex,<:AbstractMatrix{<:Number}})
-  @argumentcheck all([size(h, 1) ==  size(h, 2) for h = values(hamiltonians)])  "Hamiltonians must consists of square matrices"
+  @argumentcheck all([size(h, 1) == size(h, 2) for h = values(hamiltonians)])  "Hamiltonians must consists of square matrices"
   @assert all([v in keys(hamiltonians) for v = vlist(vset)]) "Missing hamiltonian for some vertex"
   @assert all([length(v) == size(hamiltonians[v], 1) for v = vlist(vset)]) "The vertex length and hamiltonian size do no match"
 
@@ -218,7 +245,7 @@ function nm_lind(A::AbstractMatrix{<:Number},
       L[subspace(vset[i]), k] = A[i, j]*lindbladians[length(vset[i])][:, index]
   end
   L, vset
-end
+end,
 
 function nm_lind(A::AbstractMatrix{<:Number};
                  epsilon::Real = eps())
@@ -226,12 +253,12 @@ function nm_lind(A::AbstractMatrix{<:Number};
   degrees = [length(v) for v = vlist(vset)]
 
   nm_lind(A, Dict(d =>fourier_matrix(d) for d = degrees), epsilon = epsilon)
-end
+end,
 
 function nm_lind(A::AbstractMatrix{<:Number},
                  lindbladians::Dict{Vertex, <:AbstractMatrix{<:Number}};
                  epsilon::Real = eps())
-  @argumentcheck all([size(l, 1) ==  size(l, 2) for l = values(lindbladians)])  "lindbladians should consist of square matrix"
+  @argumentcheck all([size(l, 1) == size(l, 2) for l = values(lindbladians)])  "lindbladians should consist of square matrix"
   @argumentcheck epsilon>= 0 "epsilon needs to be nonnegative"
 
   revincidence_list = reversed_incidence_list(A; epsilon = epsilon)
@@ -314,12 +341,12 @@ function nm_glob_ham(A::AbstractMatrix{<:Number},
     if index < j
       hamiltonianshape = length.((subspace(vset[index]), subspace(vset[j])))
       @argumentcheck hamiltonianshape in keys(hamiltonians) "hamiltonian of size $hamiltonianshape not found"
-      @argumentcheck hamiltonianshape ==  size(hamiltonians[hamiltonianshape]) "hamiltonian for key $(hamiltonianshape) shoud have shape $(hamiltonianshape)"
+      @argumentcheck hamiltonianshape == size(hamiltonians[hamiltonianshape]) "hamiltonian for key $(hamiltonianshape) shoud have shape $(hamiltonianshape)"
       H[subspace(vset[index]), subspace(vset[j])] = A[index, j]*hamiltonians[hamiltonianshape]
     end
   end
   H + H'
-end
+end,
 
 function nm_glob_ham(A::T;
                      epsilon::Real = eps()) where T<:AbstractMatrix{<:Number}
@@ -333,7 +360,7 @@ function nm_glob_ham(A::T;
     alloneshamiltonians[length.((w, v))] = ones(length(w), length(v))
   end
   nm_glob_ham(A, alloneshamiltonians, epsilon = epsilon)
-end
+end,
 
 function nm_glob_ham(A::AbstractMatrix{<:Number},
                      hamiltonians::Dict{Tuple{Vertex, Vertex}, <:AbstractMatrix{<:Number}};
@@ -348,7 +375,7 @@ function nm_glob_ham(A::AbstractMatrix{<:Number},
     if index < j
       key = (vset[index], vset[j])
       @argumentcheck key in keys(hamiltonians) "hamiltonian for $key not found"
-      @argumentcheck length.(key) ==  size(hamiltonians[key]) "hamiltonian for key $key shoud have shape $(length.(key))"
+      @argumentcheck length.(key) == size(hamiltonians[key]) "hamiltonian for key $key shoud have shape $(length.(key))"
       H[subspace(vset[j]), subspace(vset[index])] = A[index, j]*transpose(hamiltonians[key])
     end
   end
@@ -362,18 +389,13 @@ end
 Return joint probability of `probability`, which is real-valued probability vector
 according to `vertexset`.
 
-    nm_measurement(state, vertexset)
-
-Return joint probability of cannonical measurement of density matrix `state`,
-according to `vertexset`.
-
-*Note:* It is up to user to provide proper probability vector/density state.
+*Note:* It is up to user to provide proper probability vector.
 
 # Examples
 
 ```jldoctest
 julia> probability = [0.05, 0.1, 0.25, 0.3, 0.01, 0.20, 0.04, 0.05]
-8-element Array{Float64, 1}:
+8-element Array{Float64,1}:
  0.05
  0.1
  0.25
@@ -384,54 +406,59 @@ julia> probability = [0.05, 0.1, 0.25, 0.3, 0.01, 0.20, 0.04, 0.05]
  0.05
 
 julia> nm_measurement(probability, VertexSet([[1, 4], [2, 3, 5], [6], [7, 8]]))
-4-element Array{Float64, 1}:
+4-element Array{Float64,1}:
  0.35
  0.36
  0.2
  0.09
-
- julia> state = [1/6 0 1/6; 0 2/3 0; 1/6 0 1/6]
- 3×3 Array{Float64, 2}:
-  0.166667  0.0       0.166667
-  0.0       0.666667  0.0
-  0.166667  0.0       0.166667
-
- julia> nm_measurement(state, VertexSet([[1, 3], [2]]))
- 2-element Array{Float64, 1}:
-  0.333333
-  0.666667
 ```
 """
 function nm_measurement(probability::AbstractVector{<:Number},
                         vset::VertexSet)
-  @assert vertexsetsize(vset) ==  length(probability) "vertexset size and probability vector length do not match"
+  @assert vertexsetsize(vset) == length(probability) "vertexset size and probability vector length do not match"
 
   [sum(probability[subspace(vertex)]) for vertex = vlist(vset)]
 end
 
+"""
+
+    nm_measurement(state, vertexset)
+
+Return joint probability of cannonical measurement of density matrix `state`,
+according to `vertexset`.
+
+*Note:* It is up to user to provide proper density state.
+
+# Examples
+
+```jldoctest
+julia> state = [1/6 0 1/6; 0 2/3 0; 1/6 0 1/6]
+3×3 Array{Float64,2}:
+ 0.166667  0.0       0.166667
+ 0.0       0.666667  0.0
+ 0.166667  0.0       0.166667
+
+julia> nm_measurement(state, VertexSet([[1, 3], [2]]))
+2-element Array{Float64,1}:
+ 0.3333333333333333
+ 0.6666666666666666
+```
+"""
 function nm_measurement(state::AbstractMatrix{<:Number},
                         vset::VertexSet)
-  @argumentcheck size(state, 1) ==  size(state, 2) "state should be square matrix"
-  @assert vertexsetsize(vset) ==  size(state, 1) "vertexset size and state size do not match"
+  @argumentcheck size(state, 1) == size(state, 2) "state should be square matrix"
+  @assert vertexsetsize(vset) == size(state, 1) "vertexset size and state size do not match"
 
   nm_measurement(real.(diag(state)), vset)
 end
 
 """
-
     nm_init(init_vertices, vertexset)
-    nm_init(init_states, vertexset)
 
-Create initial state in the case of the nonmoralizing evolution. The result is
+Create initial state in the case of the nonmoralizing evolution based on
+`init_vertices` of type `Vector{Vertex}`. The result is
 a block diagonal matrix, where each block corresponds to vertex from `vertexset`.
-If the first argument is of type `Vector{Vertex}`, then default block matrix `eye`.
-
-If first argument is of type `Dict{Vertex, SparseDenseMatrix}`,
-then for each given vertex a block from dictionary is used, otherwise zero matrix
-is chosen. Each matrix from dictionary should be nonnegative and sum of all traces
-should equal one. The keys of `init_vertices` should be a subset of `vertexset()`.
-Note that matrix from `init_states` corresponding to vertex `v` should be of
-size `length(v)`×`length(v)`.
+The final state represent an uniform probability over `nm_measurement`.
 
 *Note:* The function returns sparse matrix with `ComplexF64` field type.
 
@@ -439,37 +466,16 @@ size `length(v)`×`length(v)`.
 
 ```jldoctest
 julia> vset = VertexSet([[1], [2, 3, 4], [5, 6, 7], [8, 9]])
-QSWalk.VertexSet(QSWalk.Vertex[QSWalk.Vertex([1]), QSWalk.Vertex([2, 3, 4]), QSWalk.Vertex([5, 6, 7]), QSWalk.Vertex([8, 9])])
+VertexSet(Vertex[Vertex([1]), Vertex([2, 3, 4]), Vertex([5, 6, 7]), Vertex([8, 9])])
 
 julia> nm_init(vset[[1, 3, 4]], vset)
-9×9 sparse matrix with 6 Complex{Float64} nonzero entries:
-	[1, 1] = 0.333333+0.0im
-	[5, 5] = 0.111111+0.0im
-	[6, 6] = 0.111111+0.0im
-	[7, 7] = 0.111111+0.0im
-	[8, 8] = 0.166667+0.0im
-	[9, 9] = 0.166667+0.0im
-
-julia> A1, A2, A3 = ones(ComplexF64, 1, 1)/4, [ 1/5+0im 0 1/5; 0 1/10 0 ; 1/5 0 1/5 ], [0.125 -0.125+0im; -0.125 0.125]
-(
-Complex{Float64}[0.25+0.0im],
-
-Complex{Float64}[0.2+0.0im 0.0+0.0im 0.2+0.0im; 0.0+0.0im 0.1+0.0im 0.0+0.0im; 0.2+0.0im 0.0+0.0im 0.2+0.0im],
-
-Complex{Float64}[0.125+0.0im -0.125+0.0im; -0.125+0.0im 0.125+0.0im])
-
-julia> nm_init(Dict(vset[1] =>A1, vset[3] =>A2, vset[4] =>A3), vset)
-9×9 sparse matrix with 10 Complex{Float64} nonzero entries:
-	[1, 1] = 0.25+0.0im
-	[5, 5] = 0.2+0.0im
-	[7, 5] = 0.2+0.0im
-	[6, 6] = 0.1+0.0im
-	[5, 7] = 0.2+0.0im
-	[7, 7] = 0.2+0.0im
-	[8, 8] = 0.125+0.0im
-	[9, 8] = -0.125+0.0im
-	[8, 9] = -0.125+0.0im
-	[9, 9] = 0.125+0.0im
+9×9 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 6 stored entries:
+  [1, 1]  =  0.333333+0.0im
+  [5, 5]  =  0.111111+0.0im
+  [6, 6]  =  0.111111+0.0im
+  [7, 7]  =  0.111111+0.0im
+  [8, 8]  =  0.166667+0.0im
+  [9, 9]  =  0.166667+0.0im
 ```
 """
 function nm_init(initialvertices::AbstractVector{Vertex},
@@ -484,6 +490,43 @@ function nm_init(initialvertices::AbstractVector{Vertex},
   L
 end
 
+"""
+    nm_init(init_states, vertexset)
+
+Create initial state in the case of the nonmoralizing evolution based on `init_states`
+of type `Dict{Vertex, <:AbstractMatrix{<:Number}}`. For each given vertex a block
+from dictionary is used, otherwise zero matrix is chosen. Each matrix from
+dictionary should be nonnegative and sum of all traces should equal one. The
+keys of `init_vertices` should be a vertices from `vertexset`.
+Note that matrix from `init_states` corresponding to vertex `v` should be of
+size `length(v)`×`length(v)`.
+
+*Note:* The function returns sparse matrix with `ComplexF64` field type.
+
+# Examples
+
+```jldoctest
+julia> vset = VertexSet([[1], [2, 3, 4], [5, 6, 7], [8, 9]])
+VertexSet(Vertex[Vertex([1]), Vertex([2, 3, 4]), Vertex([5, 6, 7]), Vertex([8, 9])])
+
+julia> A1, A2, A3 = ones(ComplexF64, 1, 1)/4, [ 1/5+0im 0 1/5; 0 1/10 0 ; 1/5 0 1/5 ], [0.125 -0.125+0im; -0.125 0.125]
+(Complex{Float64}[0.25+0.0im], Complex{Float64}[0.2+0.0im 0.0+0.0im 0.2+0.0im; 0.0+0.0im 0.1+0.0im 0.0+0.0im; 0.2+0.0im 0.0+0.0im 0.2+0.0im], Complex{Float64}[0.125+0.0im -0.125+0.0im; -0.125+0.0im 0.125+0.0im])
+
+julia> nm_init(Dict(vset[1] =>A1, vset[3] =>A2, vset[4] =>A3), vset)
+9×9 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 10 stored entries:
+  [1, 1]  =  0.25+0.0im
+  [5, 5]  =  0.2+0.0im
+  [7, 5]  =  0.2+0.0im
+  [6, 6]  =  0.1+0.0im
+  [5, 7]  =  0.2+0.0im
+  [7, 7]  =  0.2+0.0im
+  [8, 8]  =  0.125+0.0im
+  [9, 8]  =  -0.125+0.0im
+  [8, 9]  =  -0.125+0.0im
+  [9, 9]  =  0.125+0.0im
+
+```
+"""
 function nm_init(initial_states::Dict{Vertex, <:AbstractMatrix{<:Number}},
                  vset::VertexSet)
   @assert all([size(state, 1) == length(k) for (k, state) = initial_states]) "The size of initial state and the vertex do not match"
diff --git a/src/dirac.jl b/src/dirac.jl
index d765fc9..849c312 100644
--- a/src/dirac.jl
+++ b/src/dirac.jl
@@ -18,9 +18,8 @@ To be consistent with Julia indexing, `index` = 1, 2, ..., `size`.
 
 ```jldoctest
 julia> ket(1, 2)
-2-element SparseVector{Int64, Int64} with 1 stored entry:
-  [1] = 1
-
+2-element SparseArrays.SparseVector{Int64,Int64} with 1 stored entry:
+  [1]  =  1
 ```
 """
 function ket(index::Int, size::Int)
@@ -36,14 +35,14 @@ end
 """
     bra(index, size)
 
-Return `index`-th row vector in the `size`-dimensional vector space, with
+Return `index`-th base row vector in the `size`-dimensional vector space, with
 `index` = 1, 2, ..., `size`.
 
 # Examples
 
 ```jldoctest
 julia> bra(1, 2)
-1×2 RowVector{Int64, SparseVector{Int64, Int64}}:
+1×2 LinearAlgebra.Adjoint{Int64,SparseArrays.SparseVector{Int64,Int64}}:
  1  0
 ```
 """
@@ -62,8 +61,8 @@ equal to one, located at position (`irow`, `icol`).
 
 ```jldoctest
 julia> ketbra(1, 2, 3)
-3×3 SparseMatrixCSC{Int64, Int64} with 1 stored entry:
-  [1, 2] = 1
+3×3 SparseArrays.SparseMatrixCSC{Int64,Int64} with 1 stored entry:
+  [1, 2]  =  1
 
 ```
 """
@@ -78,33 +77,36 @@ Return projector onto `index`-th base vector in `size`-dimensional vector space,
 with `index` = 1, 2, ..., `size`. This is equivalent to `ketbra(index, index,
 size)`.
 
+```jldoctest
+julia> proj(1, 2)
+2×2 SparseArrays.SparseMatrixCSC{Int64,Int64} with 1 stored entry:
+  [1, 1]  =  1
+```
+"""
+function proj(index::Int, size::Int)
+  ketbra(index, index, size)
+end
+
+"""
     proj(vector)
 
 Return projector onto the subspace spanned by vector `vector`.
 
 # Examples
 ```jldoctest
-julia> proj(1, 2)
-2×2 SparseMatrixCSC{Int64, Int64} with 1 stored entry:
-  [1, 1] = 1
-
 julia> v = 1/sqrt(2) * (ket(1, 3)+ket(3, 3))
-3-element SparseVector{Float64, Int64} with 2 stored entries:
-  [1] = 0.707107
-  [3] = 0.707107
+3-element SparseArrays.SparseVector{Float64,Int64} with 2 stored entries:
+  [1]  =  0.707107
+  [3]  =  0.707107
 
 julia> proj(v)
-3×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
-  [1, 1] = 0.5
-  [3, 1] = 0.5
-  [1, 3] = 0.5
-  [3, 3] = 0.5
+3×3 SparseArrays.SparseMatrixCSC{Float64,Int64} with 4 stored entries:
+  [1, 1]  =  0.5
+  [3, 1]  =  0.5
+  [1, 3]  =  0.5
+  [3, 3]  =  0.5
 ```
 """
-function proj(index::Int, size::Int)
-  ketbra(index, index, size)
-end
-
 function proj(vector::AbstractVector)
   vector*vector'
 end
@@ -114,30 +116,30 @@ end
     res(matrix)
 
 Return vectorization of the `matrix` in the row order. This is equivalent to
-`Base.vec(transpose(matrix)`.
+`Base.vec(transpose(matrix))`.
 
 # Examples
 
 ```jldoctest
 julia> M = reshape(1:9, (3, 3))'*1im
-3×3 Array{Complex{Int64}, 2}:
+3×3 Array{Complex{Int64},2}:
  0+1im  0+2im  0+3im
  0+4im  0+5im  0+6im
  0+7im  0+8im  0+9im
 
 julia> v = res(M)
-9-element Array{Complex{Int64}, 1}:
- 0+1imjulia> unres(v)
- 0+2im
- 0+3im
- 0+4im
- 0+5im
- 0+6im
- 0+7im
- 0+8im
- 0+9im
-
-julia> res(unres(v)) ==  v
+9-element reshape(::LinearAlgebra.Transpose{Complex{Int64},Array{Complex{Int64},2}}, 9) with eltype Complex{Int64}:
+ 0 + 1im
+ 0 + 2im
+ 0 + 3im
+ 0 + 4im
+ 0 + 5im
+ 0 + 6im
+ 0 + 7im
+ 0 + 8im
+ 0 + 9im
+
+julia> res(unres(v)) == v
 true
 ```
 """
@@ -155,31 +157,30 @@ have perfect square number of arguments.
 # Examples
 ```jldoctest
 julia> v = collect(1:9)*im
-9-element Array{Complex{Int64}, 1}:
- 0+1im
- 0+2im
- 0+3im
- 0+4im
- 0+5im
- 0+6im
- 0+7im
- 0+8im
- 0+9im
+9-element Array{Complex{Int64},1}:
+ 0 + 1im
+ 0 + 2im
+ 0 + 3im
+ 0 + 4im
+ 0 + 5im
+ 0 + 6im
+ 0 + 7im
+ 0 + 8im
+ 0 + 9im
 
 julia> unres(v)
-3×3 Array{Complex{Int64}, 2}:
+3×3 LinearAlgebra.Transpose{Complex{Int64},Array{Complex{Int64},2}}:
  0+1im  0+2im  0+3im
  0+4im  0+5im  0+6im
  0+7im  0+8im  0+9im
 
-julia> res(unres(v)) ==  v
+julia> res(unres(v)) == v
 true
 ```
 """
 function unres(vector::AbstractVector)
-
   dim = floor(Int64, sqrt(length(vector)))
-  @argumentcheck dim*dim ==  length(vector) "Expected vector with perfect square number of elements."
+  @argumentcheck dim*dim == length(vector) "Expected vector with perfect square number of elements."
 
   transpose(reshape(vector, (dim, dim)))
 end
@@ -195,15 +196,16 @@ Returns Fourier matrix of size `size`×`size`.
 
 ```jldoctest
 julia> fourier_matrix(1)
-1×1 SparseMatrixCSC{Complex{Float64}, Int64} with 1 stored entry:
-  [1, 1] = 1.0+0.0im
+1×1 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 1 stored entry:
+  [1, 1]  =  1.0+0.0im
 
 julia> fourier_matrix(2)
-2×2 SparseMatrixCSC{Complex{Float64}, Int64} with 4 stored entries:
-  [1, 1] = 1.0+0.0im
-  [2, 1] = 1.0+0.0im
-  [1, 2] = 1.0+0.0im
-  [2, 2] = -1.0+1.22465e-16im
+2×2 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 4 stored entries:
+  [1, 1]  =  1.0+0.0im
+  [2, 1]  =  1.0+0.0im
+  [1, 2]  =  1.0+0.0im
+  [2, 2]  =  -1.0+1.22465e-16im
+
 ```
 """
 function fourier_matrix(size::Int)
diff --git a/src/evolution.jl b/src/evolution.jl
index e319ac4..35ba9ac 100644
--- a/src/evolution.jl
+++ b/src/evolution.jl
@@ -21,30 +21,28 @@ julia> H, L = [0 1; 1 0], [[0 1; 0 0], [0 0; 1 0]]
 
 julia> evolve_operator(evolve_generator(H, L), 4.0)
 4×4 Array{Complex{Float64},2}:
- 0.499815+0.0im                0.0+0.00127256im         0.0-0.00127256im  0.500185+0.0im
-      0.0+0.00127256im  0.00960957+0.0im         0.00870607+0.0im              0.0-0.00127256im
-      0.0-0.00127256im  0.00870607+0.0im         0.00960957+0.0im              0.0+0.00127256im
- 0.500185+0.0im
+ 0.499815+0.0im                0.0+0.00127256im  …  0.500185+0.0im
+      0.0+0.00127256im  0.00960957+0.0im                 0.0-0.00127256im
+      0.0-0.00127256im  0.00870607+0.0im                 0.0+0.00127256im
+ 0.500185+0.0im                0.0-0.00127256im     0.499815+0.0im
 ```
 """
 function evolve_operator(evo_gen::AbstractMatrix{<:Number}, time::Real)
   @argumentcheck time>= 0 "Time has to be nonnegative"
-  @argumentcheck size(evo_gen, 1) ==  size(evo_gen, 2) "Argument evo_gen has to be a square matrix"
+  @argumentcheck size(evo_gen, 1) == size(evo_gen, 2) "Argument evo_gen has to be a square matrix"
 
   exp(time*evo_gen)
 end
 
 """
 
-    evolve(evo_gen, init_state, time)
-
-    evolve(evo_gen, init_state, tpoints)
-
-    evolve(evo_super, init_state)
+    evolve(evogen, state, time)
+    evolve(evogen, state, tpoints)
+    evolve(evosuper, state)
 
 Simulate the GKSL master equation according to the equation
 
-``|result⟩⟩ = exp(time*evo_gen)|init_state⟩⟩``
+``|result⟩⟩ = exp(time*evogen)|state⟩⟩``
 
 where ``|⋅⟩⟩`` denotes the vectorization.
 
@@ -52,25 +50,25 @@ where ``|⋅⟩⟩`` denotes the vectorization.
 
 The evolution can be calculated using three different approaches.
 
-In the simplest case the function accepts matrix `evo_gen` specifying the
-generator of the evolution, `init_state` describing the starting point of the
+In the simplest case the function accepts matrix `evogen` specifying the
+generator of the evolution, `state` describing the starting point of the
 evolution, and `time` specifying the time of the evolution.
 
-*Note:* If `evo_gen` is of type `Matrix`, the exponent is calculated using
-`expm` function. If `evo_gen` is of type `SparseMatrixCSC`, `expmv` from
+*Note:* If `evogen` is of type `Matrix`, the exponent is calculated using
+`exp` function. If `evogen` is of type `SparseMatrixCSC`, `expmv` from
 `Expokit.jl` is used.
 
 Alternatively, a list of point of time (`tpoints`) can be given. Points of time
 needs to be non-negative (you cannot go back in time). In this case a list of
 resulting states is returned.
 
-*Note:* It is up to the user to provide `evo_gen` and `init_state` fulfilling
-the appropriate conditions. For the procedure to work correctly `evo_gen` should
-be generated by `evolve_generator` function and `iniit_state` should be a proper
+*Note:* It is up to the user to provide `evogen` and `state` fulfilling
+the appropriate conditions. For the procedure to work correctly `evogen` should
+be generated by `evolve_generator` function and `state` should be a proper
 density matrix.
 
 The third approach can be used if the superoperator is known. In this case
-argument `evo_super` can be specified. This argument can be generated by
+argument `evosuper` can be specified. This argument can be generated by
 `evolve_operator` function. This is useful to simulate a fixed model of
 evolution in the case of multiple initial states and the same time point.
 
@@ -78,76 +76,70 @@ evolution in the case of multiple initial states and the same time point.
 
 ```jldoctest
 julia> H, L = [0 1; 1 0], [[0 1; 0 0], [0 0; 1 0]]
-(
-[0 1; 1 0],
-
-Array{Int64, 2}[
-[0 1; 0 0],
-
-[0 0; 1 0]])
+([0 1; 1 0], Array{Int64,2}[[0 1; 0 0], [0 0; 1 0]])
 
-julia> evolve(evolve_generator(H, L), proj(1, 2), 4.)
-2×2 Array{Complex{Float64}, 2}:
- 0.499815-0.0im              0.0-0.00127256im
-      0.0+0.00127256im  0.500185-0.0im
+julia> Matrix(evolve(evolve_generator(H, L), proj(1, 2), 4.))
+2×2 Array{Complex{Float64},2}:
+ 0.499815+0.0im              0.0+0.00127256im
+      0.0-0.00127256im  0.500185+0.0im
 
-julia> evolve(evolve_generator(H, L), proj(1, 2), [1., 2., 3., 4.])
-4-element Array{Array{Complex{Float64}, 2}, 1}:
- Complex{Float64}[0.433203-0.0im 0.0-0.107605im; 0.0+0.107605im 0.566797-0.0im]
- Complex{Float64}[0.485766-0.0im 0.0+0.0171718im; 0.0-0.0171718im 0.514234-0.0im]
- Complex{Float64}[0.505597-0.0im 0.0+0.00261701im; 0.0-0.00261701im 0.494403-0.0im]
- Complex{Float64}[0.499815-0.0im 0.0-0.00127256im; 0.0+0.00127256im 0.500185-0.0im]
+julia> Matrix.(evolve(evolve_generator(H, L), proj(1, 2), [1., 2., 3., 4.]))
+4-element Array{Array{Complex{Float64},2},1}:
+ [0.433203+0.0im 0.0+0.107605im; 0.0-0.107605im 0.566797+0.0im]
+ [0.485766+0.0im 0.0-0.0171718im; 0.0+0.0171718im 0.514234+0.0im]
+ [0.505597+0.0im 0.0-0.00261701im; 0.0+0.00261701im 0.494403+0.0im]
+ [0.499815+0.0im 0.0+0.00127256im; 0.0-0.00127256im 0.500185+0.0im]
 
 julia> ev_op = evolve_operator(evolve_generator(H, L), 4.)
-4×4 Array{Complex{Float64}, 2}:
- 0.499815+0.0im                0.0+0.00127256im         0.0-0.00127256im  0.500185+0.0im
-      0.0+0.00127256im  0.00960957+0.0im         0.00870607+0.0im              0.0-0.00127256im
-      0.0-0.00127256im  0.00870607+0.0im         0.00960957+0.0im              0.0+0.00127256im
- 0.500185+0.0im                0.0-0.00127256im         0.0+0.00127256im  0.499815+0.0im
-
-julia> evolve(ev_op, proj(1, 2))
-2×2 Array{Complex{Float64}, 2}:
+4×4 Array{Complex{Float64},2}:
+ 0.499815+0.0im                0.0+0.00127256im  …  0.500185+0.0im
+      0.0+0.00127256im  0.00960957+0.0im                 0.0-0.00127256im
+      0.0-0.00127256im  0.00870607+0.0im                 0.0+0.00127256im
+ 0.500185+0.0im                0.0-0.00127256im     0.499815+0.0im
+
+julia> Matrix(evolve(ev_op, proj(1, 2)))
+2×2 Array{Complex{Float64},2}:
  0.499815+0.0im              0.0+0.00127256im
       0.0-0.00127256im  0.500185+0.0im
 ```
 """
 function evolve(exp_evolve_generator::AbstractMatrix{<:Number},
                 initial_state::AbstractMatrix{<:Number})
-  @argumentcheck size(exp_evolve_generator, 1) ==  size(exp_evolve_generator, 2) "Argument exp_evolve_generator should be square"
-  @argumentcheck size(initial_state, 1) ==  size(initial_state, 2) "Initial_state should be a square matrix"
-  @assert size(exp_evolve_generator, 1) ==  size(initial_state, 1)^2 "The initial state size should be square root of exp_evolve_generator size"
+  @argumentcheck size(exp_evolve_generator, 1) == size(exp_evolve_generator, 2) "Argument exp_evolve_generator should be square"
+  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Initial_state should be a square matrix"
+  @assert size(exp_evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of exp_evolve_generator size"
 
   unres(exp_evolve_generator*res(initial_state))
-end
+end,
 
 function evolve(evolve_generator::AbstractMatrix{<:Number},
                 initial_state::AbstractMatrix{<:Number},
                 timepoint::Real)
-  @argumentcheck size(evolve_generator, 1) ==  size(evolve_generator, 2) "Argument evolve_generator should be square"
-  @argumentcheck size(initial_state, 1) ==  size(initial_state, 2) "Initial_state should be a square matrix"
-  @assert size(evolve_generator, 1) ==  size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
+  @argumentcheck size(evolve_generator, 1) == size(evolve_generator, 2) "Argument evolve_generator should be square"
+  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Initial_state should be a square matrix"
+  @assert size(evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
   @argumentcheck timepoint>= 0 "Time needs to be nonnegative"
 
   unres(exp(timepoint*evolve_generator)*res(initial_state))
-end
+end,
 
 function evolve(evolve_generator::SparseMatrixCSC{<:Number},
                 initial_state::AbstractMatrix{<:Number},
                 timepoint::Real)
-  @argumentcheck size(evolve_generator, 1) ==  size(evolve_generator, 2) "Argument evolve_generator should be a square matrix"
-  @argumentcheck size(initial_state, 1) ==  size(initial_state, 2) "Argument initial_state should be a square matrix"
-  @assert size(evolve_generator, 1) ==  size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
+  @argumentcheck size(evolve_generator, 1) == size(evolve_generator, 2) "Argument evolve_generator should be a square matrix"
+  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Argument initial_state should be a square matrix"
+  @assert size(evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
   @argumentcheck timepoint>= 0 "Time needs to be nonnegative"
 
   unres(expmv(timepoint, evolve_generator, Vector(res(initial_state))))
-end
+end,
 
 function evolve(evolve_generator::AbstractMatrix{<:Number},
                 initial_state::AbstractMatrix{<:Number},
                 timepoints::AbstractVector{<:Real})
-  @argumentcheck size(evolve_generator, 1) ==  size(evolve_generator, 2) "evolve_generator should be a square matrix"
-  @argumentcheck size(initial_state, 1) ==  size(initial_state, 2) "Argument initial_state should be a square matrix"
-  @assert size(evolve_generator, 1) ==  size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
+  @argumentcheck size(evolve_generator, 1) == size(evolve_generator, 2) "evolve_generator should be a square matrix"
+  @argumentcheck size(initial_state, 1) == size(initial_state, 2) "Argument initial_state should be a square matrix"
+  @assert size(evolve_generator, 1) == size(initial_state, 1)^2 "The initial state size should be square root of evolve_generator size"
   @argumentcheck all(timepoints.>= 0) "All time points need to be nonnegative"
 
   [evolve(evolve_generator, initial_state, t) for t = timepoints]
diff --git a/src/operator.jl b/src/operator.jl
index c7af6cc..358a944 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -6,7 +6,7 @@ export
     local_lind(A[; epsilon])
 
 Split the elements of the matrix `A` into a collection of sparse matrices with
-exactly one non-zero element. Martices are added if the absolute value of the
+exactly one non-zero element. Martices are created if the absolute value of the
 nonzero element is there are not smaller than `epsilon`, where `epsilon` should
 be nonnegative. The `epsilon` defaults to `eps()` if not specified.
 
@@ -14,29 +14,30 @@ be nonnegative. The `epsilon` defaults to `eps()` if not specified.
 
 ```jldoctest
 julia> A = [1. 2.; 3. 4.]
-2×2 Array{Float64, 2}:
+2×2 Array{Float64,2}:
  1.0  2.0
  3.0  4.0
 
 julia> local_lind(A)
-4-element Array{SparseMatrixCSC{Float64, Ti<:Integer}, 1}:
+4-element Array{SparseArrays.SparseMatrixCSC{Float64,Ti} where Ti<:Integer,1}:
 
-	[1, 1] = 1.0
+  [1, 1]  =  1.0
 
-	[1, 2] = 2.0
+  [1, 2]  =  2.0
 
-	[2, 1] = 3.0
+  [2, 1]  =  3.0
 
-	[2, 2] = 4.0
+  [2, 2]  =  4.0
 
 julia> local_lind(A, epsilon = 1.5)
-3-element Array{SparseMatrixCSC{Float64, Ti<:Integer}, 1}:
+3-element Array{SparseArrays.SparseMatrixCSC{Float64,Ti} where Ti<:Integer,1}:
 
-	[1, 2] = 2.0
+  [1, 2]  =  2.0
 
-	[2, 1] = 3.0
+  [2, 1]  =  3.0
+
+  [2, 2]  =  4.0
 
-	[2, 2] = 4.0
 ```
 """
 function local_lind(A::Matrix{T}; epsilon::Real = eps()) where T<:Number
@@ -100,9 +101,9 @@ function evolve_generator_create(H::AbstractMatrix{<:Number},
                                  α::Real,
                                  β::Real)
   @argumentcheck size(H) !=  (0, 0) "Matrix H must not be sizeless"
-  @argumentcheck size(H, 1) ==  size(H, 2) "Matrix H must be square"
-  @assert all([size(lindbladian) ==  size(H) for lindbladian in L]) "Lindblad operators must be of the same size as Hamiltonian"
-  @assert size(H) ==  size(localH) "Matrix localH must be of the same size as H"
+  @argumentcheck size(H, 1) == size(H, 2) "Matrix H must be square"
+  @assert all([size(lindbladian) == size(H) for lindbladian in L]) "Lindblad operators must be of the same size as Hamiltonian"
+  @assert size(H) == size(localH) "Matrix localH must be of the same size as H"
   @argumentcheck 0 <=  α <=  1 && 0 <=  β <=  1 "Value of ω must be nonngeative and smaller than one"
 
   F = spzeros(ComplexF64, (size(H).^2)...)
@@ -136,8 +137,7 @@ operator takes the form
 # Arguments
 - `H`: Hamiltonian, must be hermitian,
 - `L`: collection of Lindblad operators, each must be of the same size as `H`,
-- `localH`: local Hamiltonian, suggested for nonmoralized QS walk, must be hermitian
-and of the size of `H`,
+- `localH`: local Hamiltonian, suggested for nonmoralized QS walk, must be hermitian and of the size of `H`,
 - `ω`: scaling parameter, should be in [0, 1].
 
 # Return
@@ -147,16 +147,11 @@ The generator matrix, which can be used in `evolve` function.
 # Examples
 
 ```jldoctest
-julia> H, L, localH = [0 1+im; 1-im 0], [0. 1; 0 0], eye(2)
-(
-Complex{Int64}[0+0im 1+1im; 1-1im 0+0im],
-
-[0.0 1.0; 0.0 0.0],
-
-[1.0 0.0; 0.0 1.0])
+julia> H, L, localH = [0 1+im; 1-im 0], [0. 1; 0 0], [1. 0.; 0. 1.]
+(Complex{Int64}[0+0im 1+1im; 1-1im 0+0im], [0.0 1.0; 0.0 0.0], [1.0 0.0; 0.0 1.0])
 
 julia> evolve_generator(H, [L], localH, 1/2)
-4×4 Array{Complex{Float64}, 2}:
+4×4 Array{Complex{Float64},2}:
   0.0+0.0im    0.5+0.5im    0.5-0.5im   0.5+0.0im
  -0.5+0.5im  -0.25+0.0im    0.0+0.0im   0.5-0.5im
  -0.5-0.5im    0.0+0.0im  -0.25+0.0im   0.5+0.5im
@@ -169,13 +164,13 @@ function evolve_generator(H::AbstractMatrix{<:Number},
                           localH::AbstractMatrix{<:Number},
                           ω::Real)
   evolve_generator_create(H, L, localH, 1-ω, ω)
-end
+end,
 
 function evolve_generator(H::AbstractMatrix{<:Number},
                           L::AbstractVector{<:AbstractMatrix{<:Number}},
                           localH::AbstractMatrix{<:Number} = spzeros(eltype(H), size(H)...))
   evolve_generator_create(H, L, localH, 1., 1.)
-end
+end,
 
 function evolve_generator(H::AbstractMatrix{T},
                           L::AbstractVector{<:AbstractMatrix{<:Number}},
diff --git a/src/utils.jl b/src/utils.jl
index 1827cdb..f03b7c1 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -5,7 +5,8 @@ export
   VertexSet,
   vertexsetsize,
   subspace,
-  vlist
+  vlist,
+  length
 
 import Base: ==, hash, getindex, length
 
@@ -44,14 +45,31 @@ struct Vertex
   Vertex(v::Vector{Int}) = all(v.>0) ? new(v) : throw(ArgumentError("Vector should consist of positive elments"))
 end
 
+"""
+    subspace(v)
+
+Returns the subspace connected to vertex `v`.
+
+```jldoctest
+julia> v = Vertex([1,2,3])
+Vertex([1, 2, 3])
+
+julia> subspace(v)
+3-element Array{Int64,1}:
+ 1
+ 2
+ 3
+```
+"""
 subspace(v::Vertex) = v.subspace
 
-==(v::Vertex, w::Vertex) = subspace(v) ==  subspace(w)
+==(v::Vertex, w::Vertex) = subspace(v) == subspace(w)
 hash(v::Vertex) = hash(subspace(v))
 
-==(v::Tuple{Vertex, Vertex}, w::Tuple{Vertex, Vertex}) = [subspace(v[1]), subspace(v[2])] ==  [subspace(w[1]), subspace(w[2])]
+==(v::Tuple{Vertex, Vertex}, w::Tuple{Vertex, Vertex}) = [subspace(v[1]), subspace(v[2])] == [subspace(w[1]), subspace(w[2])]
 hash(v::Tuple{Vertex, Vertex}) = hash([subspace(v[1]), subspace(v[2])])
 getindex(v::Vertex, i::Int) = v.subspace[i]
+
 length(v::Vertex) = length(subspace(v))
 
 
@@ -60,7 +78,7 @@ function checkvertexset(partition::Vector{Vector{Int}})
   for lin = partition
     append!(joined, lin)
   end
-  length(joined) ==  length(Set(joined))
+  length(joined) == length(Set(joined))
 end
 
 checkvertexset(vset::Vector{Vertex}) = checkvertexset([subspace(v) for v = vset])
@@ -72,7 +90,7 @@ Type consisting of a list of `Vertex` objects. It describes the partition of the
 linear subspace. Object of this type should be constructed using
 `make_vertex_set` or by `nm_lind` functions. In order to get a
 list of the vertices from an object of type `vertexset`, one should use
-`vertexset()` function, or, for a specific `Vertex`, an getindex function
+`vlist()` function, or, for a specific `Vertex`, an getindex function
 `vertexset[i]`.
 """
 struct VertexSet
@@ -83,9 +101,26 @@ end
 
 VertexSet(vset::Vector{Vector{Int}}) = VertexSet([Vertex(v) for v = vset])
 
+"""
+    vlist(vset)
+
+Returns the list of vertices for given `vset` of type `VertexSet`.
+
+```jldoctest
+julia> vset = VertexSet([[1], [2, 3, 4], [5, 6, 7], [8, 9]])
+VertexSet(Vertex[Vertex([1]), Vertex([2, 3, 4]), Vertex([5, 6, 7]), Vertex([8, 9])])
+
+julia> vlist(vset)
+4-element Array{Vertex,1}:
+ Vertex([1])
+ Vertex([2, 3, 4])
+ Vertex([5, 6, 7])
+ Vertex([8, 9])
+```
+"""
 vlist(vset::VertexSet) = vset.vertices
 
-==(v::VertexSet, w::VertexSet) = v.vertices ==  w.vertices
+==(v::VertexSet, w::VertexSet) = v.vertices == w.vertices
 
 getindex(vset::VertexSet, i::Int) = vlist(vset)[i]
 getindex(vset::VertexSet, veci::Vector{Int}) = vlist(vset)[veci]
@@ -94,7 +129,7 @@ length(vset::VertexSet) = length(vlist(vset))
 """
     vertexsetsize(vertexset)
 
-Return the dimenion of the linearspace corresponding to given `vertexset'.
+Return the dimension of the linearspace corresponding to given `vertexset`.
 
 # Examples
 
@@ -224,34 +259,33 @@ end
     make_vertex_set(A[, epsilon])
 
 Creates object of type `VertexSet` which represents how vertices are located in
-matrix. Should be used only in the nondefault use of `global_operator` function.
-It is always equal to the second element if output of `global_operator` function.
+matrix. Should be used only in the nondefault use of `evolve_generator` function.
+It is always equal to the second element if output of `evolve_generator` function.
 
 # Examples
 
 ```jldoctest
 julia> A = [1 2 3; 0 3. 4.; 0 0 5.]
-3×3 Array{Float64, 2}:
+3×3 Array{Float64,2}:
  1.0  2.0  3.0
  0.0  3.0  4.0
  0.0  0.0  5.0
 
-julia> make_vertex_set(A)()
-3-element Array{QSWalk.Vertex, 1}:
- QSWalk.Vertex([1, 2, 3])
- QSWalk.Vertex([4, 5])
- QSWalk.Vertex([6])
-
-julia> make_vertex_set(A, epsilon = 2.5)()
-3-element Array{QSWalk.Vertex, 1}:
- QSWalk.Vertex([1])
- QSWalk.Vertex([2, 3])
- QSWalk.Vertex([4])
-
+julia> vlist(make_vertex_set(A))
+3-element Array{Vertex,1}:
+ Vertex([1, 2, 3])
+ Vertex([4, 5])
+ Vertex([6])
+
+julia> vlist(make_vertex_set(A, epsilon = 2.5))
+3-element Array{Vertex,1}:
+ Vertex([1])
+ Vertex([2, 3])
+ Vertex([4])
 ```
 """
 function make_vertex_set(A::AbstractMatrix; epsilon::Real = eps())
-  @argumentcheck epsilon >=  0 "epsilon needs to be nonnegative"
-  @argumentcheck size(A, 1) ==  size(A, 2) "A matrix must be square"
+  @argumentcheck epsilon >= 0 "epsilon needs to be nonnegative"
+  @argumentcheck size(A, 1) == size(A, 2) "A matrix must be square"
   revinc_to_vertexset(reversed_incidence_list(A, epsilon = epsilon))
 end
diff --git a/test/dirac.jl b/test/dirac.jl
index 71dfa21..30b7ca3 100644
--- a/test/dirac.jl
+++ b/test/dirac.jl
@@ -17,7 +17,7 @@
 
   @testset "ketbra" begin
     #standard tests
-    @test ketbra(1, 2, 3) ==  [0 1 0; 0 0 0; 0 0 0]
+    @test ketbra(1, 2, 3) == [0 1 0; 0 0 0; 0 0 0]
     #error tests
     @test_throws AssertionError ketbra(3, 2, 2)
     @test_throws AssertionError ketbra(2, 3, 2)
@@ -29,7 +29,7 @@
     result = [0.0+0.0im 0.0+0.0im 0.0+0.0im;
               0.0+0.0im 1.0+0.0im 0.0+0.0im;
               0.0+0.0im 0.0+0.0im 0.0+0.0im]
-    @test proj(2, 3) ==  result
+    @test proj(2, 3) == result
     @test proj(1/sqrt(2) * (ket(1, 3)+ket(3, 3))) ≈
                          [0.5+0.0im 0.0+0.0im 0.5+0.0im;
                           0.0+0.0im 0.0+0.0im 0.0+0.0im;