type array<'t, 'ndims> =
    {
        dims : int array
        data : 't array
    }

type ndims2 =
    static member n with get() = 2
    static member indexer (a:array<_, ndims2>) =
        let indexImpl (i:int, j:int) = i + j*a.dims.[1]
        indexImpl
    static member eachindex (a:array<_, ndims2>) =
        let l1 = a.dims.[0]
        let l2 = a.dims.[1]
        seq { for i in 0 .. l1 - 1 do
                for j in 0 .. l2 - 1 do
                    yield (i, j)  }

let inline indexer (a:array<_, 'ndims>) =
    (^ndims : (static member indexer : array<_, ^ndims> -> _) a)

let inline eachindex (a:array<_, 'ndims>) =
    (^ndims : (static member eachindex : array<_, ^ndims> -> _) a)

let inline get (a:array<_, 'ndims>) ind =
    let i = indexer a
    a.data.[i(ind)]

let inline set (a:array<_, 'ndims>) ind value =
    let i = indexer a
    a.data.[i(ind)] <- value

let inline mapi (f:int -> 't -> 'u) (a:array<'t, 'ndims>) =
    let index = indexer a
    // let res = ...
    // ... some boilerplate
    for ind in (eachindex a) do
        let i = index ind
        res.data.[i] <- f i a.data.[i]
    res

let m1 : array<float, ndims2> =
    {
        dims = [|2; 2|]
        data = [|1.; 0.; 0.; 1.|]
    }

mapi (fun i x -> (i, x)) m1