Conditional Maximum Likelihood for Rasch Model

example.jl

using Optim
using DelimitedFiles  
using ForwardDiff

#data
X = readdlm("lsat.txt");
n, m = size(X)
r = vec(mapslices(sum, X, dims=2))
s = vec(mapslices(sum, X, dims=1))
f = [sum(r.==j) for j in 1:m]

#auxiliary functions needed
function esf_sum!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  @inbounds for col in 1:n
    for r in 1:col
        row = col - r + 1 
        S[row+1] = S[row+1] + x[col] * S[row]
    end
  end
end

function esf_ext!(S::AbstractArray{T,1}, H::AbstractArray{T,3}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  esf_sum!(S, x)
  H[:,:,1] .= zero(T)
  H[:,:,2] .= one(T)

  @inbounds for i in 3:n+1
    for j in 1:n
      H[j,j,i] = S[i-1] - x[j] * H[j,j,i-1]
      for k in j+1:n
        H[k,j,i] = S[i-1] - ((x[j]+x[k])*H[k,j,i-1] + x[j]*x[k]*H[k,j,i-2])
        H[j,k,i] = H[k,j,i]
      end
    end
  end
end


##Newton Raphson
esf_input = ones(Float64, m)
beta = zeros(Float64, m)
last_beta = fill(NaN, m)
S = zeros(Float64, m+1)
H = zeros(Float64, m, m, m+1)

function calculate_common!(x, last_x)
  if x != last_x
    copyto!(last_x, x)
    beta .= x
    esf_input .= exp.(-beta)
    esf_ext!(S, H, esf_input)
  end
end
function logLC(beta)
  calculate_common!(beta, last_beta)
  return -s'log.(esf_input) + f'log.(S[2:end])
end
function g!(storage, beta)
  calculate_common!(beta, last_beta)
  for j in 1:m
    storage[j] = s[j]
    for l in 1:m
      storage[j] -= esf_input[j] * f[l] * (H[j,j,l+1] / S[l+1]) 
    end
  end
end
function h!(storage, beta)
  calculate_common!(beta, last_beta)
  for j in 1:m
    for k in 1:m
      storage[k,j] = 0.0
      for l in 1:m
        if j == k
          storage[j,j] += f[l] * (esf_input[j]*H[j,j,l+1] / S[l+1]) * 
            (1 - esf_input[j]*H[j,j,l+1] / S[l+1])
        elseif k > j
          storage[k,j] += esf_input[j] * esf_input[k] * f[l] * 
            ((H[k,j,l] / S[l+1]) - (H[j,j,l+1] * H[k,k,l+1]) / S[l+1] ^ 2)
        else #k < j
          storage[k,j] += esf_input[j] * esf_input[k] * f[l] * 
            ((H[j,k,l] / S[l+1]) - (H[j,j,l+1] * H[k,k,l+1]) / S[l+1] ^ 2)
        end
      end
    end
  end
end
con_c!(c, x) = (c[1] = sum(x); c)
function con_jacobian!(J, x)
  J[1,:] .= ones(length(x))
end
function con_h!(h, x, λ)
    for i in 1:size(h)[1]
        for j in 1:size(h)[2]
            h[i,j] += (i == j) ? λ[1] : 0.0
        end
    end
end
lx = Float64[]; ux = Float64[]
lc = [0.0]; uc = [0.0]
df = TwiceDifferentiable(logLC, g!, h!, zeros(Float64, m))
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, zeros(Float64,m), IPNewton())



#LBFGS with autodiff
struct DiffCache{T<:AbstractArray, S<:AbstractArray}
  du::T
  dual_du::S
end

function DiffCache(T, size, ::Type{Val{chunk_size}}) where chunk_size
  DiffCache(zeros(T, size...), zeros(ForwardDiff.Dual{nothing,T,chunk_size}, size...))
end

DiffCache(u::AbstractArray) = DiffCache(eltype(u), size(u), 
  Val{ForwardDiff.pickchunksize(length(u))})
get_cache(dc::DiffCache, ::Type{T}) where T<: ForwardDiff.Dual = dc.dual_du
get_cache(dc::DiffCache, T) = dc.du

beta_0 = zeros(Float64, m-1)

#cache for ESF values 
z_cache1 = DiffCache(Float64, (m+1,), Val{ForwardDiff.pickchunksize(m-1)})
z_cache2 = DiffCache(Float64, (m,), Val{ForwardDiff.pickchunksize(m-1)})

function logLC(beta)
  esf_input = get_cache(z_cache2, eltype(beta))
  esf_input_typed = reinterpret(eltype(beta[1]), esf_input)
  esf_input_typed[1] = exp(sum(beta))
  esf_input_typed[2:end] .= exp.(-beta)

  #get correctly typed cache value and fix tag
  esf = get_cache(z_cache1, eltype(beta))
  esf_typed = reinterpret(eltype(beta[1]), esf)

  #in-place ESF 
  esf_sum!(esf_typed, esf_input_typed)
  
  return -s'log.(esf_input_typed) + f'log.(esf_typed[2:end]) 
end

df = OnceDifferentiable(logLC, beta_0; autodiff = :forward)
res = Optim.optimize(df, beta_0)

lsat.txt

0	0	0	0	0
0	0	0	0	0
0	0	0	0	0
0	0	0	0	1
0	0	0	0	1
0	0	0	0	1
0	0	0	0	1
0	0	0	0	1
0	0	0	0	1
0	0	0	1	0
0	0	0	1	0
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	0	1	1
0	0	1	0	0
0	0	1	0	1
0	0	1	1	0
0	0	1	1	0
0	0	1	1	0
0	0	1	1	1
0	0	1	1	1
0	0	1	1	1
0	0	1	1	1
0	1	0	0	0
0	1	0	0	1
0	1	0	0	1
0	1	0	0	1
0	1	0	0	1
0	1	0	0	1
0	1	0	0	1
0	1	0	0	1
0	1	0	0	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	0	1	1
0	1	1	0	1
0	1	1	0	1
0	1	1	0	1
0	1	1	1	0
0	1	1	1	0
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
0	1	1	1	1
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	0	1
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	0
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	0	1	1
1	0	1	0	0
1	0	1	0	0
1	0	1	0	0
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	0	1
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	0
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	0	1	1	1
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	0
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	0	1
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	0
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	0	1	1
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	0
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	0	1
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	0
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1
1	1	1	1	1