diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..b0506e2
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,15 @@
+using Documenter, SpinGlassNetworks
+
+_pages = [
+    "Introduction" => "index.md",
+    "User guide" => "userguide.md",
+    "API Reference" => "api.md"
+]
+
+# ============================
+
+makedocs(
+    sitename="SpinGlassNetworks",
+    modules = [SpinGlassNetworks],
+    pages = _pages
+    )
\ No newline at end of file
diff --git a/docs/src/api.md b/docs/src/api.md
new file mode 100644
index 0000000..8a1c80a
--- /dev/null
+++ b/docs/src/api.md
@@ -0,0 +1,36 @@
+# Library
+
+---
+
+```@meta
+CurrentModule = SpinGlassNetworks
+```
+
+## Ising Graphs
+
+```@docs
+ising_graph
+biases
+couplings
+cluster
+prune
+rank_vec
+```
+
+## Factor Graphs
+
+```@docs
+factor_graph
+decode_factor_graph_state
+```
+
+## Auxiliary Functions
+
+```@docs
+brute_force
+energy
+full_spectrum
+gibbs_tensor
+rank_reveal
+super_square_lattice
+```
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..7fd6448
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,6 @@
+# SpinGlassNetworks
+
+A [Julia](http://julialang.org) package for building and interacting with Ising spin glass models in context of tensor networks. Part of [SpinGlassPEPS](https://github.com/euro-hpc-pl/SpinGlassPEPS.jl) package.
+
+
+
diff --git a/docs/src/userguide.md b/docs/src/userguide.md
new file mode 100644
index 0000000..20f33b5
--- /dev/null
+++ b/docs/src/userguide.md
@@ -0,0 +1,122 @@
+# Introduction
+This section provides examples of how to use the 'ising_graph' function to generate an Ising model graph and the 'factor_graph' function to convert the Ising graph into a factor graph.
+
+# Ising Graph
+The Ising model is a mathematical model used to describe the behavior of interacting particles, such as atoms or molecules, in a magnetic field. In the Ising model, each particle is represented as a binary variable $s_i$ that can take on the values of either +1 or -1. The total energy of the system is given by the Hamiltonian:
+
+$$H =  \sum_{(i,j) \in \mathcal{E}} J_{ij} s_i s_j + \sum_{i} h_i s_i$$
+
+where $J_{ij}$ is the coupling constant between particles $i$ and $j$, $h_i$ is the external magnetic field at particle $i$, and the sum is taken over all pairs of particles and all particles in the system $\mathcal{E}$, respectively.
+
+
+## simple examle
+```@example
+using SpinGlassNetworks
+# Create Ising instance
+instance = Dict((1, 1) => 1.0, (2, 2) => 0.5, (3, 3) => -0.25, 
+(1, 2) => -1.0, (2, 3) => 1.0)
+
+# Generate Ising graph
+ig = ising_graph(instance)
+
+# View graph properties
+@show biases(ig), couplings(ig)
+
+```
+
+## Generalisation
+
+The 'ising_graph' function is quite general function. One can easly convert instance into another hamiltonian convention, ex.
+```math
+H = - \sum_{(i,j) \in \mathcal{E}} J_{ij} s_i s_j - \mu \sum_{i} h_i s_i
+```
+This conversion is directed by `sgn` parameter of `ising_graph`. It is wort doing, because then all `energy` functions will work as one intended.
+```@example
+using SpinGlassNetworks, LabelledGraphs
+
+# lets define simple Ising system
+instance = Dict((1,1) => 1.0, (2, 2) => -1.0, (1, 2) => -1.25)
+ig = ising_graph(instance)
+
+# use above convention
+μ = 0.25
+sgn = -1
+
+# absorb μ into h and construct ising graph and 
+instance2 = Dict((1,1) => instance[(1,1)]*μ, (2, 2) => instance[(2,2)]*μ, (1, 2) => -1.25)
+
+ig2 = ising_graph(instance2, sgn=sgn)
+
+# select arbitrary state
+σ = [1, 1]
+
+println("energy without conversion: ", energy(σ, ig), 
+" energy after conversion: ", energy(σ, ig2))
+```
+
+'ising_graph' also allows for building other spin system, not only spin-``\frac{1}{2}`` system. For example one may build chain where sites allows for 3 spin values.
+
+```@example
+using SpinGlassNetworks
+
+# Create spin chain instance 
+instance = Dict((1, 1) => 1.0, (2, 2) => 0.5, (3, 3) => -0.25, 
+(1, 2) => -1.0, (2, 3) => 1.0)
+
+rank_override = Dict(1 => 3, 2 => 3, 3 => 3)
+
+sg = ising_graph(instance, rank_override=rank_override)
+
+# see ranks of each site
+rank_vec(sg)
+```
+
+# factor graph
+
+A factor graph is a graphical representation that allows for a convenient and intuitive way to describe the structure of a tensor network.
+
+A factor graph consists of two types of nodes: variable nodes and factor nodes. Variable nodes represent the tensors in the tensor network, and factor nodes represent the operations that are applied to those tensors.
+
+Each variable node in the factor graph corresponds to a tensor in the tensor network, and each factor node corresponds to a tensor contraction or operation that is applied to one or more of those tensors. The edges between the nodes in the factor graph represent the indices that are shared between the corresponding tensors in the tensor network.
+
+## Simple example
+
+```@example
+using SpinGlassNetworks
+
+# Prepare simple instance
+instance = Dict((1, 1) => 1.0, (2, 2) => 0.5, (3, 3) => -0.25, 
+(1, 2) => -1.0, (2, 3) => 1.0)
+ig = ising_graph(instance)
+
+# Create factor graph.
+fg = factor_graph(
+    ig,
+    cluster_assignment_rule = super_square_lattice((3, 1, 1))
+)
+```
+
+## Chimera graphs
+The Chimera graph is a type of graph architecture used in quantum computing systems, particularly in the quantum annealing machines developed by D-Wave Systems. It is a two-dimensional lattice of unit cells, each consisting of a bipartite graph of $K_{4,4}$ complete bipartite subgraphs. Futer details can be found [here](https://docs.dwavesys.com/docs/latest/c_gs_4.html#chimera-graph).
+
+
+```@example
+using SpinGlassEngine, SpinGlassNetworks, LabelledGraphs
+
+# load Chimera instance and create Ising graph
+instance = "$(@__DIR__)/../../src/instances/chimera_droplets/128power/001.txt"
+ig = ising_graph(instance)
+
+# Loaded instance is 4x4x8 chimera graph
+m = 4
+n = 4
+t = 8
+
+fg = fg = factor_graph(
+    ig,
+    cluster_assignment_rule = super_square_lattice((m, n, t))
+)
+
+println("Number of nodes in oryginal instance: ", length(LabelledGraphs.vertices(ig)), "\n",
+        " Number of nodes in factor graph: ", length(LabelledGraphs.vertices(fg)))
+```
\ No newline at end of file
diff --git a/src/factor.jl b/src/factor.jl
index 06292f9..c52bf24 100644
--- a/src/factor.jl
+++ b/src/factor.jl
@@ -27,7 +27,28 @@ function factor_graph(
     )
 end
 
-function factor_graph(
+"""
+    factor_graph(
+        ig::IsingGraph,
+        num_states_cl::Dict;
+        spectrum::Function = full_spectrum,
+        cluster_assignment_rule::Dict) where {T}
+
+Constructs a factor graph representation of the given `IsingGraph`. 
+
+# Arguments
+- `ig::IsingGraph`: The Ising graph to convert to a factor graph.
+- `num_states_cl::Dict` : A dictionary mapping each cluster to the number of states it can take on. If given empty dictionary
+    function will try to interfere number of states for each cluster. Can be also given as `Int`, then each cluster will have 
+    specified number of states. If ommited completly, function will behave as if an empy directory was passed.
+- `spectrum::Function`: A function that computes the spectrum (i.e., list of energies of all possible states) of
+  a given cluster. The default is `full_spectrum`, which computes the spectrum exactly.
+- `cluster_assignment_rule::Dict`: A dictionary that assigns each vertex of the Ising graph to a cluster.
+
+# Output
+A `LabelledGraph` that represents the factor graph of the Ising graph.
+"""
+ function factor_graph(
     ig::IsingGraph,
     num_states_cl::Dict{T,Int};
     spectrum::Function = full_spectrum,
@@ -64,7 +85,8 @@ function factor_graph(
         end
     end
     fg
-end
+end   
+
 
 function factor_graph(
     ig::IsingGraph;
@@ -79,6 +101,15 @@ function factor_graph(
     )
 end
 
+"""
+    rank_reveal(energy, order = :PE)
+
+calculate the rank of the matrix represented by `energy`, and returns a tuple of two matrices `P` and `E`, 
+where `P` is a binary matrix that reveals the rank of `energy`, and `E` is the set of non-zero singular 
+values of `energy` in decreasing order. Argument `order=:PE` specifies the order of the output matrices. 
+The default is `:PE`, which means that the `P` matrix is returned first, followed by the `E` matrix. 
+If `order` is set to `:EP`, the order is reversed.
+"""
 function rank_reveal(energy, order = :PE)
     @assert order ∈ (:PE, :EP)
     dim = order == :PE ? 1 : 2
@@ -98,7 +129,11 @@ function rank_reveal(energy, order = :PE)
     order == :PE ? (P, E) : (E, P)
 end
 
+"""
+    decode_factor_graph_state(fg, state::Vector{Int})
 
+Convert factor graphs clusters states into state of the original Ising system.
+"""
 function decode_factor_graph_state(fg, state::Vector{Int})
     ret = Dict{Int,Int}()
     for (i, vert) ∈ zip(state, vertices(fg))
diff --git a/src/ising.jl b/src/ising.jl
index dfc3d3b..17fb1f4 100644
--- a/src/ising.jl
+++ b/src/ising.jl
@@ -9,16 +9,29 @@ const Instance = Union{String,Dict}
 unique_nodes(ising_tuples) =
     sort(collect(Set(Iterators.flatten((i, j) for (i, j, _) ∈ ising_tuples))))
 
+
+
 const IsingGraph = LabelledGraph{MetaGraph{Int64,Float64},Int64}
 
 """
-$(TYPEDSIGNATURES)
+    ising_graph(instance, sgn::Number = 1.0, rank_override::Dict{Int,Int} = Dict{Int,Int}())
+
+Create an Ising model graph from an input `instance``.
 
-Create the Ising spin glass model.
+# Arguments
+- `instance`: Instance of the Ising spin glass system. Can be given eiter as dictionary or a CSV file.
+- `sgn::Number`: decides the used convention for the energy sum
+- `rank_override`: Dict{Int,Int}
 
-# Details
+# Output
+A `LabelledGraph` that represents inputed Ising instance.
 
-Store extra information
+# Example
+```@example
+instance =  Dict((1, 1) => 0.5, (2, 2) => 0.75, (3, 3) => -0.25, (1, 2) => -1.0, (2, 3) => 1.0)
+
+ig = ising_graph(instance)
+```
 """
 function ising_graph(
     instance::Instance;
@@ -58,10 +71,27 @@ function ising_graph(
     ig
 end
 
+"""
+    rank_vec(ising_graph)
+
+Return vector of 'ranks' of nodes. In this context rank tels us how many spins values node can have. 
+"""
 rank_vec(ig::IsingGraph) = Int[get_prop((ig), v, :rank) for v ∈ vertices(ig)]
 basis_size(ig::IsingGraph) = prod(prod(rank_vec(ig)))
+
+"""
+    biases(ising_graph)
+
+Return vector of biases of ising model defined by given 'ising_graph'.
+"""
 biases(ig::IsingGraph) = get_prop.(Ref(ig), vertices(ig), :h)
 
+"""
+    couplings(ising_graph)
+
+Return couplings of ising model defined by given 'ising_graph'. 
+The couplings are presented as an uper-triangular matrix 
+"""
 function couplings(ig::IsingGraph)
     J = zeros(nv(ig), nv(ig))
     for edge in edges(ig)
@@ -71,6 +101,11 @@ function couplings(ig::IsingGraph)
     J
 end
 
+"""
+    cluster(ig::IsingGraph, verts)
+
+Return cluster specified by `verts`
+"""
 cluster(ig::IsingGraph, verts) = induced_subgraph(ig, collect(verts))
 
 function inter_cluster_edges(ig::IsingGraph, cl1::IsingGraph, cl2::IsingGraph)
@@ -85,6 +120,12 @@ function inter_cluster_edges(ig::IsingGraph, cl1::IsingGraph, cl2::IsingGraph)
     outer_edges, J
 end
 
+
+"""
+    prune(ising_graph)
+
+Return copy of 'ising_graph' with removed isolated vertices.
+"""
 function prune(ig::IsingGraph)
     to_keep = vcat(
         findall(!iszero, degree(ig)),
diff --git a/src/lattice.jl b/src/lattice.jl
index 51c604c..29a4211 100644
--- a/src/lattice.jl
+++ b/src/lattice.jl
@@ -7,6 +7,14 @@ function super_square_lattice(size::NTuple{5,Int})
     Dict(old[k, uj, j, ui, i] => (i, j) for i = 1:m, ui = 1:um, j = 1:n, uj = 1:un, k = 1:t)
 end
 
+"""
+    super_square_lattice(size::NTuple{3,Int})
+
+Create cluster assignment rule for Chimera architecture. The input is asumed to have form `(m, n, t)` where `m` is number of
+columns, `n` number of rows and `t` is the size of the shore within each Chimera tile 
+([details](https://docs.ocean.dwavesys.com/en/stable/docs_dnx/reference/generated/dwave_networkx.chimera_graph.html)). 
+If `t=1` then this assigment rule can be used for square lattice, where every site will form its own unit cell.  
+"""
 function super_square_lattice(size::NTuple{3,Int})
     m, n, t = size
     super_square_lattice((m, 1, n, 1, t))
diff --git a/src/spectrum.jl b/src/spectrum.jl
index 944037f..5a7bed6 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -4,7 +4,7 @@ export brute_force, full_spectrum, energy
 export Spectrum
 
 """
-$(TYPEDSIGNATURES)
+    gibbs_tensor(ig::IsingGraph, β::Float64 = 1.0)
 
 Calculates Gibbs state of a classical Ising Hamiltonian
 
@@ -12,11 +12,11 @@ Calculates Gibbs state of a classical Ising Hamiltonian
 
 Calculates matrix elements (probabilities) of \$\\rho\$
 ```math
-\$\\bra{\\σ}\\rho\\ket{\\sigma}\$
+\\bra{\\sigma}\\rho\\ket{\\sigma}
 ```
-for all possible configurations \$\\σ\$.
+for all possible configurations \$\\sigma\$.
 """
-function gibbs_tensor(ig::IsingGraph, β = Float64 = 1.0)
+function gibbs_tensor(ig::IsingGraph, β::Float64 = 1.0)
     states = collect.(all_states(rank_vec(ig)))
     ρ = exp.(-β .* energy.(states, Ref(ig)))
     ρ ./ sum(ρ)
@@ -24,20 +24,39 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+    energy(σ::Vector, ig::IsingGraph)
 
 Calculate the Ising energy
 ```math
-E = -\\sum_<i,j> s_i J_{ij} * s_j - \\sum_j h_i s_j.
+E = \\sum_{<i,j>} σ_i J_{ij} σ_j + \\sum_j h_i σ_j.
+```
+"""
+energy(σ::Vector, ig::IsingGraph)= energy(σ, couplings(ig)) + energy(σ, biases(ig))
+
+"""
+    energy(σ::Vector, h::Vector) = dot(h, σ)
+
+Caltulate linear part of Ising energy
+```math
+\\sum_j h_i σ_j
 ```
 """
-energy(σ::Vector, J::Matrix, η::Vector = σ) = dot(σ, J, η)
 energy(σ::Vector, h::Vector) = dot(h, σ)
-energy(σ::Vector, ig::IsingGraph) = energy(σ, couplings(ig)) + energy(σ, biases(ig))
 
+"""
+    energy(σ::Vector, J::Matrix, η::Vector = σ)
 
+Calculate quadratic part of Ising energy
+```math
+\\sum_{<i,j>} σ_i J_{ij} σ_j
+```
 """
-$(TYPEDSIGNATURES)
+energy(σ::Vector, J::Matrix, η::Vector = σ) = dot(σ, J, η)
+ 
+
+
+"""
+    brute_force(ig::IsingGraph; sorted = true, num_states::Int = 1)
 
 Return the low energy spectrum
 
@@ -47,7 +66,6 @@ Calculates \$k\$ lowest energy states
 together with the coresponding energies
 of a classical Ising Hamiltonian
 """
-
 function brute_force(ig::IsingGraph; sorted = true, num_states::Int = 1)
     if nv(ig) == 0
         return Spectrum(zeros(1), [])
@@ -65,6 +83,13 @@ function brute_force(ig::IsingGraph; sorted = true, num_states::Int = 1)
     end
 end
 
+
+"""
+    full_spectrum(ig::IsingGraph; num_states::Int = 1)
+
+Calculates full energy spectrum via the brute force method. Returns `num_states` energy 
+states together with the coresponding energies
+"""
 full_spectrum(ig::IsingGraph; num_states::Int = 1) =
     brute_force(ig, sorted = false, num_states = num_states)