text:
function clenshaw_evaluate(weights::Array{T,1},x::Array{T,1},order::Array{S,1},domain::Array{T,1}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain)

  z = zeros(T,(order[1]+1)+2)

  for i in (order[1]+1):-1:1

	  z[i] = weights[i]+2*x1*z[i+1]-z[i+2]

  end

  y = z[1]-x1*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,2},x::Array{T,1},order::Array{S,1},domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])

  p = zeros(T,order[2]+1)

	for j = 1:(order[2]+1)

    z = zeros(T,(order[1]+1)+2)

    for i in (order[1]+1):-1:1

      z[i] = weights[i,j]+2*x1*z[i+1]-z[i+2]

    end

    p[j] = z[1]-x1*z[2]

  end

  z = zeros(T,(order[2]+1)+2)

  for j in (order[2]+1):-1:1

    z[j] = p[j]+2*x2*z[j+1]-z[j+2]

  end

  y = z[1]-x2*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,3},x::Array{T,1},order::Array{S,1},domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])

	pp = zeros(T,(order[2]+1),(order[3]+1))

	for k  = 1:(order[3]+1)
		for j = 1:(order[2]+1)

      z = zeros(T,(order[1]+1)+2)

      for i in (order[1]+1):-1:1

        z[i] = weights[i,j,k]+2*x1*z[i+1]-z[i+2]

      end

      pp[j,k] = z[1]-x1*z[2]

    end
	end

	p = zeros(T,(order[3]+1))

	for k = 1:(order[3]+1)

    z = zeros(T,(order[2]+1)+2)

    for j in (order[2]+1):-1:1

      z[j] = pp[j,k]+2*x2*z[j+1]-z[j+2]

    end

    p[k] = z[1]-x2*z[2]

  end

  z = zeros(T,(order[3]+1)+2)

  for k in (order[3]+1):-1:1

    z[k] = p[k]+2*x3*z[k+1]-z[k+2]

  end

  y = z[1]-x3*z[2]

	return y

end

function clenshaw_evaluate(weights::Array{T,4},x::Array{T,1},order::Array{S,1},domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])

	ppp = zeros(T,(order[2]+1),(order[3]+1),(order[4]+1))

	for l = 1:(order[4]+1)
	  for k = 1:(order[3]+1)
		  for j = 1:(order[2]+1)

        z = zeros(T,(order[1]+1)+2)

        for i in (order[1]+1):-1:1

          z[i] = weights[i,j,k,l]+2*x1*z[i+1]-z[i+2]

        end

        ppp[j,k,l] = z[1]-x1*z[2]

      end
	  end
	end

	pp = zeros(T,(order[3]+1),(order[4]+1))

	for l = 1:(order[4]+1)
		for k = 1:(order[3]+1)

      z = zeros(T,(order[2]+1)+2)

      for j in (order[2]+1):-1:1

        z[j] = ppp[j,k,l]+2*x2*z[j+1]-z[j+2]

      end

      pp[k,l] = z[1]-x2*z[2]

    end
	end

	p = zeros(T,(order[4]+1))

	for l = 1:(order[4]+1)

    z = zeros(T,(order[3]+1)+2)

    for k in (order[3]+1):-1:1

      z[k] = pp[k,l]+2*x3*z[k+1]-z[k+2]

    end

    p[l] = z[1]-x3*z[2]

  end

  z = zeros(T,(order[4]+1)+2)

  for l in (order[4]+1):-1:1

    z[l] = p[l]+2*x4*z[l+1]-z[l+2]

  end

  y = z[1]-x4*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,5},x::Array{T,1},order::Array{S,1},domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])

  pppp = zeros(T,(order[2]+1),(order[3]+1),(order[4]+1),(order[5]+1))

  for m = 1:(order[5]+1)
	  for l = 1:(order[4]+1)
	    for k = 1:(order[3]+1)
		    for j = 1:(order[2]+1)

          z = zeros(T,(order[1]+1)+2)

          for i in (order[1]+1):-1:1

            z[i] = weights[i,j,k,l,m]+2*x1*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x1*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order[3]+1),(order[4]+1),(order[5]+1))

	for m = 1:(order[5]+1)
	  for l = 1:(order[4]+1)
		  for k = 1:(order[3]+1)

        z = zeros(T,(order[2]+1)+2)

        for j in (order[2]+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x2*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x2*z[2]

      end
	  end
	end

	pp = zeros(T,(order[4]+1),(order[5]+1))

	for m = 1:(order[5]+1)
		for l = 1:(order[4]+1)

      z = zeros(T,(order[3]+1)+2)

      for k in (order[3]+1):-1:1

        z[k] = ppp[k,l,m]+2*x3*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x3*z[2]

    end
	end

	p = zeros(T,(order[5]+1))

	for m = 1:(order[5]+1)

    z = zeros(T,(order[4]+1)+2)

    for l in (order[4]+1):-1:1

      z[l] = pp[l,m]+2*x4*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x4*z[2]

  end

  z = zeros(T,(order[5]+1)+2)

  for m in (order[5]+1):-1:1

    z[m] = p[m]+2*x5*z[m+1]-z[m+2]

  end

  y = z[1]-x5*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,6},x::Array{T,1},order::Array{S,1},domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])
  x6 = normalize_node(x[6],domain[:,6])

  ppppp = zeros(T,(order[2]+1),(order[3]+1),(order[4]+1),(order[5]+1),(order[6]+1))

  for n = 1:(order[6]+1)
    for m = 1:(order[5]+1)
	    for l = 1:(order[4]+1)
	      for k = 1:(order[3]+1)
		      for j = 1:(order[2]+1)

            z = zeros(T,(order[1]+1)+2)

            for i in (order[1]+1):-1:1

              z[i] = weights[i,j,k,l,m,n]+2*x1*z[i+1]-z[i+2]

            end

            ppppp[j,k,l,m,n] = z[1]-x1*z[2]

          end
	      end
	    end
    end
  end

  pppp = zeros(T,(order[3]+1),(order[4]+1),(order[5]+1),(order[6]+1))

  for m = 1:(order[6]+1)
	  for l = 1:(order[5]+1)
	    for k = 1:(order[4]+1)
		    for j = 1:(order[3]+1)

          z = zeros(T,(order[2]+1)+2)

          for i in (order[2]+1):-1:1

            z[i] = ppppp[i,j,k,l,m]+2*x1*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x1*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order[4]+1),(order[5]+1),(order[6]+1))

	for m = 1:(order[6]+1)
	  for l = 1:(order[5]+1)
		  for k = 1:(order[4]+1)

        z = zeros(T,(order[3]+1)+2)

        for j in (order[3]+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x2*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x2*z[2]

      end
	  end
	end

	pp = zeros(T,(order[5]+1),(order[6]+1))

	for m = 1:(order[6]+1)
		for l = 1:(order[5]+1)

      z = zeros(T,(order[4]+1)+2)

      for k in (order[4]+1):-1:1

        z[k] = ppp[k,l,m]+2*x3*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x3*z[2]

    end
	end

	p = zeros(T,(order[6]+1))

	for m = 1:(order[6]+1)

    z = zeros(T,(order[5]+1)+2)

    for l in (order[5]+1):-1:1

      z[l] = pp[l,m]+2*x4*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x4*z[2]

  end

  z = zeros(T,(order[6]+1)+2)

  for m in (order[6]+1):-1:1

    z[m] = p[m]+2*x5*z[m+1]-z[m+2]

  end

  y = z[1]-x5*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,1},x::Array{T,1},order::S,domain::Array{T,1}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain)

  z = zeros(T,(order+1)+2)

  for i in (order+1):-1:1

	  z[i] = weights[i]+2*x1*z[i+1]-z[i+2]

  end

  y = z[1]-x1*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,2},x::Array{T,1},order::S,domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])

  p = zeros(T,order+1)

	for j = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for i in (order+1):-1:1

      z[i] = weights[i,j]+2*x1*z[i+1]-z[i+2]

    end

    p[j] = z[1]-x1*z[2]

  end

  z = zeros(T,(order+1)+2)

  for j in (order+1):-1:1

    z[j] = p[j]+2*x2*z[j+1]-z[j+2]

  end

  y = z[1]-x2*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,3},x::Array{T,1},order::S,domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])

	pp = zeros(T,(order+1),(order+1))

	for k  = 1:(order+1)
		for j = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for i in (order+1):-1:1

        z[i] = weights[i,j,k]+2*x1*z[i+1]-z[i+2]

      end

      pp[j,k] = z[1]-x1*z[2]

    end
	end

	p = zeros(T,(order+1))

	for k = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for j in (order+1):-1:1

      z[j] = pp[j,k]+2*x2*z[j+1]-z[j+2]

    end

    p[k] = z[1]-x2*z[2]

  end

  z = zeros(T,(order+1)+2)

  for k in (order+1):-1:1

    z[k] = p[k]+2*x3*z[k+1]-z[k+2]

  end

  y = z[1]-x3*z[2]

	return y

end

function clenshaw_evaluate(weights::Array{T,4},x::Array{T,1},order::S,domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])

	ppp = zeros(T,(order+1),(order+1),(order+1))

	for l = 1:(order+1)
	  for k = 1:(order+1)
		  for j = 1:(order+1)

        z = zeros(T,(order+1)+2)

        for i in (order+1):-1:1

          z[i] = weights[i,j,k,l]+2*x1*z[i+1]-z[i+2]

        end

        ppp[j,k,l] = z[1]-x1*z[2]

      end
	  end
	end

	pp = zeros(T,(order+1),(order+1))

	for l = 1:(order+1)
		for k = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for j in (order+1):-1:1

        z[j] = ppp[j,k,l]+2*x2*z[j+1]-z[j+2]

      end

      pp[k,l] = z[1]-x2*z[2]

    end
	end

	p = zeros(T,(order+1))

	for l = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for k in (order+1):-1:1

      z[k] = pp[k,l]+2*x3*z[k+1]-z[k+2]

    end

    p[l] = z[1]-x3*z[2]

  end

  z = zeros(T,(order+1)+2)

  for l in (order+1):-1:1

    z[l] = p[l]+2*x4*z[l+1]-z[l+2]

  end

  y = z[1]-x4*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,5},x::Array{T,1},order::S,domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])

  pppp = zeros(T,(order+1),(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
	    for k = 1:(order+1)
		    for j = 1:(order+1)

          z = zeros(T,(order+1)+2)

          for i in (order+1):-1:1

            z[i] = weights[i,j,k,l,m]+2*x1*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x1*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
		  for k = 1:(order+1)

        z = zeros(T,(order+1)+2)

        for j in (order+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x2*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x2*z[2]

      end
	  end
	end

	pp = zeros(T,(order+1),(order+1))

	for m = 1:(order+1)
		for l = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for k in (order+1):-1:1

        z[k] = ppp[k,l,m]+2*x3*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x3*z[2]

    end
	end

	p = zeros(T,(order+1))

	for m = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for l in (order+1):-1:1

      z[l] = pp[l,m]+2*x4*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x4*z[2]

  end

  z = zeros(T,(order+1)+2)

  for m in (order+1):-1:1

    z[m] = p[m]+2*x5*z[m+1]-z[m+2]

  end

  y = z[1]-x5*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,6},x::Array{T,1},order::S,domain::Array{T,2}) where {T<:AbstractFloat,S<:Integer}

  x1 = normalize_node(x[1],domain[:,1])
  x2 = normalize_node(x[2],domain[:,2])
  x3 = normalize_node(x[3],domain[:,3])
  x4 = normalize_node(x[4],domain[:,4])
  x5 = normalize_node(x[5],domain[:,5])
  x6 = normalize_node(x[6],domain[:,6])

  ppppp = zeros(T,(order+1),(order+1),(order+1),(order+1),(order+1))

	for n = 1:(order+1)
	  for m = 1:(order+1)
	    for l = 1:(order+1)
	      for k = 1:(order+1)
		      for j = 1:(order+1)

            z = zeros(T,(order+1)+2)

            for i in (order+1):-1:1

              z[i] = weights[i,j,k,l,m,n]+2*x1*z[i+1]-z[i+2]

            end

            ppppp[j,k,l,m,n] = z[1]-x1*z[2]

          end
	      end
	    end
    end
  end

  pppp = zeros(T,(order+1),(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
	    for k = 1:(order+1)
		    for j = 1:(order+1)

          z = zeros(T,(order+1)+2)

          for i in (order+1):-1:1

            z[i] = ppppp[i,j,k,l,m]+2*x1*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x1*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
		  for k = 1:(order+1)

        z = zeros(T,(order+1)+2)

        for j in (order+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x2*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x2*z[2]

      end
	  end
	end

	pp = zeros(T,(order+1),(order+1))

	for m = 1:(order+1)
		for l = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for k in (order+1):-1:1

        z[k] = ppp[k,l,m]+2*x3*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x3*z[2]

    end
	end

	p = zeros(T,(order+1))

	for m = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for l in (order+1):-1:1

      z[l] = pp[l,m]+2*x4*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x4*z[2]

  end

  z = zeros(T,(order+1)+2)

  for m in (order+1):-1:1

    z[m] = p[m]+2*x5*z[m+1]-z[m+2]

  end

  y = z[1]-x5*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,1},x::Array{T,1},order::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

  z = zeros(T,(order[1]+1)+2)

  for i in (order[1]+1):-1:1

	  z[i] = weights[i]+2*x[1]*z[i+1]-z[i+2]

  end

  y = z[1]-x[1]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,2},x::Array{T,1},order::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

  p = zeros(T,order[2]+1)

	for j = 1:(order[2]+1)

    z = zeros(T,(order[1]+1)+2)

    for i in (order[1]+1):-1:1

      z[i] = weights[i,j]+2*x[1]*z[i+1]-z[i+2]

    end

    p[j] = z[1]-x[1]*z[2]

  end

  z = zeros(T,(order[2]+1)+2)

  for j in (order[2]+1):-1:1

    z[j] = p[j]+2*x[2]*z[j+1]-z[j+2]

  end

  y = z[1]-x[2]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,3},x::Array{T,1},order::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

	pp = zeros(T,(order[2]+1),(order[3]+1))

	for k  = 1:(order[3]+1)
		for j = 1:(order[2]+1)

      z = zeros(T,(order[1]+1)+2)

      for i in (order[1]+1):-1:1

        z[i] = weights[i,j,k]+2*x[1]*z[i+1]-z[i+2]

      end

      pp[j,k] = z[1]-x[1]*z[2]

    end
	end

	p = zeros(T,(order[3]+1))

	for k = 1:(order[3]+1)

    z = zeros(T,(order[2]+1)+2)

    for j in (order[2]+1):-1:1

      z[j] = pp[j,k]+2*x[2]*z[j+1]-z[j+2]

    end

    p[k] = z[1]-x[2]*z[2]

  end

  z = zeros(T,(order[3]+1)+2)

  for k in (order[3]+1):-1:1

    z[k] = p[k]+2*x[3]*z[k+1]-z[k+2]

  end

  y = z[1]-x[3]*z[2]

	return y

end

function clenshaw_evaluate(weights::Array{T,4},x::Array{T,1},order::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

	ppp = zeros(T,(order[2]+1),(order[3]+1),(order[4]+1))

	for l = 1:(order[4]+1)
	  for k = 1:(order[3]+1)
		  for j = 1:(order[2]+1)

        z = zeros(T,(order[1]+1)+2)

        for i in (order[1]+1):-1:1

          z[i] = weights[i,j,k,l]+2*x[1]*z[i+1]-z[i+2]

        end

        ppp[j,k,l] = z[1]-x[1]*z[2]

      end
	  end
	end

	pp = zeros(T,(order[3]+1),(order[4]+1))

	for l = 1:(order[4]+1)
		for k = 1:(order[3]+1)

      z = zeros(T,(order[2]+1)+2)

      for j in (order[2]+1):-1:1

        z[j] = ppp[j,k,l]+2*x[2]*z[j+1]-z[j+2]

      end

      pp[k,l] = z[1]-x[2]*z[2]

    end
	end

	p = zeros(T,(order[4]+1))

	for l = 1:(order[4]+1)

    z = zeros(T,(order[3]+1)+2)

    for k in (order[3]+1):-1:1

      z[k] = pp[k,l]+2*x[3]*z[k+1]-z[k+2]

    end

    p[l] = z[1]-x[3]*z[2]

  end

  z = zeros(T,(order[4]+1)+2)

  for l in (order[4]+1):-1:1

    z[l] = p[l]+2*x[4]*z[l+1]-z[l+2]

  end

  y = z[1]-x[4]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,5},x::Array{T,1},order::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

  pppp = zeros(T,(order[2]+1),(order[3]+1),(order[4]+1),(order[5]+1))

  for m = 1:(order[5]+1)
	  for l = 1:(order[4]+1)
	    for k = 1:(order[3]+1)
		    for j = 1:(order[2]+1)

          z = zeros(T,(order[1]+1)+2)

          for i in (order[1]+1):-1:1

            z[i] = weights[i,j,k,l,m]+2*x[1]*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x[1]*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order[3]+1),(order[4]+1),(order[5]+1))

	for m = 1:(order[5]+1)
	  for l = 1:(order[4]+1)
		  for k = 1:(order[3]+1)

        z = zeros(T,(order[2]+1)+2)

        for j in (order[2]+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x[2]*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x[2]*z[2]

      end
	  end
	end

	pp = zeros(T,(order[4]+1),(order[5]+1))

	for m = 1:(order[5]+1)
		for l = 1:(order[4]+1)

      z = zeros(T,(order[3]+1)+2)

      for k in (order[3]+1):-1:1

        z[k] = ppp[k,l,m]+2*x[3]*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x[3]*z[2]

    end
	end

	p = zeros(T,(order[5]+1))

	for m = 1:(order[5]+1)

    z = zeros(T,(order[4]+1)+2)

    for l in (order[4]+1):-1:1

      z[l] = pp[l,m]+2*x[4]*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x[4]*z[2]

  end

  z = zeros(T,(order[5]+1)+2)

  for m in (order[5]+1):-1:1

    z[m] = p[m]+2*x[5]*z[m+1]-z[m+2]

  end

  y = z[1]-x[5]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,6},x::Array{T,1},order::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

  ppppp = zeros(T,(order[2]+1),(order[3]+1),(order[4]+1),(order[5]+1),(order[6]+1))

  for n = 1:(order[6]+1)
    for m = 1:(order[5]+1)
	    for l = 1:(order[4]+1)
	      for k = 1:(order[3]+1)
		      for j = 1:(order[2]+1)

            z = zeros(T,(order[1]+1)+2)

            for i in (order[1]+1):-1:1

              z[i] = weights[i,j,k,l,m,n]+2*x1*z[i+1]-z[i+2]

            end

            ppppp[j,k,l,m,n] = z[1]-x1*z[2]

          end
	      end
	    end
    end
  end

  pppp = zeros(T,(order[3]+1),(order[4]+1),(order[5]+1),(order[6]+1))

  for m = 1:(order[6]+1)
	  for l = 1:(order[5]+1)
	    for k = 1:(order[4]+1)
		    for j = 1:(order[3]+1)

          z = zeros(T,(order[2]+1)+2)

          for i in (order[2]+1):-1:1

            z[i] = ppppp[i,j,k,l,m]+2*x1*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x1*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order[4]+1),(order[5]+1),(order[6]+1))

	for m = 1:(order[6]+1)
	  for l = 1:(order[5]+1)
		  for k = 1:(order[4]+1)

        z = zeros(T,(order[3]+1)+2)

        for j in (order[3]+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x2*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x2*z[2]

      end
	  end
	end

	pp = zeros(T,(order[5]+1),(order[6]+1))

	for m = 1:(order[6]+1)
		for l = 1:(order[5]+1)

      z = zeros(T,(order[4]+1)+2)

      for k in (order[4]+1):-1:1

        z[k] = ppp[k,l,m]+2*x3*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x3*z[2]

    end
	end

	p = zeros(T,(order[6]+1))

	for m = 1:(order[6]+1)

    z = zeros(T,(order[5]+1)+2)

    for l in (order[5]+1):-1:1

      z[l] = pp[l,m]+2*x4*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x4*z[2]

  end

  z = zeros(T,(order[6]+1)+2)

  for m in (order[6]+1):-1:1

    z[m] = p[m]+2*x5*z[m+1]-z[m+2]

  end

  y = z[1]-x5*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,1},x::Array{T,1},order::S) where {T<:AbstractFloat,S<:Integer}

  z = zeros(T,(order+1)+2)

  for i in (order+1):-1:1

	  z[i] = weights[i]+2*x[1]*z[i+1]-z[i+2]

  end

  y = z[1]-x[1]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,2},x::Array{T,1},order::S) where {T<:AbstractFloat,S<:Integer}

  p = zeros(T,order+1)

	for j = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for i in (order+1):-1:1

      z[i] = weights[i,j]+2*x[1]*z[i+1]-z[i+2]

    end

    p[j] = z[1]-x[1]*z[2]

  end

  z = zeros(T,(order+1)+2)

  for j in (order+1):-1:1

    z[j] = p[j]+2*x[2]*z[j+1]-z[j+2]

  end

  y = z[1]-x[2]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,3},x::Array{T,1},order::S) where {T<:AbstractFloat,S<:Integer}

	pp = zeros(T,(order+1),(order+1))

	for k  = 1:(order+1)
		for j = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for i in (order+1):-1:1

        z[i] = weights[i,j,k]+2*x[1]*z[i+1]-z[i+2]

      end

      pp[j,k] = z[1]-x[1]*z[2]

    end
	end

	p = zeros(T,(order+1))

	for k = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for j in (order+1):-1:1

      z[j] = pp[j,k]+2*x[2]*z[j+1]-z[j+2]

    end

    p[k] = z[1]-x[2]*z[2]

  end

  z = zeros(T,(order+1)+2)

  for k in (order+1):-1:1

    z[k] = p[k]+2*x[3]*z[k+1]-z[k+2]

  end

  y = z[1]-x[3]*z[2]

	return y

end

function clenshaw_evaluate(weights::Array{T,4},x::Array{T,1},order::S) where {T<:AbstractFloat,S<:Integer}

	ppp = zeros(T,(order+1),(order+1),(order+1))

	for l = 1:(order+1)
	  for k = 1:(order+1)
		  for j = 1:(order+1)

        z = zeros(T,(order+1)+2)

        for i in (order+1):-1:1

          z[i] = weights[i,j,k,l]+2*x[1]*z[i+1]-z[i+2]

        end

        ppp[j,k,l] = z[1]-x[1]*z[2]

      end
	  end
	end

	pp = zeros(T,(order+1),(order+1))

	for l = 1:(order+1)
		for k = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for j in (order+1):-1:1

        z[j] = ppp[j,k,l]+2*x[2]*z[j+1]-z[j+2]

      end

      pp[k,l] = z[1]-x[2]*z[2]

    end
	end

	p = zeros(T,(order+1))

	for l = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for k in (order+1):-1:1

      z[k] = pp[k,l]+2*x[3]*z[k+1]-z[k+2]

    end

    p[l] = z[1]-x[3]*z[2]

  end

  z = zeros(T,(order+1)+2)

  for l in (order+1):-1:1

    z[l] = p[l]+2*x[4]*z[l+1]-z[l+2]

  end

  y = z[1]-x[4]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,5},x::Array{T,1},order::S) where {T<:AbstractFloat,S<:Integer}

  pppp = zeros(T,(order+1),(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
	    for k = 1:(order+1)
		    for j = 1:(order+1)

          z = zeros(T,(order+1)+2)

          for i in (order+1):-1:1

            z[i] = weights[i,j,k,l,m]+2*x[1]*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x[1]*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
		  for k = 1:(order+1)

        z = zeros(T,(order+1)+2)

        for j in (order+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x[2]*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x[2]*z[2]

      end
	  end
	end

	pp = zeros(T,(order+1),(order+1))

	for m = 1:(order+1)
		for l = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for k in (order+1):-1:1

        z[k] = ppp[k,l,m]+2*x[3]*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x[3]*z[2]

    end
	end

	p = zeros(T,(order+1))

	for m = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for l in (order+1):-1:1

      z[l] = pp[l,m]+2*x[4]*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x[4]*z[2]

  end

  z = zeros(T,(order+1)+2)

  for m in (order+1):-1:1

    z[m] = p[m]+2*x[5]*z[m+1]-z[m+2]

  end

  y = z[1]-x[5]*z[2]

  return y

end

function clenshaw_evaluate(weights::Array{T,6},x::Array{T,1},order::S) where {T<:AbstractFloat,S<:Integer}

  ppppp = zeros(T,(order+1),(order+1),(order+1),(order+1),(order+1))

	for n = 1:(order+1)
	  for m = 1:(order+1)
	    for l = 1:(order+1)
	      for k = 1:(order+1)
		      for j = 1:(order+1)

            z = zeros(T,(order+1)+2)

            for i in (order+1):-1:1

              z[i] = weights[i,j,k,l,m,n]+2*x1*z[i+1]-z[i+2]

            end

            ppppp[j,k,l,m,n] = z[1]-x1*z[2]

          end
	      end
	    end
    end
  end

  pppp = zeros(T,(order+1),(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
	    for k = 1:(order+1)
		    for j = 1:(order+1)

          z = zeros(T,(order+1)+2)

          for i in (order+1):-1:1

            z[i] = ppppp[i,j,k,l,m]+2*x1*z[i+1]-z[i+2]

          end

          pppp[j,k,l,m] = z[1]-x1*z[2]

        end
	    end
	  end
  end

	ppp = zeros(T,(order+1),(order+1),(order+1))

	for m = 1:(order+1)
	  for l = 1:(order+1)
		  for k = 1:(order+1)

        z = zeros(T,(order+1)+2)

        for j in (order+1):-1:1

          z[j] = pppp[j,k,l,m]+2*x2*z[j+1]-z[j+2]

        end

        ppp[k,l,m] = z[1]-x2*z[2]

      end
	  end
	end

	pp = zeros(T,(order+1),(order+1))

	for m = 1:(order+1)
		for l = 1:(order+1)

      z = zeros(T,(order+1)+2)

      for k in (order+1):-1:1

        z[k] = ppp[k,l,m]+2*x3*z[k+1]-z[k+2]

      end

      pp[l,m] = z[1]-x3*z[2]

    end
	end

	p = zeros(T,(order+1))

	for m = 1:(order+1)

    z = zeros(T,(order+1)+2)

    for l in (order+1):-1:1

      z[l] = pp[l,m]+2*x4*z[l+1]-z[l+2]

    end

    p[m] = z[1]-x4*z[2]

  end

  z = zeros(T,(order+1)+2)

  for m in (order+1):-1:1

    z[m] = p[m]+2*x5*z[m+1]-z[m+2]

  end

  y = z[1]-x5*z[2]

  return y

end

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
#=
Examples of operators.
=#

println("6 + 3 = ", 6 + 3)

println("3 - 8 = ", 3 - 8)

println("4 * 5 = ", 4 * 5)

println("14 / 5 = ", 14 / 5)

println("14 % 5 = ", 14 % 5)

println("9^2 = ", 9^2)

ext: jl

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
export CartPoleEnv

struct CartPoleEnvParams{T}
    gravity::T
    masscart::T
    masspole::T
    totalmass::T
    halflength::T
    polemasslength::T
    forcemag::T
    dt::T
    thetathreshold::T
    xthreshold::T
    max_steps::Int
end

Base.show(io::IO, params::CartPoleEnvParams) = print(
    io,
    join(["$p=$(getfield(params, p))" for p in fieldnames(CartPoleEnvParams)], ","),
)

mutable struct CartPoleEnv{T,R<:AbstractRNG} <: AbstractEnv
    params::CartPoleEnvParams{T}
    state::Array{T,1}
    action::Int
    done::Bool
    t::Int
    rng::R
end

"""
    CartPoleEnv(;kwargs...)

# Keyword arguments
- `T = Float64`
- `gravity = T(9.8)`
- `masscart = T(1.0)`
- `masspole = T(0.1)`
- `halflength = T(0.5)`
- `forcemag = T(10.0)`
- `max_steps = 200`
- 'dt = 0.02'
- `rng = Random.GLOBAL_RNG`
"""
function CartPoleEnv(;
    T = Float64,
    gravity = 9.8,
    masscart = 1.0,
    masspole = 0.1,
    halflength = 0.5,
    forcemag = 10.0,
    max_steps = 200,
    dt = 0.02,
    rng = Random.GLOBAL_RNG,
)
    params = CartPoleEnvParams{T}(
        gravity,
        masscart,
        masspole,
        masscart + masspole,
        halflength,
        masspole * halflength,
        forcemag,
        dt,
        2 * 12 * π / 360,
        2.4,
        max_steps,
    )
    high = cp = CartPoleEnv(params, zeros(T, 4), 2, false, 0, rng)
    reset!(cp)
    cp
end

CartPoleEnv{T}(; kwargs...) where {T} = CartPoleEnv(; T = T, kwargs...)

function RLBase.reset!(env::CartPoleEnv{T}) where {T<:Number}
    env.state[:] = T(0.1) * rand(env.rng, T, 4) .- T(0.05)
    env.t = 0
    env.action = 2
    env.done = false
    nothing
end

RLBase.action_space(env::CartPoleEnv) = Base.OneTo(2)

RLBase.state_space(env::CartPoleEnv{T}) where {T} = Space(
    ClosedInterval{T}[
        (-2*env.params.xthreshold)..(2*env.params.xthreshold),
        -1e38..1e38,
        (-2*env.params.thetathreshold)..(2*env.params.thetathreshold),
        -1e38..1e38,
    ],
)

RLBase.reward(env::CartPoleEnv{T}) where {T} = env.done ? zero(T) : one(T)
RLBase.is_terminated(env::CartPoleEnv) = env.done
RLBase.state(env::CartPoleEnv) = env.state

function (env::CartPoleEnv)(a)
    @assert a in (1, 2)
    env.action = a
    env.t += 1
    force = a == 2 ? env.params.forcemag : -env.params.forcemag
    x, xdot, theta, thetadot = env.state
    costheta = cos(theta)
    sintheta = sin(theta)
    tmp = (force + env.params.polemasslength * thetadot^2 * sintheta) / env.params.totalmass
    thetaacc =
        (env.params.gravity * sintheta - costheta * tmp) / (
            env.params.halflength *
            (4 / 3 - env.params.masspole * costheta^2 / env.params.totalmass)
        )
    xacc = tmp - env.params.polemasslength * thetaacc * costheta / env.params.totalmass
    env.state[1] += env.params.dt * xdot
    env.state[2] += env.params.dt * xacc
    env.state[3] += env.params.dt * thetadot
    env.state[4] += env.params.dt * thetaacc
    env.done =
        abs(env.state[1]) > env.params.xthreshold ||
        abs(env.state[3]) > env.params.thetathreshold ||
        env.t > env.params.max_steps
    nothing
end

Random.seed!(env::CartPoleEnv, seed) = Random.seed!(env.rng, seed)

function plotendofepisode(x, y, d)
    if d
        setmarkercolorind(7)
        setmarkertype(-1)
        setmarkersize(6)
        polymarker([x], [y])
    end
    return nothing
end

function GR.plot(env::CartPoleEnv)
    s, a, d = env.state, env.action, env.done
    x, xdot, theta, thetadot = s
    l = 2 * env.params.halflength
    clearws()
    setviewport(0, 1, 0, 1)
    xthreshold = env.params.xthreshold
    setwindow(-xthreshold, xthreshold, -.1, l + 0.1)
    fillarea([x - 0.5, x - 0.5, x + 0.5, x + 0.5], [-.05, 0, 0, -.05])
    setlinecolorind(4)
    setlinewidth(3)
    polyline([x, x + l * sin(theta)], [0, l * cos(theta)])
    setlinecolorind(2)
    drawarrow(x + (a == 1) - 0.5, -.025, x + 1.4 * (a == 1) - 0.7, -.025)
    plotendofepisode(xthreshold - 0.2, l, d)
    updatews()
end


text:
export TR_Nwt_ls

"""
A one dimensionnal trust region algorithm using the exact second derivative in
the second order Taylor expansion. For more documentation referer to TR_generic_ls.
"""
function TR_Nwt_ls(h :: LineModel,  stop_ls :: LS_Stopping;  kwargs...)

    (state, stop_ls.meta.optimal) = TR_generic_ls(h, stop_ls; kwargs...)
    return (state, stop_ls.meta.optimal)
end


text:
"""
2. Velocity and pressure only (initialize with equilibrium)
f₀(x) = f_eq(ρ₀(x), u₀(x), p₀(x))
ρ₀ = ρ₀(x)
u₀ = u₀(x)
p₀ = p₀(x)
"""
struct AnalyticalEquilibrium <: InitializationStrategy end
function initial_condition(
    ::AnalyticalEquilibrium,
    q::Quadrature,
    problem::FluidFlowProblem,
    x::T,
    y::T,
) where {T <: Real}
    return equilibrium(q, problem, x, y)
end


text:
using GLM, DataFrames, Distributions, Plots, CSV; pyplot()

data = CSV.read("../data/examData.csv")
model = glm(@formula(Pass ~ Hours), data, Binomial(), LogitLink())
println(model)

b0, b1 = coef(model)
pred(x) = 1/(1+exp(-(b0 + b1*x)))

xGrid = 0:0.1:maximum(data.Hours)
scatter(data.Hours, data.Pass, c=:blue, msw=0)
plot!(xGrid, pred.(xGrid),
        c=:red, xlabel="Hours studied", legend=:none,
        xlims=(0, maximum(data.Hours)), yticks=([0:1;], ["Fail", "Pass"]))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Reexport

@reexport module TransferOperators

    # Abstract types and composite types for different estimators.
    include("AbstractTransferOperator.jl")

    include("composite_types/TransferOperatorRectangularBinning.jl")
    include("composite_types/TransferOperatorTriangulationExact.jl")
    include("composite_types/TransferOperatorTriangulationApproximate.jl")
    include("composite_types/plot_recipes.jl")

    # Constructors for the estimators.
    include("Estimators/Estimators.jl")

end # module

"""
    TransferOperators

A module containing types and functions to approximate the transfer operator
from numerical data. The estimators in this module are the basis for computing
invariant measures, which is done in the `InvariantMeasures` module.
"""
TransferOperators


text:
using Luxor
using Colors
using BoundaryValueDiffEq

# constants for differential equations and movie
const g = 9.81
const L = 1.0                         # pendulum length in meters
const bobd = 0.10                     # pendulum bob diameter in meters
const framerate = 50.0                # intended frame rate/sec
const t0 = 0.0                        # start time (s)
const tf = 2.3                        # end simulation time (s)
const dtframe = 1.0/framerate         # time increment per frame
const tspan = LinRange(t0, tf, Int(floor(tf*framerate)))  # array of time points in animation

const bgcolor = "black"               # gif background
const leaderhue = (0.80, 0.70, 0.20)  # gif swing arm hue light gold
const hslcolors = [HSL(col) for col in (distinguishable_colors(
                   Int(floor(tf*framerate)+3),[RGB(1,1,1)])[2:end])]
const giffilename = "pendulum.gif"    # output file

# differential equations
simplependulum(du, u, p, t) = (θ=u[1]; dθ=u[2]; du[1]=dθ; du[2]=-(g/L)*sin(θ))
bc2(residual, u, p, t) = (residual[1] = u[end÷2][1] + pi/2; residual[2] = u[end][1] - pi/2)
bvp2 = BVProblem(simplependulum, bc2, [pi/2,pi/2], (tspan[1],tspan[end]))
sol2 = solve(bvp2, MIRK4(), dt=dtframe) # use the MIRK4 solver for TwoPointBVProblem

# movie making background
backdrop(scene, framenumber) = background(bgcolor)

function frame(scene, framenumber)
    u1, u2 = sol2.u[framenumber]
    y, x = L*cos(u1), L*sin(u1)
    sethue(leaderhue)
    poly([Point(-4.0, 0.0), Point(4.0, 0.0),
          Point(160.0x,160.0y)], :fill)
    sethue(Colors.HSV(framenumber*4.0, 1, 1))
    circle(Point(160.0x,160.0y), 160bobd, :fill)
    text(string("frame $framenumber of $(scene.framerange.stop)"),
        Point(0.0, -190.0),
        halign=:center)
end

muv = Movie(400, 400, "Pendulum Demo", 1:length(tspan))
animate(muv, [Scene(muv, backdrop),
              Scene(muv, frame, easingfunction=easeinoutcubic)],
              creategif=true, pathname=giffilename)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
```
moment_tables(m; percent = 0.90, subset_inds = 1:0, subset_string = "",
    verbose = :none)
```

Computes prior and posterior parameter moments. Tabulates prior mean, posterior mean, and
bands in various LaTeX tables stored `tablespath(m)`.

### Inputs

- `m::AbstractModel`: model object

### Keyword Arguments

- `percent::AbstractFloat`: the percentage of the mass of draws from
  Metropolis-Hastings included between the bands displayed in output tables.
- `subset_inds::Range{Int64}`: indices specifying the draws we want to use
- `subset_string::String`: short string identifying the subset to be
  appended to the output filenames. If `subset_inds` is nonempty but
  `subset_string` is empty, an error is thrown
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`
- `use_mode::Bool`: Return a table with the modal parameters as opposed to the mean
"""
function moment_tables(m::AbstractModel; percent::AbstractFloat = 0.90,
                       subset_inds::Range{Int64} = 1:0, subset_string::String = "",
                       groupings::Associative{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}(),
                       verbose::Symbol = :low, use_mode::Bool = false)

    ### 1. Load parameter draws from Metropolis-Hastings

    params = if !isempty(subset_inds)
        # Use subset of draws
        if isempty(subset_string)
            error("Must supply a nonempty subset_string if subset_inds is nonempty")
        end
        load_draws(m, :subset; subset_inds = subset_inds, verbose = verbose)
    else
        # Use all draws
        load_draws(m, :full; verbose = verbose)
    end

    ### 2. Compute posterior moments

    if use_mode
        post_mode = h5read(get_forecast_input_file(m, :mode), "params")
    end
    post_means = vec(mean(params, 1))

    # Save posterior means
    basename = "paramsmean"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    filename = workpath(m, "estimate", "$basename.h5")
    h5open(filename, "w") do file
        write(file, "post_means", post_means)
    end
    post_bands = find_density_bands(params, percent; minimize = true)'

    ### 3. Produce TeX tables

    if use_mode
        prior_posterior_table(m, post_mode; subset_string = subset_string, groupings = groupings, use_mode = use_mode)
    else
        prior_posterior_table(m, post_means; subset_string = subset_string, groupings = groupings, use_mode = use_mode)
    end
    prior_posterior_moments_table(m, post_means, post_bands;
                                  percent = percent, subset_string = subset_string, groupings = groupings)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        @printf "Tables are saved as %s.\n" tablespath(m, "estimate", "*.tex")
    end
end

"""
```
moments(θ::Parameter)
```

If θ's prior is a `DSGE.RootInverseGamma`, τ and ν. Otherwise, returns the mean
and standard deviation of the prior. If θ is fixed, returns `(θ.value, 0.0)`.
"""
function moments(θ::Parameter)
    if θ.fixed
        return θ.value, 0.0
    else
        prior = get(θ.prior)
        if isa(prior, RootInverseGamma)
            return prior.τ, prior.ν
        else
            return mean(prior), std(prior)
        end
    end
end

"""
```
prior_table(m; subset_string = "", groupings = Dict{String, Vector{Parameter}}())
```
"""
function prior_table(m::AbstractModel; subset_string::String = "",
             groupings::Associative{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}())

    if isempty(groupings)
        sorted_parameters = sort(m.parameters, by = (x -> x.key))
        groupings[""] = sorted_parameters
    end

    # Open the TeX file
    basename = "priors"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    priors_out = tablespath(m, "estimate", "$basename.tex")
    priors_fid = open(priors_out, "w")

    # Write header
    write_table_preamble(priors_fid)

    @printf priors_fid "\\renewcommand*\\footnoterule{}"
    @printf priors_fid "\\vspace*{.5cm}\n"
    @printf priors_fid "{\\small\n"
    @printf priors_fid "\\begin{longtable}{rlrr@{\\hspace{1in}}rlrr}\n"
    @printf priors_fid "\\caption{Priors}\n"
    @printf priors_fid "\\\\ \\hline\n"

    @printf priors_fid "& Dist & Mean & Std Dev & & Dist & Mean & Std Dev \\\\ \\hline\n"
    @printf priors_fid "\\endhead\n"

    @printf priors_fid "\\hline \\\\\n"
    @printf priors_fid "\\multicolumn{8}{c}{\\footnotesize Note: For Inverse Gamma prior mean and SD, \$\\tau\$ and \$\\nu\$ reported.}\n"
    @printf priors_fid "\\endfoot\n"

    # Map prior distributions to identifying strings
    distid(::Distributions.Beta)    = "Beta"
    distid(::Distributions.Gamma)   = "Gamma"
    distid(::Distributions.Normal)  = "Normal"
    distid(::DSGE.RootInverseGamma) = "InvG"

    # Write priors
    for group_desc in keys(groupings)
        params = groupings[group_desc]
        n_params = length(params)
        n_rows = convert(Int, ceil(n_params/2))

        # Write grouping description if not empty
        if !isempty(group_desc)
            @printf priors_fid "\\multicolumn{8}{l}{\\textit{%s}} \\\\[3pt]\n" group_desc
        end

        # Write footnote about standard deviations of anticipated policy shocks
        function anticipated_shock_footnote(θ::Parameter)
            if n_anticipated_shocks(m) > 0 && θ.key == :σ_r_m1
                nantpad          = n_anticipated_shocks_padding(m)
                all_sigmas       = [m[Symbol("σ_r_m$i")]::Parameter for i = 1:nantpad]
                nonzero_sigmas   = Base.filter(θ -> !(θ.fixed && θ.value == 0), all_sigmas)
                n_nonzero_sigmas = length(nonzero_sigmas)

                if n_nonzero_sigmas > 1
                    text = "\$\\sigma_{ant1}\$ through \$\\sigma_{ant$(n_nonzero_sigmas)}\$ all have the same distribution."
                    @printf priors_fid "\\let\\thefootnote\\relax\\footnote{\\centering %s}" text
                end
            end
        end

        for i = 1:n_rows
            # Write left column
            θ = params[i]
            (prior_mean, prior_std) = moments(θ)
            @printf priors_fid "\$%s\$ & " θ.tex_label
            @printf priors_fid "%s & " (θ.fixed ? "-" : distid(get(θ.prior)))
            @printf priors_fid "%0.2f & " prior_mean
            @printf priors_fid "%0.2f & " prior_std
            anticipated_shock_footnote(θ)

            # Write right column if it exists
            if n_rows + i <= n_params
                θ = params[n_rows + i]
                (prior_mean, prior_std) = moments(θ)
                @printf priors_fid "\$%s\$ & " θ.tex_label
                @printf priors_fid "%s & " (θ.fixed ? "-" : distid(get(θ.prior)))
                @printf priors_fid "%0.2f & " prior_mean
                @printf priors_fid "%0.2f" prior_std
                anticipated_shock_footnote(θ)
            else
                @printf priors_fid "& & &"
            end

            # Add padding after last row in a grouping
            if i == n_rows
                @printf priors_fid " \\\\[3pt]\n"
            else
                @printf priors_fid " \\\\\n"
            end
        end
    end

    # Write footer
    write_table_postamble(priors_fid; small = true)

    # Close file
    close(priors_fid)
end

"""
```
prior_posterior_moments_table(m, post_means, post_bands; percent = 0.9,
    subset_string = "")
```

Produces a table of prior means, prior standard deviations, posterior means, and
90% bands for posterior draws. Saves to:

```
tablespath(m, \"estimate\", \"moments[_sub=\$subset_string].tex\")
```
"""
function prior_posterior_moments_table(m::AbstractModel,
                 post_means::Vector, post_bands::Matrix;
                 percent::AbstractFloat = 0.9,
                 subset_string::String = "",
                 groupings::Associative{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}())

    if isempty(groupings)
        sorted_parameters = sort(m.parameters, by = (x -> x.key))
        groupings[""] = sorted_parameters
    end

    # Open the TeX file
    basename = "moments"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    moments_out = tablespath(m, "estimate", "$basename.tex")
    moments_fid = open(moments_out, "w")

    # Write header
    write_table_preamble(moments_fid)

    @printf moments_fid "\\vspace*{.5cm}\n"
    @printf moments_fid "{\\small\n"
    @printf moments_fid "\\begin{longtable}{lcccccc}\n"
    @printf moments_fid "\\caption{Parameter Estimates}\n"
    @printf moments_fid "\\\\ \\hline\n"

    # Two-row column names. First row is multicolumn, where entries (i,str) are `i` columns
    # with content `str`.
    colnames0 = [(1,""), (3,"Prior"),(3,"Posterior")]
    colnames = ["Parameter", "Type", "Mean", "SD", "Mean",
                "$(100*percent)\\% {\\tiny Lower Band}",
                "$(100*percent)\\% {\\tiny Upper Band}"]

    @printf moments_fid "\\multicolumn{%d}{c}{%s}" colnames0[1][1] colnames0[1][2]
    for col in colnames0[2:end]
        @printf moments_fid " & \\multicolumn{%d}{c}{%s}" col[1] col[2]
    end
    @printf moments_fid " \\\\\n"
    @printf moments_fid "%s" colnames[1]
    for col in colnames[2:end]
        @printf moments_fid " & %s" col
    end
    @printf moments_fid " \\\\\n"
    @printf moments_fid "\\cmidrule(lr){1-1} \\cmidrule(lr){2-4} \\cmidrule(lr){5-7}\n"
    @printf moments_fid "\\endhead\n"

    @printf moments_fid "\\hline\n"
    @printf moments_fid "\\\\ \\multicolumn{7}{c}{\\footnotesize Note: For Inverse Gamma (IG) prior mean and SD, \$\\tau\$ and \$\\nu\$ reported.}\n"
    @printf moments_fid "\\endfoot\n"

    # Map prior distributions to identifying strings
    distid(::Distributions.Beta)    = "B"
    distid(::Distributions.Gamma)   = "G"
    distid(::Distributions.Normal)  = "N"
    distid(::DSGE.RootInverseGamma) = "IG"

    # Write parameter moments
    # sorted_parameters = sort(m.parameters, by = (x -> x.key))
    for group_desc in keys(groupings)
        params = groupings[group_desc]

        # Write grouping description if not empty
        if !isempty(group_desc)
            @printf moments_fid "\\multicolumn{7}{c}{\\textit{%s}} \\\\[3pt]\n" group_desc
        end

        for param in params
            index = m.keys[param.key]
            (prior_mean, prior_std) = moments(param)

            @printf moments_fid "\$\%4.99s\$ & " param.tex_label
            @printf moments_fid "%s & " (param.fixed ? "-" : distid(get(param.prior)))
            @printf moments_fid "%8.3f & " prior_mean
            @printf moments_fid "%8.3f & " prior_std
            @printf moments_fid "%8.3f & " post_means[index]
            @printf moments_fid "%8.3f & %8.3f \\\\\n" post_bands[index, :]...
        end
    end

    # Write footer
    write_table_postamble(moments_fid; small=true)

    # Close file
    close(moments_fid)
end

"""
```
prior_posterior_table(m, post_values; subset_string = "")
```

Produces a table of prior means and posterior means or mode. Saves to:

```
tablespath(m, \"estimate\", \"prior_posterior_means[_sub=\$subset_string].tex\")
```
or
```
tablespath(m, \"estimate\", \"prior_posterior_mode[_sub=\$subset_string].tex\")
```
"""
function prior_posterior_table(m::AbstractModel, post_values::Vector;
                 subset_string::String = "",
                 groupings::Associative{String, Vector{Parameter}} = Dict{String, Vector{Parameter}}(),
                 use_mode::Bool = false)

    if isempty(groupings)
        sorted_parameters = sort(m.parameters, by = (x -> x.key))
        groupings[""] = sorted_parameters
    end

    # Open the TeX file
    basename = use_mode ? "prior_posterior_mode" : "prior_posterior_means"
    if !isempty(subset_string)
        basename *= "_sub=$(subset_string)"
    end
    table_out = tablespath(m, "estimate", "$basename.tex")
    fid = open(table_out, "w")

    # Write header
    write_table_preamble(fid)
    @printf fid "\\vspace*{.5cm}\n"
    @printf fid "\\begin{longtable}{ccc}\n"
    if use_mode
        @printf fid "\\caption{Parameter Estimates: Prior Mean and Posterior Mode}\n"
    else
        @printf fid "\\caption{Parameter Estimates: Prior and Posterior Means}\n"
    end
    @printf fid "\\\\ \\hline\n"
    @printf fid "Parameter & Prior & Posterior\n"
    @printf fid "\\\\ \\hline\n"
    @printf fid "\\endhead\n"
    @printf fid "\\hline\n"
    @printf fid "\\endfoot\n"

    # Write results
    # sorted_parameters = sort(m.parameters, by = (x -> x.key))

    for group_desc in keys(groupings)
        params = groupings[group_desc]

        # Write grouping description if not empty
        if !isempty(group_desc)
            @printf means_fid "\\multicolumn{7}{c}{\\textit{%s}} \\\\[3pt]\n" group_desc
        end

        for param in params
            index = m.keys[param.key]

            post_value = if param.fixed
                param.value
            else
                prior = get(param.prior)
                isa(prior, RootInverseGamma) ? prior.τ : mean(prior)
            end

            @printf fid "\$\%4.99s\$ & " param.tex_label
            @printf fid "%8.3f & " post_value
            @printf fid "\%8.3f \\\\\n" post_values[index]
        end
    end

    # Write footer
    write_table_postamble(fid)

    # Close file
    close(fid)
end

"""
```
find_density_bands(draws::Matrix, percent::AbstractFloat; minimize::Bool=true)
```

Returns a `2` x `cols(draws)` matrix `bands` such that `percent` of the mass of `draws[:,i]`
is above `bands[1,i]` and below `bands[2,i]`.

### Arguments

- `draws`: `ndraws` by `nperiods` matrix of parameter draws (from Metropolis-Hastings, for example)
- `percent`: percent of data within bands (e.g. .9 to get 90% of mass within bands)

### Optional Arguments

- `minimize`: if `true`, choose shortest interval, otherwise just chop off lowest and
  highest (percent/2)
"""
function find_density_bands{T<:AbstractFloat}(draws::Matrix, percent::T; minimize::Bool = true)

    if !(0 <= percent <= 1)
        error("percent must be between 0 and 1")
    end

    ndraws, nperiods = size(draws)

    if ndraws == 1
        band = repmat(draws, 2, 1)
        return band
    end

    band = zeros(2, nperiods)
    n_in_band  = round(Int, percent * ndraws)  # number of draws in the band

    for i in 1:nperiods

        # Sort response for parameter i such that 1st element is largest
        draw_variable_i = draws[:,i]
        sort!(draw_variable_i, rev=true)

        # Search to find the shortest interval containing `percent` of
        # the mass `low` is the index of the largest draw in the band
        # (but the first index to take in `draw_variable_i`, `high` is
        # the smallest (but the highest index to take)

        low = if minimize

            low         = 1
            done        = 0
            j           = 2
            minwidth = draw_variable_i[1] - draw_variable_i[n_in_band]

            while j <= (ndraws - n_in_band + 1)

                newwidth = draw_variable_i[j] - draw_variable_i[j + n_in_band - 1]

                if newwidth < minwidth
                    low = j
                    minwidth = newwidth
                end

                j += 1
            end

            low
        else
            # Chop off lowest and highest percent/2
            ndraws - n_in_band - round(Int, floor(.5*(ndraws-n_in_band)))
        end

        high = low + n_in_band - 1

        band[2,i] = draw_variable_i[low]
        band[1,i] = draw_variable_i[high]
    end

    return band
end

"""
```
find_density_bands{T<:AbstractFloat}(draws::Matrix, percents::Vector{T}; minimize::Bool=true)
```

Returns a `2` x `cols(draws)` matrix `bands` such that `percent` of the mass of `draws[:,i]`
is above `bands[1,i]` and below `bands[2,i]`.

### Arguments

- `draws`: Matrix of parameter draws (from Metropolis-Hastings, for example)
- `percent`: percent of data within bands (e.g. .9 to get 90% of mass within bands)

### Optional Arguments

- `minimize`: if `true`, choose shortest interval, otherwise just chop off lowest and
  highest (percent/2)
"""
function find_density_bands{T<:AbstractFloat}(draws::Matrix, percents::Vector{T}; minimize::Bool = true)

    bands = DataFrame()

    for p in percents
        out = find_density_bands(draws, p, minimize = minimize)

        bands[Symbol("$(100*p)\% UB")] = vec(out[2,:])
        bands[Symbol("$(100*p)\% LB")] = vec(out[1,:])
    end

    bands
end

function write_table_preamble(fid::IOStream)
    @printf fid "\\documentclass[12pt]{article}\n"
    @printf fid "\\usepackage{booktabs}\n"
    @printf fid "\\usepackage[justification=centering]{caption}\n"
    @printf fid "\\usepackage[margin=1in]{geometry}\n"
    @printf fid "\\usepackage{longtable}\n"
    @printf fid "\\begin{document}\n"
    @printf fid "\\pagestyle{empty}\n"
end

# `small`: Whether to print an additional curly bracket after "\end{longtable}" (necessary if
# the table is enclosed by "\small{}")
function write_table_postamble(fid::IOStream; small::Bool=false)
    if small
        @printf fid "\\end{longtable}}\n"
    else
        @printf fid "\\end{longtable}\n"
    end

    @printf fid "\\end{document}"
end


text:
export variable

function convertconstant end
Base.convert(::Type{P}, α) where P<:APL = convertconstant(P, α)
function Base.convert(::Type{P}, p::P) where P<:AbstractPolynomial
    return p
end
function Base.convert(::Type{P}, p::APL) where {T, P<:AbstractPolynomial{T}}
    return convert(P, polynomial(p, T))
end

Base.convert(::Type{Any}, p::APL) = p
# Conversion polynomial -> scalar
function Base.convert(S::Type{<:Union{Number, T}}, p::APL{T}) where T
    s = zero(S)
    for t in terms(p)
        if !isconstant(t)
            # The polynomial is not constant
            throw(InexactError(:convert, S, p))
        end
        s += S(coefficient(t))
    end
    s
end

# Fix ambiguity caused by above conversions
Base.convert(::Type{P}, p::APL) where P<:APL = P(p)

Base.convert(::Type{PT}, p::PT) where {PT<:APL} = p
function Base.convert(::Type{MT}, t::AbstractTerm) where {MT<:AbstractMonomial}
    if isone(coefficient(t))
        monomial(t)
    else
        error("Cannot convert a term with a coefficient that is not one into a monomial")
    end
end


text:
module CombinatorialBandits
  using LinearAlgebra

  using DataStructures
  using IterTools
  using Distributions

  import Munkres # Avoid clashing with Hungarian.
  import Hungarian
  using LightGraphs
  using JuMP

  import Base: push!, copy, hash, isequal
  import JuMP: value

  export Policy, CombinatorialInstance, State,
         initial_state, initial_trace, simulate, choose_action, pull, update!, solve_linear, solve_budgeted_linear, solve_all_budgeted_linear, get_lp_formulation, has_lp_formulation, is_feasible, is_partially_acceptable,
         ThompsonSampling, LLR, CUCB, ESCB2, OLSUCB,
         ThompsonSamplingDetails, LLRDetails, CUCBDetails, ESCB2Details, OLSUCBDetails,
         ESCB2OptimisationAlgorithm, ESCB2Exact, ESCB2Greedy, ESCB2Budgeted, OLSUCBOptimisationAlgorithm, OLSUCBGreedy,
         PerfectBipartiteMatching, UncorrelatedPerfectBipartiteMatching, CorrelatedPerfectBipartiteMatching, PerfectBipartiteMatchingSolver, PerfectBipartiteMatchingLPSolver, PerfectBipartiteMatchingMunkresSolver, PerfectBipartiteMatchingHungarianSolver,
         ElementaryPath, ElementaryPathSolver, ElementaryPathLightGraphsDijkstraSolver, ElementaryPathLPSolver, ElementaryPathAlgosSolver,
         SpanningTree, SpanningTreeSolver, SpanningTreeLightGraphsPrimSolver, SpanningTreeAlgosSolver, SpanningTreeLPSolver,
         MSet, MSetSolver, MSetAlgosSolver, MSetLPSolver,
         # Algos.
         MSetInstance, MSetSolution, dimension, m, value, values, msets_greedy, msets_dp, msets_lp,
         BudgetedMSetInstance, BudgetedMSetSolution, weight, weights, budget, max_weight, items, items_all_budgets, budgeted_msets_dp, budgeted_msets_lp, budgeted_msets_lp_select, budgeted_msets_lp_all,
         ElementaryPathInstance, ElementaryPathSolution, graph, costs, src, dst, cost, lp_dp,
         BudgetedLongestPathInstance, BudgetedLongestPathSolution, rewards, reward, budgeted_lp_dp,
         SpanningTreeInstance, SpanningTreeSolution, st_prim,
         BudgetedSpanningTreeInstance, BudgetedSpanningTreeSolution, BudgetedSpanningTreeLagrangianSolution, SimpleBudgetedSpanningTreeSolution, _budgeted_spanning_tree_compute_value, _budgeted_spanning_tree_compute_weight, st_prim_budgeted_lagrangian, st_prim_budgeted_lagrangian_search, _solution_symmetric_difference, _solution_symmetric_difference_size, st_prim_budgeted_lagrangian_refinement

  # General algorithm.
  abstract type Policy end
  abstract type PolicyDetails end
  abstract type CombinatorialInstance{T} end

  # Define the state of a bandit (evolves at each round).
  mutable struct State{T}
    round::Int
    regret::Float64
    reward::Float64
    arm_counts::Dict{T, Int}
    arm_reward::Dict{T, Float64}
    arm_average_reward::Dict{T, Float64}
  end

  copy(s::State{T}) where T = State{T}(s.round, s.regret, s.reward, s.arm_counts, s.arm_reward, s.arm_average_reward)

  # Define the trace of the execution throughout the rounds.
  struct Trace{T}
    states::Vector{State{T}}
    arms::Vector{Vector{T}}
    reward::Vector{Vector{Float64}}
    policy_details::Vector{PolicyDetails}
    time_choose_action::Vector{Int}
  end

  """
      push!(trace::Trace{T}, state::State{T}, arms::Vector{T}, reward::Vector{Float64}, policy_details::PolicyDetails, time_choose_action::Int) where T

  Appends the arguments to the execution trace of the bandit algorithm. More specifically, `trace`'s data structures are
  updated to also include `state`, `arms`, `reward`, `policy_details`, and `time_choose_action` (expressed in milliseconds).
  All of these arguments are copied, *except* `policy_details`.
  (Indeed, the usual scenario is to keep updating the state, the arms and the rewards, but to build the details at each round from the ground up.)
  """
  function push!(trace::Trace{T}, state::State{T}, arms::Vector{T}, reward::Vector{Float64}, policy_details::PolicyDetails, time_choose_action::Int) where T
    push!(trace.states, copy(state))
    push!(trace.arms, copy(arms))
    push!(trace.reward, copy(reward))
    push!(trace.policy_details, policy_details)
    push!(trace.time_choose_action, time_choose_action)
  end

  # Interface for combinatorial instances.
  function initial_state(instance::CombinatorialInstance{T}) where T end # TODO: Can't this be provided by default based on template types?
  function initial_trace(instance::CombinatorialInstance{T}) where T end # TODO: Can't this be provided by default based on template types?
  function is_feasible(instance::CombinatorialInstance{T}, arms::Vector{T}) where T end
  function is_partially_acceptable(instance::CombinatorialInstance{T}, arms::Vector{T}) where T end

  function all_arm_indices(reward::Matrix{Distribution})
    reward_indices_cartesian = eachindex(view(reward, [1:s for s in size(reward)]...))
    return [Tuple(i) for i in reward_indices_cartesian]
  end

  function all_arm_indices(reward::Vector{Distribution})
    return eachindex(view(reward, [1:s for s in size(reward)]...))
  end

  function all_arm_indices(reward::Dict{T, Distribution}) where T
    return collect(keys(reward))
  end

  function all_arm_indices(instance::CombinatorialInstance{T}) where T
    if isa(instance.reward, MultivariateDistribution) # Correlated arms.
      # Nothing as generic as the uncorrelated case is available, due to
      # the fact that only a vector is known, for any kind of arms (unlike
      # the uncorrelated case, where there is a distinction between vectors
      # and matrices of reward distributions).
      return instance.all_arm_indices
    else # Uncorrelated arms.
      return all_arm_indices(instance.reward)
    end
  end

  function pull(instance::CombinatorialInstance{T}, arms::Vector{T}) where T
    # Draw the rewards for this round. If T is a tuple, the reward distributions
    # are stored in a matrix, hence the splatting.
    arm_indices = all_arm_indices(instance)
    if isa(instance.reward, MultivariateDistribution) # Correlated arms.
      true_rewards_vector = mean(instance.reward)
      true_rewards = Dict{T, Float64}(arm => true_rewards_vector[i] for (i, arm) in enumerate(arm_indices))

      drawn_rewards_vector = rand(instance.reward)
      drawn_rewards = Dict{T, Float64}(arm => true_rewards_vector[i] for (i, arm) in enumerate(arm_indices))
    else # Uncorrelated arms.
      true_rewards = Dict{T, Float64}(i => mean(instance.reward[i...]) for i in arm_indices)
      drawn_rewards = Dict{T, Float64}(i => rand(instance.reward[i...]) for i in arm_indices)
    end

    # Select the information that will be provided back to the bandit policy.
    # Here is implemented the semi-bandit setting.
    true_reward = sum(true_rewards[arm] for arm in arms)
    reward = Float64[drawn_rewards[arm] for arm in arms]

    # Compute the incurred regret from the provided solution.
    incurred_regret = instance.optimal_average_reward - sum(true_reward)

    return reward, incurred_regret
  end

  solve_linear(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}) where T = solve_linear(instance.solver, rewards)
  solve_budgeted_linear(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}, weights::Dict{T, Int}, budget::Int) where T =
    solve_budgeted_linear(instance.solver, rewards, weights, budget)
  solve_all_budgeted_linear(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}, weights::Dict{T, Int}, max_budget::Int) where T =
    solve_all_budgeted_linear(instance.solver, rewards, weights, max_budget)
  has_lp_formulation(instance::CombinatorialInstance{T}) where T = has_lp_formulation(instance.solver)
  get_lp_formulation(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}) where T = has_lp_formulation(instance) ?
    get_lp_formulation(instance.solver, rewards) :
    error("The chosen solver uses no LP formulation.")

  # Implement the most common case. For tuples, the number vary more widely.
  function initial_state(instance::CombinatorialInstance{Int})
    n = instance.n_arms
    zero_counts = Dict(i => 0 for i in 1:n)
    zero_rewards = Dict(i => 0.0 for i in 1:n)
    return State{Int}(0, 0.0, 0.0, zero_counts, zero_rewards, copy(zero_rewards))
  end

  function initial_trace(instance::CombinatorialInstance{Int})
    return Trace{Int}(State{Int}[], Vector{Int}[], Vector{Float64}[], PolicyDetails[], Int[])
  end

  function initial_trace(instance::CombinatorialInstance{Tuple{Int, Int}})
    return Trace{Tuple{Int, Int}}(State{Tuple{Int, Int}}[], Vector{Tuple{Int, Int}}[], Vector{Float64}[], PolicyDetails[], Int[])
  end

  # Interface for policies.
  function choose_action(instance::CombinatorialInstance{T}, policy::Policy, state::State{T}) where T end

  # Update the state before the new round.
  function update!(state::State{T}, instance::CombinatorialInstance{T}, arms::Vector{T}, reward::Vector{Float64}, incurred_regret::Float64) where T
    state.round += 1

    # One reward per arm, i.e. semi-bandit feedback (not bandit feedback, where there would be only one reward for all arms).
    for i in 1:length(arms)
      state.arm_counts[arms[i]] += 1
      state.arm_reward[arms[i]] += reward[i]
      state.arm_average_reward[arms[i]] = state.arm_reward[arms[i]] / state.arm_counts[arms[i]]
      state.reward += reward[i]
    end

    state.regret += incurred_regret
  end

  # Use the bandit for the given number of steps.
  function simulate(instance::CombinatorialInstance{T}, policy::Policy, steps::Int; with_trace::Bool=false) where T
    state = initial_state(instance)
    if with_trace
      trace = initial_trace(instance)
    end

    for i in 1:steps
      t0 = time_ns()
      if with_trace
        arms, run_details = choose_action(instance, policy, state, with_trace=true)
      else
        arms = choose_action(instance, policy, state, with_trace=false)
      end
      t1 = time_ns()

      if length(arms) == 0
        error("No arms have been chosen at round $(i)!")
      end

      reward, incurred_regret = pull(instance, arms)
      update!(state, instance, arms, reward, incurred_regret)

      if with_trace
        push!(trace, state, arms, reward, run_details, round(Int, (t1 - t0) / 1_000_000_000))
      end

      if i % 100 == 0
        println(i)
      end
    end

    if ! with_trace
      return state
    else
      return state, trace
    end
  end

  ## Combinatorial algorithms. TODO: put this in another package.
  include("algos/msets.jl")
  include("algos/msets_budgeted.jl")
  include("algos/lp.jl")
  include("algos/lp_budgeted.jl")
  include("algos/st.jl")
  include("algos/st_budgeted.jl")

  ## Bandit policies.
  include("policies/thompson.jl")
  include("policies/llr.jl")
  include("policies/cucb.jl")
  include("policies/escb2.jl")
  include("policies/olsucb.jl")

  include("policies/escb2_exact.jl")
  include("policies/escb2_greedy.jl")
  include("policies/escb2_budgeted.jl")

  include("policies/olsucb_greedy.jl")

  ## Potential problems to solve.
  include("instances/perfectbipartitematching.jl")
  include("instances/perfectbipartitematching_lp.jl")
  include("instances/perfectbipartitematching_munkres.jl")
  include("instances/perfectbipartitematching_hungarian.jl")
  # Not using LightGraphsMatching, due to BlossomV dependency (hard to get to work…)

  include("instances/elementarypath.jl")
  include("instances/elementarypath_algos.jl")
  include("instances/elementarypath_lightgraphsdijkstra.jl")
  include("instances/elementarypath_lp.jl")

  include("instances/spanningtree.jl")
  include("instances/spanningtree_algos.jl")
  include("instances/spanningtree_lightgraphsprim.jl")
  include("instances/spanningtree_lp.jl")

  include("instances/mset.jl")
  include("instances/mset_algos.jl")
  include("instances/mset_lp.jl")
end


text:
module CreateTableSearches

import Migration: create_sequence, create_table, column_id, column,
                  column_id_sequence, add_index, drop_table, drop_sequence,
                  constraint, nextval

function up()
  create_sequence(:searches, :id)
  create_table(:searches) do
    [
      column_id(constraint = constraint(:searches, :id), nextval = nextval(:searches, :id))
      column(:search, :string, limit = 50)
    ]
  end
  column_id_sequence(:searches, :id)
  add_index(:searches, :search)
end

function down()
  drop_table(:searches)
end

end

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
#=
    Electrolyzer modelling
 =#

mutable struct Electrolyzer <: AbstractConverter
     # Paramètres
     α_p::Float64
     η_E_H2::Float64
     η_E_H::Float64
     lifetime::Int64
     nHoursMax::Float64
     bounds::NamedTuple{(:lb, :ub), Tuple{Float64, Float64}}
     # Initial conditions
     powerMax_ini::Float64
     soh_ini::Float64
     # Variables
     powerMax::AbstractArray{Float64,2}
     carrier::Vector{EnergyCarrier}
     soh::AbstractArray{Float64,3}
     # Eco
     cost::AbstractArray{Float64,2}
     # Inner constructor
     Electrolyzer(; α_p = 5/100,
                 η_E_H2 = 0.5,
                 η_E_H = 0.3,
                 lifetime = 15,
                 nHoursMax = 26000.,
                 bounds = (lb = 0., ub = 50.),
                 powerMax_ini = 0.,
                 soh_ini = 1.) =
                 new(α_p, η_E_H2, η_E_H, lifetime, nHoursMax, bounds, powerMax_ini, soh_ini)
end

### Preallocation
function preallocate!(elyz::Electrolyzer, nh::Int64, ny::Int64, ns::Int64)
     elyz.powerMax = convert(SharedArray,zeros(ny+1, ns)) ; elyz.powerMax[1,:] .= elyz.powerMax_ini
     elyz.carrier = [Electricity(), Heat(), Hydrogen()]
     elyz.carrier[1].power = convert(SharedArray,zeros(nh, ny, ns))
     elyz.carrier[2].power = convert(SharedArray,zeros(nh, ny, ns))
     elyz.carrier[3].power = convert(SharedArray,zeros(nh, ny, ns))
     elyz.soh = convert(SharedArray,zeros(nh+1, ny+1, ns)) ; elyz.soh[1,1,:] .= elyz.soh_ini
     elyz.cost = convert(SharedArray,zeros(ny, ns))
     return elyz
end

### Operation dynamic
function compute_operation_dynamics!(h::Int64, y::Int64, s::Int64, elyz::Electrolyzer, decision::Float64, Δh::Int64)
    elyz.soh[h+1,y,s], elyz.carrier[1].power[h,y,s], elyz.carrier[2].power[h,y,s], elyz.carrier[3].power[h,y,s] =
    compute_operation_dynamics(elyz, (powerMax = elyz.powerMax[y,s], soh = elyz.soh[h,y,s]), decision, Δh)
end

function compute_operation_dynamics(elyz::Electrolyzer, state::NamedTuple{(:powerMax, :soh), Tuple{Float64, Float64}}, decision::Float64, Δh::Int64)
    # Power constraint and correction
    elyz.α_p * state.powerMax >= decision && state.soh * elyz.nHoursMax / Δh > 1. ? power_E = max(decision, -state.powerMax) : power_E = 0.
    # Power conversion
    power_H2 = - power_E * elyz.η_E_H2
    power_H = - power_E * elyz.η_E_H
    # SoH computation
    soh_next = state.soh - (power_E > 0.) * Δh / elyz.nHoursMax
    return soh_next, power_E, power_H, power_H2
end
### Investment dynamic
function compute_investment_dynamics!(y::Int64, s::Int64, elyz::Electrolyzer, decision::Union{Float64, Int64})
    elyz.powerMax[y+1,s], elyz.soh[1,y+1,s] = compute_investment_dynamics(elyz, (powerMax = elyz.powerMax[y,s], soh = elyz.soh[end,y,s]), decision)
end

function compute_investment_dynamics(elyz::Electrolyzer, state::NamedTuple{(:powerMax, :soh), Tuple{Float64, Float64}}, decision::Union{Float64, Int64})
    if decision > 1e-2
        powerMax_next = decision
        soh_next = 1.
    else
        powerMax_next = state.powerMax
        soh_next = state.soh
    end
    return powerMax_next, soh_next
end


text:
export BrickletServoV2Status
struct BrickletServoV2Status
    enabled::Vector{Bool}
    current_position::Vector{Integer}
    current_velocity::Vector{Integer}
    current::Vector{Integer}
    input_voltage::Integer
end

export BrickletServoV2MotionConfiguration
struct BrickletServoV2MotionConfiguration
    velocity::Integer
    acceleration::Integer
    deceleration::Integer
end

export BrickletServoV2PulseWidth
struct BrickletServoV2PulseWidth
    min::Integer
    max::Integer
end

export BrickletServoV2Degree
struct BrickletServoV2Degree
    min::Integer
    max::Integer
end

export BrickletServoV2SPITFPErrorCount
struct BrickletServoV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletServoV2Identity
struct BrickletServoV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletServoV2
"""
Drives up to 10 RC Servos
"""
mutable struct BrickletServoV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletServoV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_servo_v2")
        deviceInternal = package.BrickletServoV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_status
"""
    $(SIGNATURES)

Returns the status information of the Servo Bricklet 2.0.

The status includes

* for each channel if it is enabled or disabled,
* for each channel the current position,
* for each channel the current velocity,
* for each channel the current usage and
* the input voltage.

Please note that the position and the velocity is a snapshot of the
current position and velocity of the servo in motion.
"""
function get_status(device::BrickletServoV2)
    return device.deviceInternal.get_status()
end

export set_enable
"""
    $(SIGNATURES)

Enables a servo channel (0 to 9). If a servo is enabled, the configured position,
velocity, acceleration, etc. are applied immediately.
"""
function set_enable(device::BrickletServoV2, servo_channel, enable)
    device.deviceInternal.set_enable(servo_channel, enable)
end

export get_enabled
"""
    $(SIGNATURES)

Returns *true* if the specified servo channel is enabled, *false* otherwise.
"""
function get_enabled(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_enabled(servo_channel)
end

export set_position
"""
    $(SIGNATURES)

Sets the position in °/100 for the specified servo channel.

The default range of the position is -9000 to 9000, but it can be specified
according to your servo with :func:`Set Degree`.

If you want to control a linear servo or RC brushless motor controller or
similar with the Servo Brick, you can also define lengths or speeds with
:func:`Set Degree`.
"""
function set_position(device::BrickletServoV2, servo_channel, position)
    device.deviceInternal.set_position(servo_channel, position)
end

export get_position
"""
    $(SIGNATURES)

Returns the position of the specified servo channel as set by :func:`Set Position`.
"""
function get_position(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_position(servo_channel)
end

export get_current_position
"""
    $(SIGNATURES)

Returns the *current* position of the specified servo channel. This may not be the
value of :func:`Set Position` if the servo is currently approaching a
position goal.
"""
function get_current_position(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_current_position(servo_channel)
end

export get_current_velocity
"""
    $(SIGNATURES)

Returns the *current* velocity of the specified servo channel. This may not be the
velocity specified by :func:`Set Motion Configuration`. if the servo is
currently approaching a velocity goal.
"""
function get_current_velocity(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_current_velocity(servo_channel)
end

export set_motion_configuration
"""
    $(SIGNATURES)

Sets the maximum velocity of the specified servo channel in °/100s as well as
the acceleration and deceleration in °/100s²

With a velocity of 0 °/100s the position will be set immediately (no velocity).

With an acc-/deceleration of 0 °/100s² the velocity will be set immediately
(no acc-/deceleration).
"""
function set_motion_configuration(device::BrickletServoV2, servo_channel, velocity, acceleration, deceleration)
    device.deviceInternal.set_motion_configuration(servo_channel, velocity, acceleration, deceleration)
end

export get_motion_configuration
"""
    $(SIGNATURES)

Returns the motion configuration as set by :func:`Set Motion Configuration`.
"""
function get_motion_configuration(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_motion_configuration(servo_channel)
end

export set_pulse_width
"""
    $(SIGNATURES)

Sets the minimum and maximum pulse width of the specified servo channel in µs.

Usually, servos are controlled with a
`PWM <https://en.wikipedia.org/wiki/Pulse-width_modulation>`__, whereby the
length of the pulse controls the position of the servo. Every servo has
different minimum and maximum pulse widths, these can be specified with
this function.

If you have a datasheet for your servo that specifies the minimum and
maximum pulse width, you should set the values accordingly. If your servo
comes without any datasheet you have to find the values via trial and error.

Both values have a range from 1 to 65535 (unsigned 16-bit integer). The
minimum must be smaller than the maximum.

The default values are 1000µs (1ms) and 2000µs (2ms) for minimum and
maximum pulse width.
"""
function set_pulse_width(device::BrickletServoV2, servo_channel, min, max)
    device.deviceInternal.set_pulse_width(servo_channel, min, max)
end

export get_pulse_width
"""
    $(SIGNATURES)

Returns the minimum and maximum pulse width for the specified servo channel as set by
:func:`Set Pulse Width`.
"""
function get_pulse_width(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_pulse_width(servo_channel)
end

export set_degree
"""
    $(SIGNATURES)

Sets the minimum and maximum degree for the specified servo channel (by default
given as °/100).

This only specifies the abstract values between which the minimum and maximum
pulse width is scaled. For example: If you specify a pulse width of 1000µs
to 2000µs and a degree range of -90° to 90°, a call of :func:`Set Position`
with 0 will result in a pulse width of 1500µs
(-90° = 1000µs, 90° = 2000µs, etc.).

Possible usage:

* The datasheet of your servo specifies a range of 200° with the middle position
  at 110°. In this case you can set the minimum to -9000 and the maximum to 11000.
* You measure a range of 220° on your servo and you don't have or need a middle
  position. In this case you can set the minimum to 0 and the maximum to 22000.
* You have a linear servo with a drive length of 20cm, In this case you could
  set the minimum to 0 and the maximum to 20000. Now you can set the Position
  with :func:`Set Position` with a resolution of cm/100. Also the velocity will
  have a resolution of cm/100s and the acceleration will have a resolution of
  cm/100s².
* You don't care about units and just want the highest possible resolution. In
  this case you should set the minimum to -32767 and the maximum to 32767.
* You have a brushless motor with a maximum speed of 10000 rpm and want to
  control it with a RC brushless motor controller. In this case you can set the
  minimum to 0 and the maximum to 10000. :func:`Set Position` now controls the rpm.

Both values have a possible range from -32767 to 32767
(signed 16-bit integer). The minimum must be smaller than the maximum.

The default values are -9000 and 9000 for the minimum and maximum degree.
"""
function set_degree(device::BrickletServoV2, servo_channel, min, max)
    device.deviceInternal.set_degree(servo_channel, min, max)
end

export get_degree
"""
    $(SIGNATURES)

Returns the minimum and maximum degree for the specified servo channel as set by
:func:`Set Degree`.
"""
function get_degree(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_degree(servo_channel)
end

export set_period
"""
    $(SIGNATURES)

Sets the period of the specified servo channel in µs.

Usually, servos are controlled with a
`PWM <https://en.wikipedia.org/wiki/Pulse-width_modulation>`__. Different
servos expect PWMs with different periods. Most servos run well with a
period of about 20ms.

If your servo comes with a datasheet that specifies a period, you should
set it accordingly. If you don't have a datasheet and you have no idea
what the correct period is, the default value (19.5ms) will most likely
work fine.

The minimum possible period is 1µs and the maximum is 1000000µs.

The default value is 19.5ms (19500µs).
"""
function set_period(device::BrickletServoV2, servo_channel, period)
    device.deviceInternal.set_period(servo_channel, period)
end

export get_period
"""
    $(SIGNATURES)

Returns the period for the specified servo channel as set by :func:`Set Period`.
"""
function get_period(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_period(servo_channel)
end

export get_servo_current
"""
    $(SIGNATURES)

Returns the current consumption of the specified servo channel in mA.
"""
function get_servo_current(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_servo_current(servo_channel)
end

export set_servo_current_configuration
"""
    $(SIGNATURES)

Sets the averaging duration of the current measurement for the specified servo channel in ms.
"""
function set_servo_current_configuration(device::BrickletServoV2, servo_channel, averaging_duration)
    device.deviceInternal.set_servo_current_configuration(servo_channel, averaging_duration)
end

export get_servo_current_configuration
"""
    $(SIGNATURES)

Returns the servo current configuration for the specified servo channel as set
by :func:`Set Servo Current Configuration`.
"""
function get_servo_current_configuration(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_servo_current_configuration(servo_channel)
end

export set_input_voltage_configuration
"""
    $(SIGNATURES)

Sets the averaging duration of the input voltage measurement for the specified servo channel in ms.
"""
function set_input_voltage_configuration(device::BrickletServoV2, averaging_duration)
    device.deviceInternal.set_input_voltage_configuration(averaging_duration)
end

export get_input_voltage_configuration
"""
    $(SIGNATURES)

Returns the input voltage configuration as set by :func:`Set Input Voltage Configuration`.
"""
function get_input_voltage_configuration(device::BrickletServoV2)
    return device.deviceInternal.get_input_voltage_configuration()
end

export get_overall_current
"""
    $(SIGNATURES)

Returns the current consumption of all servos together in mA.
"""
function get_overall_current(device::BrickletServoV2)
    return device.deviceInternal.get_overall_current()
end

export get_input_voltage
"""
    $(SIGNATURES)

Returns the input voltage in mV. The input voltage is
given via the black power input connector on the Servo Brick.
"""
function get_input_voltage(device::BrickletServoV2)
    return device.deviceInternal.get_input_voltage()
end

export set_current_calibration
"""
    $(SIGNATURES)

Sets an offset value (in mA) for each channel.

Note: On delivery the Servo Bricklet 2.0 is already calibrated.
"""
function set_current_calibration(device::BrickletServoV2, offset)
    device.deviceInternal.set_current_calibration(offset)
end

export get_current_calibration
"""
    $(SIGNATURES)

Returns the current calibration as set by :func:`Set Current Calibration`.
"""
function get_current_calibration(device::BrickletServoV2)
    return device.deviceInternal.get_current_calibration()
end

export set_position_reached_callback_configuration
"""
    $(SIGNATURES)

Enable/Disable :cb:`Position Reached` callback.
"""
function set_position_reached_callback_configuration(device::BrickletServoV2, servo_channel, enabled)
    device.deviceInternal.set_position_reached_callback_configuration(servo_channel, enabled)
end

export get_position_reached_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set Position Reached Callback Configuration`.
"""
function get_position_reached_callback_configuration(device::BrickletServoV2, servo_channel)
    return device.deviceInternal.get_position_reached_callback_configuration(servo_channel)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletServoV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletServoV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletServoV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletServoV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletServoV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletServoV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletServoV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletServoV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletServoV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletServoV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletServoV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletServoV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletServoV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 5abfca08-6935-11eb-37b5-374ea3b4ddf0
begin
	using Revise
	using BSplineKit
	using SpecialFunctions: airyai
	using LinearAlgebra
	using BandedMatrices
	using QuadGK
	using Plots
end

# ╔═╡ 48367b6e-6955-11eb-3ea4-5b1f8099f8ab
md"""
# Airy equation

Example taken from Olver & Townsend, SIAM Rev. 2013.
"""

# ╔═╡ ab3e2c3c-6952-11eb-2c48-d5dea019ea78
N = 500

# ╔═╡ d1540532-6951-11eb-2df0-858510e4d960
breaks =
	# vcat(
	# 	range(-1, 0; length = N + 1),
	# 	range(0, 1; length = (N ÷ 8) + 1)[2:end],
	# )
	collect(range(-1, 1; length = N + 1))
	# [-cos(π * i / N) for i = 0:N]

# ╔═╡ ce14d81c-6936-11eb-2876-fd539b26f77c
B = BSplineBasis(BSplineOrder(6), copy(breaks))

# ╔═╡ e4b55060-6936-11eb-336d-7f6d75edabff
R = RecombinedBSplineBasis(Derivative(0), B)

# ╔═╡ 06d6d560-6937-11eb-03e3-5b2be0c46741
xs = collocation_points(B);

# ╔═╡ 947045b2-6937-11eb-0501-139f38e7783b
Cfact = lu!(collocation_matrix(B, xs));

# ╔═╡ 763b4c5a-6938-11eb-2969-958078a01dea
md"## Variable coefficient"

# ╔═╡ 86f7485a-6938-11eb-2775-3997a69d16dc
cs = Cfact \ xs

# ╔═╡ b9eb0044-6938-11eb-0ece-fb3c8466596d
Tc = galerkin_tensor((R, R, B), Derivative.((0, 0, 0)));

# ╔═╡ 2f232272-6945-11eb-3c84-4bf241954db3
Mc = Hermitian(Tc * cs);

# ╔═╡ 390c6438-694a-11eb-05ee-9da9e9f26d18
issymmetric(parent(Mc))

# ╔═╡ ca4a99fa-6946-11eb-3931-a7cb96c58684
summary(Mc)

# ╔═╡ ea8e4e44-6939-11eb-0a3c-4b669c5645c8
md"## LHS"

# ╔═╡ 69982892-6938-11eb-2972-333ba2799cf4
md"## RHS"

# ╔═╡ 32ba8d7e-6938-11eb-3056-b543422fa4f5
# RHS operator
Mf = galerkin_matrix((R, B));

# ╔═╡ b343c62e-694a-11eb-37bd-d14a3d58f8d8
md"## Solution"

# ╔═╡ 70f79dde-6938-11eb-13b3-f9775864234c
md"## Plots"

# ╔═╡ b71a4ed6-6935-11eb-2db1-c9bbdb97b9cc
u_exact(x; ε) = airyai(x / cbrt(ε))

# ╔═╡ 19b426e4-6937-11eb-2d24-436700fa442e
# Affine part of the solution
u_0(x; ε) = ((1 - x) * u_exact(-1; ε) + (1 + x) * u_exact(1; ε)) / 2

# ╔═╡ a9e506de-6935-11eb-39f0-fb44ea95ffd3
ε = 10^(-5)

# ╔═╡ f05d7e94-6939-11eb-3c0e-af2d31d8eee0
L = let
	L = galerkin_matrix(R, Derivative.((1, 1))) :: Hermitian
	A = parent(L)
	A .*= ε
	@assert L.uplo == 'U'
	A .+= UpperTriangular(Mc)
	BandedMatrix(L)  # convert to "full" banded matrix, for LU
end;

# ╔═╡ e345757a-694a-11eb-0c84-71197d96685d
Lfact = lu(L);

# ╔═╡ 8a214fa4-6937-11eb-382e-bf64ff050dbb
# RHS coefficients
fs = ldiv!(Cfact, xs .* u_0.(xs; ε))

# ╔═╡ 59760524-6938-11eb-0766-a193e2632d53
rhs = -Mf * fs

# ╔═╡ b9344f10-694a-11eb-21c5-c1bf6c0733e1
vR = Lfact \ rhs

# ╔═╡ ff706bd6-694a-11eb-10a5-01a9d40686f3
vB = recombination_matrix(R) * vR

# ╔═╡ 1b6b3ba6-694b-11eb-0f4e-a18cc5b6e991
Sv = Spline(B, vB);

# ╔═╡ 761215d4-694b-11eb-1089-a7ccde18df5b
usol(x; ε) = u_0(x; ε) + Sv(x)

# ╔═╡ 4a7dcc82-6951-11eb-1c74-ddfda1515953
uerr(x; ε) = (usol(x; ε) - u_exact(x; ε))^2

# ╔═╡ 079e6a38-6952-11eb-2d5f-a9d39f6ca4e0
quadgk(x -> uerr(x; ε), -1, 1; rtol = 1e-4)

# ╔═╡ 4b056762-6935-11eb-116e-47a5a46dbfff
let
	fig = plot(x -> u_exact(x; ε), -1, 1)
	# plot!(ax, -1..1, x -> u_0(x; ε); color = :red)
	# plot!(ax, -1..1, x -> Spline(B, fs)(x); color = :blue)
	plot!(fig, x -> usol(x; ε), -1, 1; color = :orange)
	fig
end

# ╔═╡ ff09558c-694b-11eb-379e-0f3b985db80a
let
	fig = plot(x -> log10.(uerr.(x; ε)), -1, 1; ylims = (-20, 1))
end

# ╔═╡ Cell order:
# ╟─48367b6e-6955-11eb-3ea4-5b1f8099f8ab
# ╠═ab3e2c3c-6952-11eb-2c48-d5dea019ea78
# ╠═d1540532-6951-11eb-2df0-858510e4d960
# ╠═ce14d81c-6936-11eb-2876-fd539b26f77c
# ╠═e4b55060-6936-11eb-336d-7f6d75edabff
# ╠═06d6d560-6937-11eb-03e3-5b2be0c46741
# ╠═947045b2-6937-11eb-0501-139f38e7783b
# ╟─763b4c5a-6938-11eb-2969-958078a01dea
# ╠═86f7485a-6938-11eb-2775-3997a69d16dc
# ╠═b9eb0044-6938-11eb-0ece-fb3c8466596d
# ╠═2f232272-6945-11eb-3c84-4bf241954db3
# ╠═390c6438-694a-11eb-05ee-9da9e9f26d18
# ╠═ca4a99fa-6946-11eb-3931-a7cb96c58684
# ╟─ea8e4e44-6939-11eb-0a3c-4b669c5645c8
# ╠═f05d7e94-6939-11eb-3c0e-af2d31d8eee0
# ╠═e345757a-694a-11eb-0c84-71197d96685d
# ╟─69982892-6938-11eb-2972-333ba2799cf4
# ╠═8a214fa4-6937-11eb-382e-bf64ff050dbb
# ╠═32ba8d7e-6938-11eb-3056-b543422fa4f5
# ╠═59760524-6938-11eb-0766-a193e2632d53
# ╟─b343c62e-694a-11eb-37bd-d14a3d58f8d8
# ╠═b9344f10-694a-11eb-21c5-c1bf6c0733e1
# ╠═ff706bd6-694a-11eb-10a5-01a9d40686f3
# ╠═1b6b3ba6-694b-11eb-0f4e-a18cc5b6e991
# ╠═761215d4-694b-11eb-1089-a7ccde18df5b
# ╠═4a7dcc82-6951-11eb-1c74-ddfda1515953
# ╠═079e6a38-6952-11eb-2d5f-a9d39f6ca4e0
# ╟─70f79dde-6938-11eb-13b3-f9775864234c
# ╠═4b056762-6935-11eb-116e-47a5a46dbfff
# ╠═ff09558c-694b-11eb-379e-0f3b985db80a
# ╠═b71a4ed6-6935-11eb-2db1-c9bbdb97b9cc
# ╠═19b426e4-6937-11eb-2d24-436700fa442e
# ╠═a9e506de-6935-11eb-39f0-fb44ea95ffd3
# ╠═5abfca08-6935-11eb-37b5-374ea3b4ddf0

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
function generate_random_specification(LSM_ratio, porosity)
  v = Array{Int16}(undef, 100)

  hole_length = Int32(round(porosity*100))
  LSM_length = Int32(round((1 - porosity)*LSM_ratio*100))

  v .= i_YSZ
  v[1 : hole_length] .= i_hole
  v[hole_length + 1 : hole_length + LSM_length] .= i_LSM
  return v
end

# generate random matrix of dimensions
function generate_matrix(dimensions::Tuple{T, T} where T <: Integer, porosity::Float64, LSM_ratio::Float64)
  return rand(
                          generate_random_specification(LSM_ratio, porosity),
                          dimensions...
         )
end

function generate_matrix(dimensions::Tuple{T, T, T} where T <: Integer, porosity::Float64, LSM_ratio::Float64)
  res = Array{Int64}(undef, dimensions...)
  for layer in 1:dimensions[3]
    res[:, :, layer] = rand(
                          generate_random_specification(LSM_ratio, porosity),
                          dimensions[1], dimensions[2]
                      )
  end
  return res
end

function generate_submatrix_to_matrix(matrix, left_upper::Union{Tuple, Array}, right_lower::Union{Tuple, Array}, porosity::Float64, LSM_ratio::Float64)
  submatrix = generate_matrix(right_lower .- left_upper .+ (1,1), porosity, LSM_ratio)
  matrix[left_upper[1] : right_lower[1], left_upper[2] : right_lower[2]] = deepcopy(submatrix)
  return
end

function generate_submatrix_to_matrix(matrix, right_lower::Union{Tuple, Array}, porosity::Float64, LSM_ratio::Float64)
  generate_submatrix_to_matrix(matrix, (1,1), right_lower::Union{Tuple, Array}, porosity::Float64, LSM_ratio::Float64)
end

function get_default_submatrix_list()
  # 1 item in the resulting list is
  #
  # left upper, right lower corner, porosity, LSM_ratio

  first_block_column = 1
  second_block_column = 21
  third_block_column = 41
  width = 42
  height = 80
  return [
    [(1,1), (height, first_block_column), 1.0, 0.5],
    [(1,first_block_column + 1), (height, second_block_column), 0.2, 0.5],
    [(1,second_block_column + 1), (height, third_block_column), 0.2, 1.0],
    [(1,third_block_column + 1), (height, width), 0.2, 1.0],
    #
    [(10, 1), (15, first_block_column), 0.0, 1.0],
    [(40, 1), (45, first_block_column), 0.0, 1.0],
    [(60, 1), (65, first_block_column), 0.0, 1.0]
  ]
end

function generate_matrix(submatrix_list::Array=get_default_submatrix_list())
  max_height, max_width = -1, -1
  for submatrix in submatrix_list
    s = submatrix[2]
    if s[1] > max_height
      max_height = s[1]
    end
    if s[2] > max_width
      max_width = s[2]
    end
  end

  output_matrix = Matrix{Int16}(undef, max_height, max_width)
  output_matrix .= -1

  for submatrix in submatrix_list
    generate_submatrix_to_matrix(output_matrix, submatrix...)
  end
  invalid_submatrix_list = false
  foreach(x -> x==-1 ?  invalid_submatrix_list = true : false, output_matrix)
  if invalid_submatrix_list
    println("ERROR: invalid_submatrix_list input")
    return throw(Exception)
  else
    return output_matrix
  end
end

function three_column_domain_template(LSM_ratio1, LSM_ratio2, LSM_ratio3;
                              #
                              porosity1=0.5, porosity2=0.5, porosity3=0.5,
                              #
                              positions_of_contacts=[15, 50], height_of_contacts=5,
                              #
                              column_width = 5,
                              #
                              height = 70
                              )
  # 1 item in the resulting list is
  #
  # left upper, right lower corner, porosity, LSM_ratio
  first_block_column = 1
  second_block_column = column_width + 1
  third_block_column = 2*column_width + 1
  fourth_block_column = 3*column_width + 1
  width = fourth_block_column + 1

  output = [
    [(1,1), (height, first_block_column), 1.0, 0.5],
    [(1,first_block_column + 1), (height, second_block_column), porosity1, LSM_ratio1],
    [(1,second_block_column + 1), (height, third_block_column), porosity2, LSM_ratio2],
    [(1,third_block_column + 1), (height, fourth_block_column), porosity3, LSM_ratio3],
    [(1, fourth_block_column + 1), (height, width), 0.0, 1.0]
  ]

  # contacts
  for p_of_contact in positions_of_contacts
    push!(output,
        [(p_of_contact, 1), (p_of_contact + height_of_contacts - 1, first_block_column), 0.0, 1.0]
    )
  end

  return output
end

function chess_matrix(m, n)
  A = Matrix{Int64}(undef, m,n)
  A .= 0
  for i in 1:m, j in 1:n
      if mod(i+j, 2) == 0
        A[i,j] = 1
      end
  end
  return A
end

max_stars_repo_name: Masicko/ImageToEIS

max_issues_repo_name: Masicko/ImageToEIS

max_forks_repo_name: Masicko/ImageToEIS


text:
using StaticArrays
export IterControl, itercontrol, controldo, group_shift!, lmove

# NOTE: use SortedVector in Blocks would help benchmarks

"""
    IterControl{S}
    IterControl(n::Int, base::Int, masks, ks) -> IterControl

Iterator to iterate through controlled subspace. See also [`itercontrol`](@ref).
 `S` is the number of shifts,
 `n` is the size of Hilbert space,
 `base` is the base of counter,
 `masks` and `ks` are helpers for enumerating over the target Hilbert Space.
"""
struct IterControl{S}
    n::Int
    base::Int
    masks::NTuple{S, Int}
    ks::NTuple{S, Int}

    function IterControl(n::Int, base::Int, masks::NTuple{S, Int}, ks::NTuple{S, Int}) where {S}
        new{S}(n, base, masks, ks)
    end
end

"""
    itercontrol([T=Int], nbits, positions, bit_configs)

Returns an iterator which iterate through controlled subspace of bits.

# Example

To iterate through all the bits satisfy `0xx10x1` where `x` means an arbitrary bit.

```jldoctest
julia> for each in itercontrol(7, [1, 3, 4, 7], (1, 0, 1, 0))
           println(string(each, base=2, pad=7))
       end
0001001
0001011
0011001
0011011
0101001
0101011
0111001
0111011

```
"""
# NOTE: positions should be vector (MVector is the best), since it need to be sorted
#       do not use Tuple, or other immutables, it increases the sorting time.
function itercontrol(nbits::Int, positions::AbstractVector, bit_configs) where T
    base = bmask(Int, positions[i] for (i, u) in enumerate(bit_configs) if u != 0)
    masks, ks = group_shift!(nbits, positions)
    S = length(masks)
    return IterControl(1<<(nbits - length(positions)), base, Tuple(masks), Tuple(ks))
end

"""
    controldo(f, itr::IterControl)

Execute `f` while iterating `itr`.

!!! note

    this is faster but equivalent than using `itr` as an iterator.
    See also [`itercontrol`](@ref).
"""
function controldo(f::Base.Callable, ic::IterControl{S}) where {S}
    for i in 0:ic.n-1
        @simd for s in 1:S
            @inbounds i = lmove(i, ic.masks[s], ic.ks[s])
        end
        f(i+ic.base)
    end
    return nothing
end

Base.length(it::IterControl) = it.n
Base.eltype(it::IterControl) = Int

function Base.getindex(it::IterControl{S}, k::Int) where {S}
    out = k - 1
    @simd for s in 1:S
        @inbounds out = lmove(out, it.masks[s], it.ks[s])
    end
    return out + it.base
end

function Base.iterate(it::IterControl{S}, state = 1) where {S}
    if state > length(it)
        return nothing
    else
        return it[state], state + 1
    end
end

lmove(b::Int, mask::Int, k::Int)::Int = (b&~mask)<<k + (b&mask)

"""
    group_shift!(nbits, positions)

Shift bits on `positions` together.
"""
function group_shift!(nbits::Int, positions::AbstractVector{Int})
    sort!(positions)
    masks = Int[]; ns = Int[]
    k_prv = -1
    for k in positions
        if k == k_prv+1
            ns[end] += 1
        else
            push!(masks, bmask(0:k-1))
            push!(ns, 1)
        end
        k_prv = k
    end
    return masks, ns
end


text:
"""
    rand_init(x::Array; range = (0, 1))

Generate a random initialization of the same size and type as `x`.
If `x` is an integer array, then `range` specifies a set of integers.
Otherwise, `range` specifies the set of real numbers.
"""
rand_init(x::AbstractArray{<:Integer}; range = (0, 255)) = rand(range[1]:range[2], size(x)...)
rand_init(x::AbstractArray; range = (0, 1)) = (range[2] - range[1]) .* rand(size(x)...) .+ range[1]

"""
    proj_lball!(x, δ; ϵ, ϵnorm)

Project the perturbed sample (`x + δ`) onto an `ϵnorm`-ball of size `ϵ` around `x`.

# Arguments
- `x`: the adversarial (perturbed) sample/batch
- `δ`: the pertubations
- `ϵ`: the size of the L-norm ball
- `ϵnorm`: the type of L-norm used
"""
function proj_lball!(x, δ; ϵ, ϵnorm)
    if isinf(ϵnorm)
        x .= x .+ max.(min.(δ, ϵ), -ϵ)
    else
        map!((xi, δi) -> xi .+ δi * (ϵ / max.(norm(reshape(δi, :), ϵnorm), ϵ)),
             x, eachslice(x; dims = ndims(x)), eachslice(δ; dims = ndims(δ)))
    end

    return x
end


text:
using Test

@testset "Images" begin
include("arrays.jl")
include("algorithms.jl")
include("exposure.jl")
include("edge.jl")
include("corner.jl")
include("bwdist.jl")
include("writemime.jl")
# include("deprecated.jl")

end