using LinearMapsAA
using ImagePhantoms: shepp_logan, SheppLoganToft
using MIRTjim: jim, prompt
using InteractiveUtils: versioninfo

isinteractive() ? jim(:prompt, true) : prompt(:draw);

down1 = (x) -> (x[1:2:end,:] + x[2:2:end,:])/2 # 1D down-sampling by 2×
down2 = (x) -> down1(down1(x)')'; # 2D down-sampling by factor of 2×

down2_adj(y::AbstractMatrix{<:Number}) = kron(y, fill(0.25, (2,2)));

nx, ny = 200, 256
A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny);
    idim = (nx,ny), odim = (nx,ny) .÷ 2)

image = shepp_logan(ny, SheppLoganToft())[(ny-nx)÷2 .+ (1:nx),:]
jim(image, "SheppLogan")

down = A * image
jim(down, title="Down-sampled image")

up = A' * down
jim(up, title="Adjoint: A' * y")

B = LinearMapAA(
        x -> vec(down2(reshape(x,nx,ny))),
        y -> vec(down2_adj(reshape(y,Int(nx/2),Int(ny/2)))),
        ((nx÷2)*(ny÷2), nx*ny),
    )

y = B * vec(image) # This would fail here without the `vec`!
jim(reshape(y, nx÷2, ny÷2)) # Annoying reshape needed here!

nx, ny = 8,6
idim = (nx,ny)
odim = (nx,ny) .÷ 2
A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny); idim, odim)

jim(Matrix(A)', "A")

jim(Matrix(A')', "A'")

@assert Matrix(A)' == Matrix(A')

@assert Matrix(A)' ≈ Matrix(A')

T = eltype(A)
x = randn(T, idim)
y = randn(T, odim)

@assert sum((A*x) .* y) ≈ sum(x .* (A'*y)) # <Ax,y> = <x,A'y>