text:
module TemporalGPs

    using AbstractGPs
    using BlockDiagonals
    using ChainRulesCore
    using FillArrays
    using LinearAlgebra
    using KernelFunctions
    using Random
    using StaticArrays
    using StructArrays
    using Zygote
    using ZygoteRules

    using FillArrays: AbstractFill
    using Zygote: _pullback, AContext

    import AbstractGPs: mean, cov, logpdf, FiniteGP, AbstractGP, posterior, dtc, elbo

    using KernelFunctions:
        SimpleKernel,
        KernelSum,
        ScaleTransform,
        ScaledKernel,
        TransformedKernel

    export
        to_sde,
        SArrayStorage,
        ArrayStorage,
        RegularSpacing,
        checkpointed,
        posterior,
        logpdf_and_rand,
        Separable

    # Various bits-and-bobs. Often commiting some type piracy.
    include(joinpath("util", "harmonise.jl"))
    include(joinpath("util", "linear_algebra.jl"))
    include(joinpath("util", "scan.jl"))
    include(joinpath("util", "zygote_friendly_map.jl"))
    include(joinpath("util", "zygote_rules.jl"))
    include(joinpath("util", "gaussian.jl"))
    include(joinpath("util", "mul.jl"))
    include(joinpath("util", "storage_types.jl"))
    include(joinpath("util", "regular_data.jl"))

    # Linear-Gaussian State Space Models.
    include(joinpath("models", "linear_gaussian_conditionals.jl"))
    include(joinpath("models", "gauss_markov_model.jl"))
    include(joinpath("models", "lgssm.jl"))
    include(joinpath("models", "missings.jl"))

    # Converting GPs to Linear-Gaussian SSMs.
    include(joinpath("gp", "data_representations.jl"))
    include(joinpath("gp", "lti_sde.jl"))
    include(joinpath("gp", "posterior_lti_sde.jl"))

    # Converting space-time GPs to Linear-Gaussian SSMs.
    include(joinpath("space_time", "rectilinear_grid.jl"))
    include(joinpath("space_time", "regular_in_time.jl"))
    include(joinpath("space_time", "separable_kernel.jl"))
    include(joinpath("space_time", "to_gauss_markov.jl"))
    include(joinpath("space_time", "pseudo_point.jl"))
end # module


text:
using JuMP, Gurobi, Distributions, StatsFuns, StatsBase

include("resource_nonlinear_20_30.jl")
include("resource_initialization_template.jl")

const gurobi_env = Gurobi.Env()

const numVariables = numResources

# generate samples of the random variables
function generateRandomSamples(numSamples::Int64)

	lambda = rand(MvNormal(lambda_mean[1,:],lambda_covariance),numSamples)
	rho = Array{Float64}(numResources,numSamples)
	for i = 1:numResources
		rho[i,:] = rand(Normal(rho_mean[i],rho_stdev),numSamples)
		for samp = 1:numSamples
			rho[i,samp] = min(rho[i,samp],rho_max)
		end
	end

	return lambda, rho
end

srand(1234)

# generate samples to estimate actual risk level of solution
const numAnalyticalSamples = 10000
const lambda_analytical, rho_analytical = generateRandomSamples(numAnalyticalSamples)

srand()

# compute conservative bound on the probability that
# constraints aren't satisfied
function computeRiskLevel(x::Array{Float64},reliabilityLevel::Float64=1E-06)

	const numSamples = size(lambda_analytical,2)
	simpleEst::Bool = true

	numViolated::Int64 = checkScenarioConstraints(x,lambda_analytical,rho_analytical)

	riskLevel::Float64 = Inf

	if(simpleEst)
		riskLevel = numViolated*1.0/numSamples
		return riskLevel
	else
		if(numViolated == 0)
			riskLevel = -log(reliabilityLevel)/numSamples
			return riskLevel
		end

		 # termination criterion for bisection
		gamma_tolerance::Float64 = 0.1/numSamples

		# use tail bounds on binomial distribution to estimate lower and upper bounds
		gamma_low::Float64 = (numViolated - sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_low < 0.0)
			gamma_low = 0.0
		end
		gamma_up::Float64 = (numViolated + sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_up > 1.0)
			gamma_up = 1.0
		end

		gamma::Float64 = 0.0
		func_value::Float64 = 0.0

		# use bisection to solve the nonlinear equation
		while (gamma_up - gamma_low > gamma_tolerance)

			gamma = (gamma_low + gamma_up)/2.0
			func_value = -reliabilityLevel
			for i = 1:numViolated
				tmp_func = i*log(gamma) + (numSamples-i)*log(1-gamma) + lfact(numSamples) - lfact(numSamples-i) - lfact(i)
				func_value += exp(tmp_func)
				if(func_value > 0)
					break
				end
			end
			if(func_value > 0)
				gamma_low = (gamma_low + gamma_up)/2.0
			else
				gamma_up = (gamma_low + gamma_up)/2.0
			end

		end

		riskLevel = gamma_up
		return riskLevel
	end

end

# determines if all constraints are satisfied for a given single scenario
# returns number of violated scenarios
function solveBasicModel(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64})

	const numSamples = size(lambda,2)
	obj = zeros(Float64,numSamples)
	dualmult = zeros(Float64,numResources,numSamples)

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, y[1:numResources,1:numCustomers] >= 0)
	@variable(mod, slack)

	@objective(mod, Min, slack)
	@constraint(mod, con1[i=1:numResources], sum(y[i,j] for j = 1:numCustomers) <= 0)
	@constraint(mod, con2[j=1:numCustomers], sum(nu[i,j]*y[i,j] for i = 1:numResources) + slack >= 0)

	for iter = 1:numSamples
		for i = 1:numResources
			JuMP.setRHS(con1[i],rho[i,iter]*(x[i])^2)
		end

		for j = 1:numCustomers
			JuMP.setRHS(con2[j],lambda[j,iter])
		end

		status = solve(mod)

		obj[iter] = getvalue(slack)
		dualmult[:,iter] = getdual(con1)
	end

	mod = 0

	return obj, dualmult
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables, and constraint scalings
function evaluateConstraint(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64},scalings::Float64)

	con, alpha_opt = solveBasicModel(x,lambda,rho)

	con /= scalings

	return con
end

# evaluate the constraint values for given decision vector x and
# realization of the random variables
function evaluateConstraint(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64})

	return evaluateConstraint(x,lambda,rho,1.0)
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64})

	return evaluateConstraintGradients(x,lambda,rho,1.0)
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64},scalings::Float64)

	const numSamples = size(lambda,2)
	grad_con = zeros(Float64,numResources,numSamples)

	con, alpha_opt = solveBasicModel(x,lambda,rho)

	for iter = 1:numSamples
		grad_con[:,iter] = 2*alpha_opt[:,iter].*rho[:,iter].*x/scalings
	end

	return grad_con
end

# evaluate constraint gradients
function evaluateConstraintAndGradient(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64},scalings::Float64)

	const numSamples = size(lambda,2)
	con = zeros(Float64,numSamples)
	grad_con = zeros(Float64,numResources,numSamples)

	con, alpha_opt = solveBasicModel(x,lambda,rho)
	con /= scalings

	for iter = 1:numSamples
		grad_con[:,iter] = 2*alpha_opt[:,iter].*rho[:,iter].*x/scalings
	end

	return con, grad_con
end

# determines if all constraints are satisfied for a given single scenario
# returns number of violated scenarios
function checkScenarioConstraints(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64})

	const numSamples = size(lambda,2)
	numViolated::Int64 = 0

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, y[1:numResources,1:numCustomers] >= 0)

	@objective(mod, Min, 0)
	@constraint(mod, con1[i=1:numResources], sum(y[i,j] for j = 1:numCustomers) <= 0)
	@constraint(mod, con2[j=1:numCustomers], sum(nu[i,j]*y[i,j] for i = 1:numResources) >= 0)

	for iter = 1:numSamples
		for i = 1:numResources
			JuMP.setRHS(con1[i],rho[i,iter]*(x[i])^2)
		end

		for j = 1:numCustomers
			JuMP.setRHS(con2[j],lambda[j,iter])
		end

		status = solve(mod,suppress_warnings=true)

		if(status != Symbol("Optimal"))
			numViolated += 1
		end
	end

	mod = 0

	return numViolated
end

# project the current iterate onto the feasible set
function projectOntoFeasibleSet(x_curr::Array{Float64},objectiveBound::Float64,proj_params::Array{Float64})

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, y[1:numResources] >= 0)

	@objective(mod, Min, sum((y[i]-x_curr[i])^2 for i = 1:numResources))
	@constraint(mod, sum(costVector[i]*y[i] for i=1:numResources) <= objectiveBound)

	solve(mod)
	x_proj = getvalue(y)

	return x_proj, proj_params
end

# get the stochastic gradient of the approximation
function getStochasticGradient(x::Array{Float64},numSamples::Int64,mu::Float64,tau::Float64,scalings::Float64)

	stochgrad = zeros(numResources)

	lambda, rho = generateRandomSamples(numSamples)

	con, grad_con = evaluateConstraintAndGradient(x,lambda,rho,scalings)

	baseFactor1 = tau./((sqrt(mu)*exp.(0.5*tau*con) + sqrt(1.0/mu)*exp.(-0.5*tau*con)).^2)

	stochgrad = (grad_con*baseFactor1)/numSamples

	return stochgrad

end

# get the stochastic gradient of the approximation for fixed random sample
function getConstraintStochasticGradient(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64},mu::Float64,tau::Float64,scalings::Float64)

	numSamples::Int64 = size(lambda,2)

	stochgrad = zeros(Float64,numResources)

	grad_con = evaluateConstraintGradients(x,lambda,rho,scalings)

	stochgrad = sum(grad_con,2)/numSamples

	return stochgrad

end

# estimate constraint scalings
function estimateScalings(x::Array{Float64},numSamples::Int64)

	scaling_tol::Float64 = 1E-06

	lambda, rho = generateRandomSamples(numSamples)

	constraint_values = evaluateConstraint(x,lambda,rho)

	abs_constraint_values = abs.(constraint_values)

	scalings::Float64 = max(median(abs_constraint_values),scaling_tol)

	return scalings
end

# estimate weak convexity parameter of the approximation
function estimateConstraintLipschitzConstant(x_ref::Array{Float64},numSamplesForLipVarEst::Int64,numGradientSamples::Int64,batchSize::Int64,mu::Float64,tau::Float64,scalings::Float64,objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	lambda, rho = generateRandomSamples(batchSize)

	sample_radius::Float64 = norm(x_ref)/10.0

	L_max::Float64 = 0.0
	L_max_sqr::Float64 = 0.0

	for batch = 1:batchSize
		xi_start::Int64 = batch
		xi_end::Int64 = batch

		L_max_xi::Float64 = 0.0
		for iter = 1:numSamplesForLipVarEst
			x_1 = pickPointOnSphere(x_ref,sample_radius)
			x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

			x_2 = pickPointOnSphere(x_ref,sample_radius)
			x_2,~ = projectOntoFeasibleSet(x_2,objectiveBound,proj_params)

			grad_1 = getConstraintStochasticGradient(x_1,lambda[:,xi_start:xi_end],rho[:,xi_start:xi_end],mu,tau,scalings)
			grad_2 = getConstraintStochasticGradient(x_2,lambda[:,xi_start:xi_end],rho[:,xi_start:xi_end],mu,tau,scalings)

			L_est::Float64 = norm(grad_2 - grad_1)/norm(x_2 - x_1)
			if(L_est > L_max_xi)
				L_max_xi = L_est
			end
		end

		L_max += L_max_xi/batchSize
		L_max_sqr += (L_max_xi)^2/batchSize
	end

	return L_max, L_max_sqr
end

# estimate step length factor
function estimateStepLength(x_ref::Array{Float64},numGradientSamples::Int64,numSamplesForLipVarEst::Int64,
							batchSize::Int64,maxNumIterations::Int64,minNumReplicates::Int64,mu::Float64,tau::Float64,
							scalings::Float64,objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	sample_radius::Float64 = norm(x_ref)/10.0

	grad_phi_up::Float64 = tau/4.0
	L_grad_phi::Float64 = 0.1*tau^2

	variance_phi = tau^2/16.0
	variance_avg::Float64 = 0.0
	for iter = 1:numSamplesForLipVarEst
		x_1 = pickPointOnSphere(x_ref,sample_radius)
		x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

		lambda, rho = generateRandomSamples(numGradientSamples*batchSize)

		variance_1::Float64 = 0.0

		for i = 1:batchSize
			xi_start::Int64 = (i-1)*numGradientSamples + 1
			xi_end::Int64 = i*numGradientSamples

			grad_1 = getConstraintStochasticGradient(x_1,lambda[:,xi_start:xi_end],rho[:,xi_start:xi_end],mu,tau,scalings)

			variance_1 += norm(grad_1)^2/batchSize
		end

		if(variance_1 > variance_avg)
			variance_avg = variance_1
		end
	end

	variance_est::Float64 = variance_avg*variance_phi

	L_max::Float64, L_max_sqr::Float64 = estimateConstraintLipschitzConstant(x_ref,numSamplesForLipVarEst,numGradientSamples,batchSize,mu,tau,scalings,objectiveBound)

	rho_est::Float64 = grad_phi_up*L_max + L_grad_phi*(2*sample_radius)^2*L_max_sqr/4.0 + L_grad_phi*L_max_sqr

	gamma::Float64 = 1/sqrt(rho_est*variance_est)

	stepLength::Float64 = gamma/sqrt((maxNumIterations+1.0)*minNumReplicates)

	return stepLength
end

# given a reference point, pick a random point within a given radius
function pickPointOnSphere(x_orig::Array{Float64},radius::Float64)

	if(radius < 1E-09)
		error("Sampling radius in pickPointOnSphere is almost zero! Set manually.")
	end

	n::Int64 = size(x_orig,1)
	unifRand = rand(1)
	factor1::Float64 = radius*(unifRand[1])^(1.0/n)
	normRand = rand(Normal(0,1), n)
	factor2::Float64 = norm(normRand)

	x_samp = x_orig + factor1*normRand/factor2

	return x_samp
end

ext: jl

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
module StopIon
using Unitful
using SpecialFunctions: erf

export Particle, bethe_bloch, fixed_background, self_consistent

const mc² = Unitful.me*Unitful.c0*Unitful.c0
const e⁴_mc² = (Unitful.q^2/(4π*Unitful.ɛ0))^2/mc²
const mc³_e² = mc²*1.0u"c"*4π*Unitful.ɛ0/(Unitful.q^2)
const mc²_ħ = mc²/Unitful.ħ

@derived_dimension Concentration inv(Unitful.𝐋^3)
@derived_dimension ParticleDensityFlux inv(Unitful.𝐋^2)

const shell_correction_coefficient = 0.5

struct Particle
    z::Float64
    m::Unitful.Mass
end

function γ²(ɛ::Unitful.Energy, p::Particle)
    return (1 + ɛ/(p.m*1u"c^2"))^2
end

function β²(ɛ::Unitful.Energy, p::Particle)
    return (γ²(ɛ,p)-1)/γ²(ɛ,p)
end

function MIP(z) :: Unitful.Energy
    return 13.6u"eV"*z
end

function EffectiveMIP(Z, Z_eff) :: Unitful.Energy
    return 13.6u"eV"*Z^2/(Z-Z_eff)
end

function z_p(ɛ::Unitful.Energy, p::Particle)
    return p.z*(1 - 1.034*exp(-137.04*sqrt(β²(ɛ, p))*p.z^(-0.688)))
end

function G(ɛ::Unitful.Energy, p::Particle, Te::Unitful.Energy)
    x = sqrt(mc²*β²(ɛ, p)*0.5/Te)
    return erf(ustrip(uconvert(Unitful.NoUnits, x))) - 2*x/sqrt(π)*exp(-x*x)
end

function L(ɛ::Unitful.Energy, p::Particle, target::Particle)
    return log(2*mc²/MIP(target.z)*(γ²(ɛ,p)-1)) - β²(ɛ,p) - shell_correction_coefficient/target.z
end

function L_bound(ɛ::Unitful.Energy, p::Particle, target::Particle, Z_eff)
    return log(2*mc²/EffectiveMIP(target.z, Z_eff)*(γ²(ɛ,p)-1)) - β²(ɛ,p) - shell_correction_coefficient/(target.z - Z_eff)
end

function L_free(ɛ::Unitful.Energy, p::Particle, N_a::Concentration, Te::Unitful.Energy, Z_eff)
    ω_p = sqrt(Z_eff*N_a*Unitful.q^2/(Unitful.ɛ0*Unitful.me))
    u = sqrt(2*Te/mc² + β²(ɛ,p))
    b_min = max(
        z_p(ɛ,p)/(mc³_e²*u^2),
        0.5/(mc²_ħ*u)
    )
    #b_min = z_p(ɛ,p)/(mc³_e²*u^2)
    η = (b_min*ω_p)/sqrt(β²(ɛ,p))
    return 0.764*log((1.0+η)/η)
end

function bethe_bloch(
    projectile::Particle,
    target::Particle,
    density::Unitful.Density,
    energy_samplings::Vector{EnergyT},
    energy_distribution_function,
    x_samplings::Vector{LengthT}) where EnergyT<:Unitful.Energy where LengthT<:Unitful.Length

    @assert length(energy_samplings) == length(energy_distribution_function) "energy samplings and energy distribution function must be of the same length!"

    KNaZ = 4π*(projectile.z)^2*e⁴_mc²*density/target.m*target.z

    x_diff = diff(x_samplings)
    deposition = fill(0.0u"MeV", length(x_diff))

    for (ɛ0, f) in zip(energy_samplings, energy_distribution_function)
        ɛ = ɛ0
        for (i, (_, Δx)) in enumerate(zip(deposition, x_diff))
            if (ɛ < 0u"MeV")
                break
            end
            if (L(ɛ, projectile, target) < 0)
                break
            end
            Δɛ = KNaZ/β²(ɛ, projectile)*L(ɛ, projectile, target)*Δx
            deposition[i] += Δɛ*f
            ɛ -= Δɛ
        end
    end

    return deposition
end

function fixed_background(
    projectile::Particle,
    target::Particle,
    density::Unitful.Density,
    electron_temperature::Unitful.Energy,
    effective_charge,
    energy_samplings::Vector{EnergyT},
    energy_distribution_function,
    x_samplings::Vector{LengthT}) where EnergyT<:Unitful.Energy where LengthT<:Unitful.Length

    @assert length(energy_samplings) == length(energy_distribution_function) "energy samplings and energy distribution function must be of the same length!"
    @assert effective_charge ≤ target.z "effective charge must not be greater than target charge!"

    KNa = 4π*e⁴_mc²*density/target.m
    KNaZ = KNa*(target.z-effective_charge)
    KNaZ_eff = KNa*effective_charge

    x_diff = diff(x_samplings)
    deposition = fill(0.0u"MeV", length(x_diff))

    for (ɛ0, f) in zip(energy_samplings, energy_distribution_function)
        ɛ = ɛ0
        for (i, (_, Δx)) in enumerate(zip(deposition, x_diff))
            if (ɛ < 0u"MeV")
                break
            end
            Δɛ_free = KNaZ_eff*z_p(ɛ,projectile)^2/β²(ɛ, projectile)*G(ɛ, projectile, electron_temperature)*L_free(ɛ, projectile, density/target.m, electron_temperature, effective_charge)*Δx
            Δɛ_bound = KNaZ*z_p(ɛ,projectile)^2/β²(ɛ, projectile)*L_bound(ɛ, projectile, target, effective_charge)*Δx
            Δɛ = Δɛ_free+Δɛ_bound
            deposition[i] += Δɛ*f
            ɛ -= Δɛ
        end
    end

    return deposition
end

function self_consistent(
    projectile::Particle,
    particle_density_flux::ParticleDensityFlux,
    target::Particle,
    density::Unitful.Density,
    energy_samplings::Vector{EnergyT},
    energy_distribution_function,
    x_samplings::Vector{LengthT}) where EnergyT<:Unitful.Energy where LengthT<:Unitful.Length

    @assert length(energy_samplings) == length(energy_distribution_function) "energy samplings and energy distribution function must be of the same length!"

    Na = density/target.m
    KNa = 4π*e⁴_mc²*Na

    x_diff = diff(x_samplings)
    deposition = fill(0.0u"MJ/g", length(x_diff))
    T_e = fill(0.01u"eV", length(x_diff))
    Z_eff = fill(0.01, length(x_diff))

    for (ɛ0, f) in zip(reverse(energy_samplings), reverse(energy_distribution_function))
        ɛ = ɛ0
        for (i, Δx) in enumerate(x_diff)
            if (ɛ < 0u"MeV")
                break
            end

            KNaZ = KNa*(target.z-Z_eff[i])
            KNaZ_eff = KNa*Z_eff[i]

            Δɛ = 0u"eV"

            Δɛ_free = KNaZ_eff*z_p(ɛ,projectile)^2/β²(ɛ, projectile)*G(ɛ, projectile, T_e[i])*L_free(ɛ, projectile, density/target.m, T_e[i], Z_eff[i])*Δx
            if Δɛ_free > 0u"MeV"
                T_e[i] += f*particle_density_flux*Δɛ_free/(Na*Δx*Z_eff[i])/1.5
                Δɛ += Δɛ_free
            end

            Z_eff_old = Z_eff[i]
            Δɛ_bound = KNaZ*z_p(ɛ,projectile)^2/β²(ɛ, projectile)*L_bound(ɛ, projectile, target, Z_eff[i])*Δx
            if (Δɛ_bound > 0u"MeV") && (Z_eff[i] < target.z)
                ionization_energy = f*particle_density_flux*Δɛ_bound/(Na*Δx)
                while ionization_energy > 0u"eV"
                    ionization_energy_limit = (floor(Z_eff[i]+1)-Z_eff[i])*MIP(floor(target.z-Z_eff[i]))*(target.z/(target.z-Z_eff[i]))^2
                    if ionization_energy > ionization_energy_limit
                        Z_eff[i] = floor(Z_eff[i]+1)
                        ionization_energy -= ionization_energy_limit
                    else
                        Z_eff[i] += ionization_energy/ionization_energy_limit * (floor(Z_eff[i]+1)-Z_eff[i])
                        break
                    end
                end

                T_e[i] *= Z_eff_old/Z_eff[i]

                Δɛ += Δɛ_bound
            end

            deposition[i] += Δɛ*f*particle_density_flux/(density*Δx)
            ɛ -= Δɛ
        end
    end

    return deposition, T_e, Z_eff
end

end


text:
export macau_hmc

using Compat

type HMCModel
  momentum::Matrix{Float64}
  G       ::Matrix{Float64} ## diagonal mass matrix
end

deepcopy(m::HMCModel) = HMCModel(copy(m.momentum))

function HMCModel(num_latent::Int, N::Int, Ldiag::Vector{Float64})
  return HMCModel(
    zeros(num_latent, N),  ## momentum
    repmat(Ldiag, 1, N)    ## mass matrix (G)
  )
end

function macau_hmc(data::RelationData;
                   num_latent::Int = 10,
                   verbose::Bool   = true,
                   burnin::Int     = 100,
                   psamples::Int   = 100,
                   L::Int          = 10,
                   L_inner::Int    = 1,
                   prior_freq::Int = 8,  ## how often to update prior
                   eps::Float64    = 0.01,
                   reset_model     = true,
                   clamp::Vector{Float64} = Float64[])
  ## initialization
  verbose     && println("Model setup")
  reset_model && reset!(data, num_latent)

  rmse = NaN
  rmse_sample = NaN
  rmse_train  = NaN
  Umodel = HMCModel(num_latent, data.entities[1].count, diag(data.entities[1].model.Lambda))
  Vmodel = HMCModel(num_latent, data.entities[2].count, diag(data.entities[2].model.Lambda))

  ## data
  df  = data.relations[1].data.df
  mean_value = mean(df[:,end])
  uid = convert(Vector{Int32}, df[:,1])
  vid = convert(Vector{Int32}, df[:,2])
  val = convert(Array, df[:,3]) - mean_value
  Udata = sparse(vid, uid, val, data.entities[2].count, data.entities[1].count)
  Vdata = Udata'

  test_uid = convert(Vector{Int}, data.relations[1].test_vec[:,1])
  test_vid = convert(Vector{Int}, data.relations[1].test_vec[:,2])
  test_val = convert(Vector, data.relations[1].test_vec[:,3])
  test_idx = hcat(test_uid, test_vid)

  alpha = data.relations[1].model.alpha
  yhat_post = zeros(length(test_val))

  for i in 1 : burnin + psamples
    time0 = time()

    if i == burnin + 1
      verbose && print("================== Burnin complete ===================\n")
    end

    verbose && @printf("======= Step %d =======\n", i)
    verbose && @printf("eps = %.2e\n", eps)

    # HMC sampling momentum
    sample!(Umodel)
    sample!(Vmodel)
    kinetic_start   = computeKinetic(Umodel) + computeKinetic(Vmodel)
    potential_start = computePotential(uid, vid, val, data.relations[1])
    Ustart = copy(data.entities[1].model.sample)
    Vstart = copy(data.entities[2].model.sample)

    # Follow Hamiltonians
    hmc_update_u!(Umodel, Vmodel, Udata, alpha, data.entities[1], data.entities[2], L_inner, eps / 2)
    for l in 1:L
      hmc_update_u!(Vmodel, Umodel, Vdata, alpha, data.entities[2], data.entities[1], L_inner, eps)
      if l < L
        hmc_update_u!(Umodel, Vmodel, Udata, alpha, data.entities[1], data.entities[2], L_inner, eps)
      end
      verbose && @printf("  Momentum %d: |r_U| = %.4e, |r_V| = %.4e\n", l, vecnorm(Umodel.momentum), vecnorm(Vmodel.momentum))
    end
    hmc_update_u!(Umodel, Vmodel, Udata, alpha, data.entities[1], data.entities[2], L_inner, eps / 2)
    verbose && @printf("  Momentum L: |r_U| = %.4e, |r_V| = %.4e\n", vecnorm(Umodel.momentum), vecnorm(Vmodel.momentum))

    # Verify p(U,V) is fine
    kinetic_final   = computeKinetic(Umodel) + computeKinetic(Vmodel)
    potential_final = computePotential(uid, vid, val, data.relations[1])

    dH = potential_start - potential_final + kinetic_start - kinetic_final
    verbose && @printf("  ΔH = %.4e  ΔKin = %.4e  ΔPot = %.4e\n",
                       -dH,
                       kinetic_final-kinetic_start, potential_final - potential_start)
    if rand() < exp(dH)
      ## accept
      verbose && print("-> ACCEPTED!\n")
    else
      ## reject
      verbose && print("-> REJECTED!\n")
      Base.copy!(data.entities[1].model.sample, Ustart)
      Base.copy!(data.entities[2].model.sample, Vstart)
      if dH < -6
        ## decreasing eps by 2
        neweps = eps / 2
        newL   = ceil(Int, L*1.6)
        verbose && @printf("Reducing eps from %.2e to %.2e.\n", eps, neweps)
        verbose && @printf("Increasing L from %d to %d.\n", L, newL)
        eps = neweps
        L   = newL
      end
    end

    # Sampling prior for latents
    if i % prior_freq == 0
      verbose && println("Updating priors...")
      for en in data.entities
        update_latent_prior!(en, true)
      end
    end

    rel  = data.relations[1]
    yhat_raw = pred(rel, test_idx, rel.test_F)
    yhat     = clamp!(yhat_raw, clamp)
    update_yhat_post!(yhat_post, yhat_raw, i, burnin)

    rmse = sqrt( mean((yhat - test_val).^2) )
    rmse_post = sqrt( mean((makeClamped(yhat_post, clamp) - test_val).^2) )
    #yhat_train = clamp!(pred(Umodel, Vmodel, mean_value, uid, vid), clamp)
    #rmse_train = sqrt( mean((yhat_train - mean_value - val).^2) )
    time1 = time()

    if verbose
      @printf("% 3d: |U|=%.4e  |V|=%.4e  RMSE=%.4f  RMSE(avg)=%.4f [took %.2fs]\n",
        i,
        vecnorm(data.entities[1].model.sample),
        vecnorm(data.entities[2].model.sample),
        rmse,
        rmse_post,
        time1 - time0)
    end
  end

  return @compat Dict(
    "rmse"   => rmse,
    "rmse_train" => rmse_train,
    "alpha"  => alpha)
end

function sample!(m::HMCModel)
  ## sampling from Normal distribution with identity cov
  for i = 1:size(m.momentum,1), j = 1:size(m.momentum,2)
    m.momentum[i,j] = randn() / sqrt(m.G[i,j])
  end
  nothing
end

function hmc_update_u!(U_hmcmodel::HMCModel,
                       V_hmcModel::HMCModel,
                       Udata,
                       alpha::Float64,
                       en::Entity,
                       enV::Entity,
                       L::Int,
                       eps::Float64)

  model      = en.model
  sample     = model.sample
  Vsample    = enV.model.sample
  momentum   = U_hmcmodel.momentum
  num_latent = size(U_hmcmodel.momentum, 1)

  # 1) make half step in momentum
  subtract_grad!(momentum, model, sample, Vsample, Udata, alpha, eps / 2)

  # 2) L times: make full steps with U, full step with momentum (except for last step)
  for i in 1:L
    add!(sample, momentum, eps)
    if i < L
      subtract_grad!(momentum, model, sample, Vsample, Udata, alpha, eps)
    end
  end

  # 3) make half step in momentum
  subtract_grad!(momentum, model, sample, Vsample, Udata, alpha, eps / 2)
  nothing
end

## computes X += mult * Y
function add!(X::Matrix, Y::Matrix, mult::Float64)
  size(X) == size(Y) || error("X and Y must have the same size.")
  @inbounds @simd for i in 1:length(X)
    X[i] += mult * Y[i]
  end
  nothing
end

function subtract_grad!(momentum::Matrix{Float64},
                        model   ::EntityModel,
                        sample  ::Matrix{Float64},
                        Vsample ::Matrix{Float64},
                        Udata   ::SparseMatrixCSC,
                        alpha   ::Float64,
                        eps     ::Float64)
  num_latent = size(momentum, 1)

  for n in 1:size(sample, 2)
    tmp = grad(n, sample, Vsample, Udata, model.Lambda, model.mu, alpha)
    @inbounds for k in 1:num_latent
      momentum[k,n] -= eps * tmp[k]
    end
  end
  nothing
end

function grad(n,
              sample  ::Matrix{Float64},
              Vsample ::Matrix{Float64},
              Udata   ::SparseMatrixCSC,
              Lambda  ::Matrix{Float64},
              mu      ::Vector{Float64},
              alpha   ::Float64)
  un  = sample[:, n]
  idx = Udata.colptr[n] : Udata.colptr[n+1]-1
  ff  = Udata.rowval[ idx ]
  rr  = Udata.nzval[ idx ]

  MM = Vsample[:, ff]
  ## add d/du A(r_V | U)
  ## add d/du 1/2 log|G_V|

  return - alpha * (MM * rr - (MM*MM') * un) - Lambda * (mu - un)
end

function computeKinetic(m::HMCModel)
  kin = 0.0
  momentum = m.momentum
  G        = m.G
  ## r'*G*r + log|G|
  @inbounds @simd for i = 1:length(momentum)
    kin += momentum[i] * momentum[i] * G[i] + log(G[i])
  end
  return 0.5 * kin
end

## computes - log [p(Y | U,V) p(U, V)]
function computePotential(uid::Vector, vid::Vector, val::Vector, r::Relation)
  alpha::Float64 = r.model.alpha
  Umodel = r.entities[1].model
  Vmodel = r.entities[2].model
  Usample = Umodel.sample
  Vsample = Vmodel.sample
  energy  = 0.0

  ## data
  for i in 1:length(uid)
    energy += (column_dot(Usample, Vsample, uid[i], vid[i]) - val[i]) ^ 2
  end
  energy *= alpha / 2

  ## Priors
  ## sum(u_i' * Lambda_u * u_i) = trace(Lambda_u * U * U')
  energy += sum(Umodel.Lambda .* (Usample * Usample')) / 2
  energy += sum(Vmodel.Lambda .* (Vsample * Vsample')) / 2

  ## -2 * sum(mu_u' * Lambda_u * u_i)
  energy -= first(Umodel.mu' * Umodel.Lambda * sum(Usample, 2))
  energy -= first(Vmodel.mu' * Vmodel.Lambda * sum(Vsample, 2))

  return energy
end

function column_dot(X, Y, i, j)
  size(X,1) == size(Y,1) || error("X and Y must have the same number of rows.")
  size(X,2) >= i         || error("X must have at least i columns.")
  size(Y,2) >= j         || error("Y must have at least j columns.")
  d = 0.0
  @inbounds for k = 1:size(X,1)
    d += X[k,i] * Y[k,j]
  end
  return d
end

function update_yhat_post!(yhat_post, yhat_raw, i, burnin)
  ## also takes care of the first sample of posterior
  if i <= burnin + 1
    Base.copy!(yhat_post, yhat_raw)
    return yhat_post
  end

  ## averaging
  n = i - burnin - 1
  Base.copy!(yhat_post, (n*yhat_post + yhat_raw) / (n + 1))
  return yhat_post
end


text:
using ITensors,
      Test

@testset "emptyITensor (Empty)" begin

  @testset "emptyITensor set elements" begin
    i = Index(2; tags="i")

    E = emptyITensor(i', dag(i))

    @test hassameinds(E, (i', i))
    @test order(E) == 2
    @test_broken E[i' => 1, i => 1] == 0

    E[i' => 1, i => 2] = 2.3

    @test E[i' => 1, i => 1] == 0
    @test E[i' => 2, i => 1] == 0
    @test E[i' => 1, i => 2] == 2.3
    @test E[i' => 2, i => 2] == 0
  end

  @testset "emptyITensor set elements (QN)" begin
    i = Index(QN(0) => 2,
              QN(1) => 2; tags="i")

    E = emptyITensor(i', dag(i))

    @test hassameinds(E, (i', i))
    @test order(E) == 2
    @test isnothing(flux(E))
    @test_broken E[i' => 1, i => 3] == 0

    E[i' => 3, i => 2] = 2.3

    @test flux(E) == QN(1)

    @test E[i' => 1, i => 1] == 0
    @test E[i' => 2, i => 1] == 0
    @test E[i' => 3, i => 2] == 2.3
    @test E[i' => 2, i => 3] == 0
    @test_throws ErrorException E[i' => 2, i => 3] = 3.2
  end

  @testset "emptyITensor(Any)" begin
    i = Index(QN(0) => 2,
              QN(1) => 2; tags="i")

    E = emptyITensor(Any)

    @test isnothing(flux(E))
    @test order(E) == Any
    @test_throws ErrorException E[i' => 1, i => 3] = 0

    A = randomITensor(i', dag(i))
    E += A

    @test E == A
  end
end

nothing


text:
function compute_tax(tax_parameters)
"""
Calculate optimal tax
"""

    space = tax_parameters.space
    c_store = tax_parameters.coefficients
    mparams = tax_parameters.mparams
    opts = tax_parameters.opts
    states = tax_parameters.states
    exog_states = tax_parameters.exog_states
    smax = tax_parameters.smax
    smin = tax_parameters.smin
    loc_next = tax_parameters.loc_next
    dam_next = tax_parameters.dam_next
    controls = tax_parameters.controls

    # Vectors denoting derivatives to take with respect to each state
    if opts[:learn_e]
        k_deriv = [1 0 0 0 0]
        m_deriv = [0 1 0 0 0]
        mu_deriv = [0 0 1 0 0]
        mu_deriv_2 = [0 0 1 0 0]
        sigma_deriv = [0 0 0 1 0]
        l_deriv = [0 0 0 0 1]
        mll_deriv = [0 1 0 0 2]
    else
        k_deriv = [1 0 0]
        m_deriv = [0 1 0]
        l_deriv = [0 0 1]
        mll_deriv = [0 1 2]
    end

    # d cumulative emissions/d emission
    dm_de = 1.

    # Reshaped quadrature weights for expectations over vectors
    mparams[:weights2] = mparams[:weights]*mparams[:weights]'
    mparams[:weights2] = reshape(mparams[:weights2], mparams[:nqnodes]^2)

    # Reshaped quadrature for damage and prior_loc
    dam_next = [reshape(dam_next[i], mparams[:nqnodes]^2) for i = 1:length(dam_next)]
    loc_next = [reshape(loc_next[i], mparams[:nqnodes]^2) for i = 1:length(loc_next)]

    # Expected damage and prior_loc next period
    dam_next_exp = [(mparams[:weights2]'*dam_next[i])[1] for i = 1:length(dam_next)]
    loc_next_exp = [(mparams[:weights2]'*loc_next[i])[1] for i = 1:length(loc_next)]

    # Variance of damage and prior_loc next period, why is it zero for first period?
    dam_var = [(mparams[:weights2]'*dam_next[i].^2 - (mparams[:weights2]'*dam_next[i]).^2)[1] for i = 1:length(dam_next)]
    loc_var = [(mparams[:weights2]'*loc_next[i].^2 - (mparams[:weights2]'*loc_next[i]).^2)[1] for i = 1:length(loc_next)]

    # Covariance between fractional net output and d_2 mean
    dl_covar = [(mparams[:weights2]'*(dam_next[i].*loc_next[i])-
        (dam_next_exp[i].*loc_next_exp[i]))[1] for i = 1:length(loc_next)]

    # States for simulation
    s_base = Array{Array}(undef, size(states, 1) - 1)
    s_base_exp = Array{Array}(undef, size(states, 1) - 1)
    s_base_ce = Array{Array}(undef, size(states, 1) - 1)

    # Compute next period's states, expected states, and certainty states
    if opts[:learn_e]
        for t = 2:size(states, 1)
            s_base[t - 1] = [states[t,[1, 2]]'.*ones(mparams[:nqnodes]^2, 2) loc_next[t - 1] states[t, 4].*ones(mparams[:nqnodes]^2, 1) dam_next[t - 1]]
            s_base_exp[t - 1] = [states[t,[1, 2]]' loc_next_exp[t - 1] states[t, 4] dam_next_exp[t - 1]]
            s_base_ce[t - 1] = [states[t,[1, 2]]' states[t, 3] 0 1 ./ (1+mparams[:coeff_mean].*(mparams[:ccr]*states[t, 2]).^states[t, 3])]
        end
    else
        for t = 2:size(states, 1)
            s_base[t - 1] = [states[t,[1, 2]]'.*ones(mparams[:nqnodes]^2, 2)  dam_next[t - 1]]
            s_base_exp[t - 1] = [states[t,[1, 2]]' dam_next_exp[t - 1]]
            s_base_ce[t - 1] = [states[t,[1, 2]]' 1 ./ (1+mparams[:coeff_mean].*(mparams[:ccr]*states[t, 2]).^states[t, 3])]
        end
    end

    # Initialize vectors for value partials and channels
    Vstar_beta = Vector{Array}(undef, size(states, 1) - 1)
    dvstar_dk = Vector{Array}(undef, size(states, 1) - 1)
    dvstar_dm = Vector{Array}(undef, size(states, 1) - 1)
    dvstar_dl = Vector{Array}(undef, size(states, 1) - 1)
    dvstar_dm_ce = Vector{Float64}(undef, size(states, 1) - 1)
    dvstar_dl_ce = Vector{Float64}(undef, size(states, 1) - 1)
    dvstar_dm_exp = Vector{Float64}(undef, size(states, 1) - 1)
    dvstar_dl_exp = Vector{Float64}(undef, size(states, 1) - 1)
    dvstar_dl3_exp = Vector{Float64}(undef, size(states, 1) - 1)
    dvstar_dl2_dm_exp = Vector{Float64}(undef, size(states, 1) - 1)
    if opts[:learn_e]
        dvstar_dsigma_ce = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dsigma_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dsigma = Vector{Array}(undef, size(states, 1) - 1)
        dvstar_dexp = Vector{Array}(undef, size(states, 1) - 1)
        dvstar_dl2_dsigma_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dl_dm_dmu_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dl_dsigma_dmu_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dl2_dmu_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dexp2_dm_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dexp2_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dexp2_dsigma_exp = Vector{Float64}(undef, size(states, 1) - 1)
        dvstar_dexp2_dl_exp = Vector{Float64}(undef, size(states, 1) - 1)
    end

    # Compute values and value partials for channels
    for t = 1:size(states, 1) - 1

        if opts[:adapted]
            space = fundefn(:cheb, mparams[:approximation_level],smin[t + 1, :],smax[t + 1, :])
        end

        # Set year
        space[:year] = t - 1

        # Set value function coefficient vector
        space[:c_cont] = c_store[t + 1]

        # beta_t*V for misspecification insurance
        Vstar_beta[t] = exog_states[:betaeff][t].*reshape(funeval(c_store[t + 1], space, s_base[t])[1], mparams[:nqnodes], mparams[:nqnodes])

        # dv/dcapital
        dvstar_dk[t] = funeval(c_store[t + 1], space, s_base[t], k_deriv)[1]
        # dv/dcumulative emissions
        dvstar_dm[t] = funeval(c_store[t + 1], space, s_base[t],  m_deriv)[1]
        # dv/dfractional net output
        dvstar_dl[t] = funeval(c_store[t + 1], space, s_base[t], l_deriv)[1]

        if opts[:learn_e]
            # dv/dmean
            dvstar_dexp[t] = funeval(c_store[t + 1], space, s_base[t],  mu_deriv)[1]
            # dv/dvariance
            dvstar_dsigma[t] = funeval(c_store[t + 1], space, s_base[t], sigma_deriv)[1]
        end

        ## Decomposition components

        # Cumulative emissions
        # Certainty-equivalent
        dvstar_dm_ce[t] = funeval(c_store[t + 1], space, s_base_ce[t], m_deriv)[1][1]
        # At expected state
        dvstar_dm_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], m_deriv)[1][1]
        # Precaution: CO2 interaction
        dvstar_dl2_dm_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], mll_deriv)[1][1]

        if opts[:learn_e]
            # Signal smoothing
            dvstar_dexp2_dm_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 1, 2, 0, 0])[1][1]
            dvstar_dexp2_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 0, 2, 0, 0])[1][1]
        end

        if opts[:learn_e]
            # Variance
            # Certainty-equivalent
            dvstar_dsigma_ce[t] = funeval(c_store[t + 1], space, s_base_ce[t], sigma_deriv)[1][1]
            # At expected state
            dvstar_dsigma_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], sigma_deriv)[1][1]
            # Precaution: scale interaction
            dvstar_dl2_dsigma_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 0, 0, 1, 2])[1][1]
            # Signal smoothing
            dvstar_dl_dsigma_dmu_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 0, 1, 1, 1])[1][1]
            # Signal smoothing
            dvstar_dexp2_dsigma_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 0, 2, 1, 0])[1][1]
        end

        # Fractional net output
        # Certainty-equivalent
        dvstar_dl_ce[t] = funeval(c_store[t + 1], space, s_base_ce[t], l_deriv)[1][1]
        # At Expected-state
        dvstar_dl_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], l_deriv)[1][1]
        # Precaution
        dvstar_dl3_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], l_deriv*3)[1][1]

        # Precaution: belief interaction
        if opts[:learn_e]
            dvstar_dl2_dmu_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 0, 1, 0, 2])[1][1]
            dvstar_dl_dm_dmu_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 1, 1, 0, 1])[1][1]
        end

        # Signal smoothing
        if opts[:learn_e]
            dvstar_dexp2_dl_exp[t] = funeval(c_store[t + 1], space, s_base_exp[t], [0, 0, 2, 0, 1])[1][1]
        end
    end

    # Fractional net output next period quadrature
    dam_next = [reshape(dam_next[i], mparams[:nqnodes], mparams[:nqnodes]) for i = 1:length(dam_next)]

    # dfractional net output/de quadrature
    if opts[:learn_e]
        dl_dm = [-qnwnorm(mparams[:nqnodes], states[i, 3], states[i, 4])[1]./
            s_base[i][1, 2].*(dam_next[i] - dam_next[i].^2) for i = 1:length(s_base)]
    else
        dl_dm = [-qnwnorm(mparams[:nqnodes], mparams[:exp_mean], mparams[:exp_var])[1]./
            s_base[i][1, 2].*(dam_next[i] - dam_next[i].^2) for i = 1:length(s_base)]
    end

    if opts[:learn_e]
        # dfractional net output/de exp
        dl_dm_exp = [(-mparams[:weights]'*(qnwnorm(mparams[:nqnodes], states[i, 3], states[i, 4])[1]./
            s_base[i][1, 2].*(dam_next[i] - dam_next[i].^2))*mparams[:weights])[1] for i = 1:length(s_base)]

        # dfractional net output/de ce
        dl_dm_ce =
            -[s_base_ce[i][3]/s_base_ce[i][2]*(s_base_ce[i][5] - s_base_ce[i][5]^2) for i = 1:length(s_base)]
    else
        # dfractional net output/de exp
        dl_dm_exp = [(-mparams[:weights]'*(qnwnorm(mparams[:nqnodes], mparams[:exp_mean], mparams[:exp_var])[1]./
            s_base[i][1, 2].*(dam_next[i] - dam_next[i].^2))*mparams[:weights])[1] for i = 1:length(s_base)]

        # dfractional net output/de ce
        dl_dm_ce =
            -[mparams[:exp_mean]./s_base_ce[i][2]*(s_base_ce[i][3] - s_base_ce[i][3]^2) for i = 1:length(s_base)]

    end

    # dfractional net output/de
    dl_dm = [reshape(dl_dm[i], mparams[:nqnodes]^2) for i = 1:length(dl_dm)]

    dam_next = [reshape(dam_next[i], mparams[:nqnodes]^2) for i = 1:length(dam_next)]

    # dvariance/de
    if opts[:learn_e]
        dsigma_de =
            [-2.0.*log.(mparams[:ccr]*s_base[i][1, 2])*states[i, 4]*(states[i, 4]*mparams[:coeff_scale])./
            s_base[i][1, 2]./(mparams[:coeff_scale] .+ states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])^2).^2 for i = 1:length(s_base_ce)]
    end

    if opts[:learn_e]
        # dmean/de
        dmu_dm =
            [(((states[i, 4]*(1 ./ s_base[i][1, 2])*(log.(1 ./ dam_next[i] .- 1) .- mparams[:coeff_loc]) .+
            states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])*(1 ./ (1 ./ dam_next[i] .- 1)).*( -1 ./ (dam_next[i].^2)).*dl_dm[i]).*
            (mparams[:coeff_scale] .+ states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2]).^2)) -
            (2*states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])./s_base[i][1, 2]*
            (mparams[:coeff_scale]*states[i, 3] .+
            states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])*(log.(1 ./ dam_next[i].- 1) .- mparams[:coeff_loc]))))./
            (mparams[:coeff_scale] .+ states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])^2).^2 for i = 1:length(s_base_ce)]

        # dmean/de exp
        dmu_dm_exp =
            [mparams[:weights2]'*(((states[i, 4]*(1 ./ s_base[i][1, 2])*(log.(1 ./ dam_next[i] .- 1) .- mparams[:coeff_loc]) .+
            states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])*(1 ./ (1 ./ dam_next[i] .- 1)).*( -1 ./ (dam_next[i].^2)).*dl_dm[i]).*
            (mparams[:coeff_scale] .+ states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2]).^2)) -
            (2*states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])./s_base[i][1, 2]*
            (mparams[:coeff_scale]*states[i, 3] .+
            states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])*(log.(1 ./ dam_next[i] .- 1) .- mparams[:coeff_loc]))))./
            (mparams[:coeff_scale] .+ states[i, 4]*log.(mparams[:ccr]*s_base[i][1, 2])^2).^2 for i = 1:length(s_base_ce)]

        # dmu/de
        dmu_dm = [reshape(dmu_dm[i], mparams[:nqnodes]^2) for i = 1:length(dmu_dm)]
    end

    # Insurance covariance term cov(dv_dl, dl_de)
    insurance_l =
        -[(mparams[:weights2]'*((dvstar_dl[i].-mparams[:weights2]'*dvstar_dl[i])
        .*(dl_dm[i].-dl_dm_exp[i])))[1] for i = 1:length(dvstar_dl)]

    # Insurance covariance term cov(dv_dmu, dmu_de)
    if opts[:learn_e]
        insurance_mu =
            -[(mparams[:weights2]'*((dvstar_dexp[i].-mparams[:weights2]'*dvstar_dexp[i])
            .*(dmu_dm[i].-dmu_dm_exp[i])))[1] for i = 1:length(dvstar_dl)]
    end

    # Define channels dictionary
    channels = Dict{Symbol,Any}()

    # Marginal utility function
    MU(x,t) = x.^(mparams[:rho] - 1)/mparams[:scaling_factor]

    # NOTE FOR ALL CHANNELS BELOW:
    # 3/11 translates tons of carbon to tons of CO2
    # 1000 translates billions to millions between capital and labor
    # At.*Lt adjusts effective terms into real terms
    if opts[:learn_e]
        channels[:ce_tax] =
            [-(1000*3/11*(dvstar_dm_ce[i]
            .+ dvstar_dsigma_ce[i]*dsigma_de[i]
            .+ dvstar_dl_ce[i]*dl_dm_ce[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:exp_tax] =
            [-(1000*3/11*(dvstar_dm_exp[i]
            .+ dvstar_dsigma_exp[i]*dsigma_de[i] .+
            dvstar_dl_exp[i]*dl_dm_exp[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:unc_adjust] =
            [-(1000*3/11*((dvstar_dm_exp[i] - dvstar_dm_ce[i])*dm_de
            .+ (dvstar_dsigma_exp[i] - dvstar_dsigma_ce[i])*dsigma_de[i]
            .+ (dvstar_dl_exp[i]*dl_dm_exp[i]- dvstar_dl_ce[i]*dl_dm_ce[i]))
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:output_insurance] =
            [(1000*3/11*(insurance_l[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:learning_insurance] =
            [(1000*3/11*(insurance_mu[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:precaution] =
            [-(1000*3/11*(1/2*(dvstar_dl3_exp[i]*dl_dm_exp[i] .+ dvstar_dl2_dm_exp[i] .+ dvstar_dl2_dsigma_exp*dsigma_de[i])*dam_var[i]
            .+ (dvstar_dl2_dmu_exp[i]*dl_dm_exp[i] .+ dvstar_dl_dm_dmu_exp[i] .+ dvstar_dl_dsigma_dmu_exp*dsigma_de[i])*dl_covar[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:smoothing] =
            [-(1000*3/11*(1/2*(dvstar_dexp2_dl_exp[i]*dl_dm_exp[i] .+ dvstar_dexp2_dm_exp[i] .+ dvstar_dexp2_dsigma_exp[i]*dsigma_de[i])*loc_var[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

    else

        channels[:ce_tax] =
            [-(1000*3/11*(dvstar_dm_ce[i]
            .+ dvstar_dl_ce[i]*dl_dm_ce[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:exp_tax] =
            [-(1000*3/11*(dvstar_dm_exp[i]
            .+ dvstar_dl_exp[i]*dl_dm_exp[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:unc_adjust] =
            [-(1000*3/11*((dvstar_dm_exp[i] - dvstar_dm_ce[i])*dm_de
            .+ (dvstar_dl_exp[i]*dl_dm_exp[i]- dvstar_dl_ce[i]*dl_dm_ce[i]))
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:output_insurance] =
            [(1000*3/11*(insurance_l[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:learning_insurance] = [0. for i = 1:length(dvstar_dk)]

        channels[:precaution] =
            [-(1000*3/11*(1/2*(dvstar_dl3_exp[i]*dl_dm_exp[i] .+ dvstar_dl2_dm_exp[i])*dam_var[i])
            .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
            exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        channels[:smoothing] = [0. for i = 1:length(dvstar_dk)]

    end

    # Compute misspecification insurance
    if opts[:robust]

        if opts[:learn_e]

            channels[:scc_direct] = [-(1000*3/11*mparams[:weights2]'*
                (dvstar_dm[i]
                .+ dvstar_dexp[i].*dmu_dm[i]
                .+ dvstar_dsigma[i].*dsigma_de[i]
                .+ dvstar_dl[i].*dl_dm[i])
                .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
                exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

            # dV/de with cols being coeff quad, rows being exponent quad
            deriv_cont_value_quad = [(dvstar_dm[i]
                + dvstar_dexp[i].*dmu_dm[i]
                + dvstar_dsigma[i].*dsigma_de[i]
                + dvstar_dl[i].*dl_dm[i]) for i = 1:length(dvstar_dk)]

            deriv_cont_value_quad =
                [reshape(deriv_cont_value_quad[i], mparams[:nqnodes], mparams[:nqnodes]) for i = 1:length(deriv_cont_value_quad)]

            # Taking expectation over model dimension: E_{d_1,omega}[dV/de]
            deriv_cont_value_exp = [(mparams[:weights]'*deriv_cont_value_quad[i]')' for i = 1:length(dvstar_dk)]

            # Shock to transition distribution from RC operator: exp(-beta*V/theta)
            rc_shock_quad = [exp.(-Vstar_beta[i]/mparams[:theta_robust]) for i = 1:length(dvstar_dk)]

            # Expectation over the model dimension: E_{d_1,omega}[exp(-beta*V/theta)]
            rc_shock_quad_exp = [(mparams[:weights]'*rc_shock_quad[i]')' for i = 1:length(dvstar_dk)]

            # E_{d_1,omega}[XY]
            rc_adj_val_shock = [deriv_cont_value_quad[i].*rc_shock_quad[i] for i = 1:length(dvstar_dk)]
            rc_adj_val_shock_exp = [(mparams[:weights]'*rc_adj_val_shock[i]')' for i = 1:length(dvstar_dk)]

            # E_{d_1,omega}[X]E_{d_1,omega}[Y]
            rc_adj_val_exp_shock_exp = [rc_shock_quad_exp[i].*deriv_cont_value_exp[i] for i = 1:length(dvstar_dk)]

            # E_{d_2}[E_{d_1,omega}[XY] - E_{d_1,omega}[X]E_{d_1,omega}[Y]] <=> E_{d_2}[cov_{d_1,omega}(X,Y)]
            rc_adjustment = [mparams[:weights]'*(rc_adj_val_shock_exp[i] - rc_adj_val_exp_shock_exp[i]) for i = 1:length(dvstar_dk)]

            channels[:robust_insurance] = [-(1000*3/11*rc_adjustment
                .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
                exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        else

            channels[:scc_direct] = [-(1000*3/11*mparams[:weights2]'*
                (dvstar_dm[i]
                .+ dvstar_dl[i].*dl_dm[i])
                .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
                exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

            # dV/de with cols being coeff quad, rows being exponent quad
            deriv_cont_value_quad = [(dvstar_dm[i]
                + dvstar_dl[i].*dl_dm[i]) for i = 1:length(dvstar_dk)]

            deriv_cont_value_quad =
            [reshape(deriv_cont_value_quad[i], mparams[:nqnodes], mparams[:nqnodes]) for i=1:length(deriv_cont_value_quad)]

            # Taking expectation over model dimension: E_{d_1,omega}[dV/de]
            deriv_cont_value_exp = [(mparams[:weights]'*deriv_cont_value_quad[i])' for i = 1:length(dvstar_dk)]

            # Shock to transition distribution from RC operator: exp(-beta*V/theta)
            rc_shock_quad = [exp.(-Vstar_beta[i]/mparams[:theta_robust]) for i = 1:length(dvstar_dk)]

            # Expectation over the model dimension: E_{d_1,omega}[exp(-beta*V/theta)]
            rc_shock_quad_exp = [(mparams[:weights]'*rc_shock_quad[i])' for i = 1:length(dvstar_dk)]

            # E_{d_1,omega}[XY]
            rc_adj_val_shock = [deriv_cont_value_quad[i].*rc_shock_quad[i] for i = 1:length(dvstar_dk)]
            rc_adj_val_shock_exp = [(mparams[:weights]'*rc_adj_val_shock[i])' for i = 1:length(dvstar_dk)]

            # E_{d_1,omega}[X]E_{d_1,omega}[Y]
            rc_adj_val_exp_shock_exp = [rc_shock_quad_exp[i].*deriv_cont_value_exp[i] for i = 1:length(dvstar_dk)]

            # E_{d_2}[E_{d_1,omega}[XY] - E_{d_1,omega}[X]E_{d_1,omega}[Y]] <=> E_{d_2}[cov_{d_1,omega}(X,Y)]
            rc_adjustment = [mparams[:weights]'*(rc_adj_val_shock_exp[i] - rc_adj_val_exp_shock_exp[i]) for i = 1:length(dvstar_dk)]

            channels[:robust_insurance] = [-(1000*3/11*rc_adjustment
                .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
                exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]

        end

        # Add up taylor expansion channels
        channels[:total_taylor] = channels[:ce_tax] .+ channels[:unc_adjust] .+
            channels[:smoothing] .+ channels[:precaution] .+
            channels[:output_insurance]

    else

        channels[:robust_insurance] = zeros(size(channels[:smoothing]))

        if opts[:learn_e]
            channels[:scc_direct] = [-(1000*3/11*mparams[:weights2]'*
                (dvstar_dm[i]
                .+ dvstar_dexp[i].*dmu_dm[i]
                .+ dvstar_dsigma[i].*dsigma_de[i]
                .+ dvstar_dl[i].*dl_dm[i])
                .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
                exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]
        else
            channels[:scc_direct] = [-(1000*3/11*mparams[:weights2]'*
                (dvstar_dm[i]
                .+ dvstar_dl[i].*dl_dm[i])
                .*(exog_states[:betaeff][i]/MU(controls[i][1],i)).*
                exog_states[:At][i].*exog_states[:Lt][i])[1] for i = 1:length(dvstar_dk)]
        end

    end

    channels[:scc_direct] = channels[:scc_direct] + channels[:robust_insurance]

    return channels

end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


text:
# Model state variables and choices for model combinations

struct FsmType

    albs::Array{Float64,1}
    Ds::Array{Float64,1}
    Nsnow::Array{Int64,1}
    Sice::Array{Float64,1}
    Sliq::Array{Float64,1}
    theta::Array{Float64,1}
    Tsnow::Array{Float64,1}
    Tsoil::Array{Float64,1}
    Tsurf::Array{Float64,1}

    am::Int64
    cm::Int64
    dm::Int64
    em::Int64
    hm::Int64

    dt::Int64

end

# Model input variables

struct FsmInput

    year::Float64
    month::Float64
    day::Float64
    hour::Float64
    SW::Float64
    LW::Float64
    Sf::Float64
    Rf::Float64
    Ta::Float64
    RH::Float64
    Ua::Float64
    Ps::Float64

end

# Outer constructor for initilizing model state variables

function FsmType(am, cm, dm, em, hm, dt = 3600)

    # Define state variables

    rho_wat = 1000.0                  # Density of water (kg/m^3)
    Tm = 273.15                       # Melting point (K)

    Nsmax = 3                         # Maximum number of snow layers
    Nsoil = 4                         # Number of soil layers

    albs  = Array{Float64}(1)         # Snow albedo
    Ds    = Array{Float64}(Nsmax)     # Snow layer thicknesses (m)
    Nsnow = Array{Int64}(1)           # Number of snow layers
    Sice  = Array{Float64}(Nsmax)     # Ice content of snow layers (kg/m^2)
    Sliq  = Array{Float64}(Nsmax)     # Liquid content of snow layers (kg/m^2)
    theta = Array{Float64}(Nsoil)     # Volumetric moisture content of soil layers
    Tsnow = Array{Float64}(Nsmax)     # Snow layer temperatures (K)
    Tsoil = Array{Float64}(Nsoil)     # Soil layer temperatures (K)
    Tsurf = Array{Float64}(1)         # Surface skin temperature (K)

    # No snow in initial state

    albs[:] = 0.8
    Ds[:] = 0
    Nsnow[:] = 0
    Sice[:] = 0
    Sliq[:] = 0
    Tsnow[:] = Tm

    # Initial soil profiles

    fcly = 0.3
    fsnd = 0.6
    Vsat = 0.505 - 0.037*fcly - 0.142*fsnd

    fsat = 0.5 * ones(Float64, Nsoil)  # Initial moisture content of soil layers as fractions of saturation
    Tsoil[:] = 285.0
    Tsurf[1] = Tsoil[1]
    for k = 1:Nsoil
	theta[k] = fsat[k]*Vsat
    end

    FsmType(albs, Ds, Nsnow, Sice, Sliq, theta, Tsnow, Tsoil, Tsurf, am, cm, dm, em, hm, dt)

end

# Run the model for a time series of meteorological input data

function run_fsm!(hs, md::FsmType, metdata)

    for itime = 1:size(metdata,1)

	year  = metdata[itime, 1]
	month = metdata[itime, 2]
	day   = metdata[itime, 3]
	hour  = metdata[itime, 4]
	SW    = metdata[itime, 5]
	LW    = metdata[itime, 6]
	Sf    = metdata[itime, 7]
	Rf    = metdata[itime, 8]
	Ta    = metdata[itime, 9]
	RH    = metdata[itime, 10]
	Ua    = metdata[itime, 11]
	Ps    = metdata[itime, 12]

	ccall((:fsm_, fsm),
              Void,
              (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},
	       Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},
	       Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Int64},
               Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},
               Ptr{Float64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64}),
	      &year, &month, &day, &hour, &SW, &LW, &Sf, &Rf, &Ta, &RH, &Ua, &Ps,
	      md.albs, md.Ds, md.Nsnow, md.Sice, md.Sliq, md.theta, md.Tsnow, md.Tsoil, md.Tsurf,
	      &md.am, &md.cm, &md.dm, &md.em, &md.hm, &md.dt)

	hs[itime] = sum(md.Ds)

    end

    return nothing

end

function run_fsm(md::FsmType, metdata)

    hs = zeros(Float64, size(metdata,1))

    run_fsm!(hs, md::FsmType, metdata)

    return hs

end

# Run the model for only one time step using FsmInput type

function run_fsm(md::FsmType, id::FsmInput)

    # Call fsm

    ccall((:fsm_, fsm),
          Void,
          (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},
	   Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},
	   Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Int64},
           Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},
           Ptr{Float64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64},Ptr{Int64}),
	  &id.year, &id.month, &id.day, &id.hour, &id.SW, &id.LW, &id.Sf, &id.Rf, &id.Ta, &id.RH, &id.Ua, &id.Ps,
	  md.albs, md.Ds, md.Nsnow, md.Sice, md.Sliq, md.theta, md.Tsnow, md.Tsoil, md.Tsurf,
	  &md.am, &md.cm, &md.dm, &md.em, &md.hm, &md.dt)

    # Save results

    hs = sum(md.Ds)

    return hs

end


text:
const pyjlmoduletype = pynew()

function pyjlmodule_dir(self::Module)
    ks = Symbol[]
    append!(ks, names(self, all = true, imported = true))
    for m in ccall(:jl_module_usings, Any, (Any,), self)::Vector
        append!(ks, names(m))
    end
    pylist(pyjl_attr_jl2py(string(k)) for k in ks)
end

function pyjlmodule_seval(self::Module, expr::Py)
    Py(self.eval(Meta.parse(pyconvert(String, expr))))
end

function init_jlwrap_module()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class ModuleValue(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        def __dir__(self):
            return ValueBase.__dir__(self) + self._jl_callmethod($(pyjl_methodnum(pyjlmodule_dir)))
        def seval(self, expr):
            return self._jl_callmethod($(pyjl_methodnum(pyjlmodule_seval)), expr)
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlmoduletype, jl.ModuleValue)
end

pyjl(v::Module) = pyjl(pyjlmoduletype, v)


text:
MODULE_NAMES = ["LsoBase", "Obj", "Ls", "GdOpt", "Plotting", "Tst", "Eval"]

# import modules
println("Importing modules...")
push!(LOAD_PATH, pwd()*"/mod")
for name in MODULE_NAMES
    try
        eval(Expr(:import, symbol(name))) # metaprogramming: import <name>
        println("Imported $(name)")
    catch e
        println("ERROR: Could not import $(name):\n$e")
    end
end

# ease reloading
reload(name::AbstractString) = Base.reload(name)
function reload()
    for name in MODULE_NAMES
        reload(name)
    end
end

# run tests
#Tst.all()

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LoopTools"
version = v"2.15"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://www.feynarts.de/looptools/LoopTools-2.15.tar.gz", "a065ffdc4fe6882aa3bb926134ba8ec875d6c0a633c3d4aa5f70db26542713f2"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/LoopTools-*
atomic_patch -p1 ../patches/simplify-configure.patch
export AR=ar
export RANLIB=ranlib
./configure --prefix=$prefix --host=${target}
make -j${nproc}
make install EXE="${exeext}"
# Now let's build the shared library
cd "${prefix}/lib"
if [[ "${target}" == *-apple-* ]]; then
    whole_archive="-all_load"
    no_whole_archive="-noall_load"
else
    whole_archive="--whole-archive"
    no_whole_archive="--no-whole-archive"
fi
mkdir -p "${libdir}"
gfortran -fPIC -shared -Wl,${whole_archive} libooptools.a -Wl,${no_whole_archive} -o "${libdir}/libooptools.${dlext}"
rm  libooptools.a
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libooptools", :libooptools),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
module JuliaApp

    export JuliaAppGUI
    export NewtonR

    include("JuliaAppGUI.jl")
    include("Method2.jl")

end # module


text:
"""
    evaluate_proposed(initxy, dx, dy) -> Tuple{Array{Float64,1},String}

Identify location of xy coordinate and compute correct steplengths/displacement vectors.
"""
function evaluate_proposed(initxy::Array{Float64,1}, dx::Float64, dy::Float64)
    if inwater(initxy)
        return watercalc(initxy, dx, dy)
    elseif inenz(initxy)
        return enzcalc(initxy, dx, dy)
    elseif inppd(initxy)
        return ppdcalc(initxy, dx, dy)
    else
        println("evaluate_proposed Error", initxy)
        return [-1.0, -1.0], "evaluate_proposed Error"
    end
end

"Helper when xy is initially in water"
function watercalc(initxy::Array{Float64,1}, dx::Float64, dy::Float64)::Tuple{Array{Float64,1},String}
    initx, inity = initxy
    flowBias = flow_arlett(initxy)
    proposedxy = [initx + flowBias + WATER_STEP_SIZE * dx, inity + WATER_STEP_SIZE * dy]

    if inwater(proposedxy)
        return proposedxy, "water"
    elseif inenz(proposedxy)
        directionOfEntry = whereoutsidespawn(initxy) # tail in water
        if directionOfEntry == "whereoutsidespawn Error"
            println(directionOfEntry, " in watercalc ", initxy, proposedxy)
        elseif directionOfEntry == "N"
            confirmxy = north!(proposedxy, initxy, dx, dy, "water", "enz")
        elseif directionOfEntry == "NE"
            confirmxy = WATER_TO_ENZ_NORTHEAST
        elseif directionOfEntry == "NW"
            confirmxy = WATER_TO_ENZ_NORTHWEST
        else # E or W
            confirmxy = eastwest!(proposedxy, initxy, dx, dy, "water", "enz", directionOfEntry)
        end
        return confirmxy, "enz"
    elseif inppd(proposedxy)
        confirmxy = north!(proposedxy, initxy, dx, dy, "water", "ppd")
        return confirmxy, "ppd"
    else # collision
        return proposedxy, "water"
    end
end

"Helper when xy is initially in the enzyme layer"
function enzcalc(initxy::Array{Float64,1}, dx::Float64, dy::Float64)::Tuple{Array{Float64,1},String}
    initx, inity = initxy
    flowBias = flow_arlett(initxy) * ENZ_STEP_SIZE / WATER_STEP_SIZE
    proposedxy = [initx + flowBias + ENZ_STEP_SIZE * dx, inity + ENZ_STEP_SIZE * dy]

    if inenz(proposedxy)
        return proposedxy, "enz"
    elseif inwater(proposedxy)
        directionOfEntry = whereoutsidespawn(proposedxy) # head in water
        if directionOfEntry == "whereoutsidespawn Error"
            println(directionOfEntry, " in enzcalc ", initxy, proposedxy)
        elseif directionOfEntry == "N"
            confirmxy = north!(proposedxy, initxy, dx, dy, "enz", "water")
        elseif directionOfEntry == "NE"
            confirmxy = ENZ_TO_WATER_NORTHEAST
        elseif directionOfEntry == "NW"
            confirmxy = ENZ_TO_WATER_NORTHWEST
        else # E or W
            confirmxy = eastwest!(proposedxy, initxy, dx, dy, "enz", "water", directionOfEntry)
        end
        return confirmxy, "water"
    elseif inppd(proposedxy)
        confirmxy = north!(proposedxy, initxy, dx, dy, "enz", "ppd")
        return confirmxy, "ppd"
    else # collision
        return proposedxy, "enz"
    end
end

"Helper when xy is initially in the ppd layer"
function ppdcalc(initxy::Array{Float64,1}, dx::Float64, dy::Float64)::Tuple{Array{Float64,1},String}
    initx, inity = initxy
    PPD_STEP_SIZE = STEP_SIZE_DICT["ppd"]
    proposedxy = [initx + PPD_STEP_SIZE * dx, inity + PPD_STEP_SIZE * dy] # no flow

    if inppd(proposedxy)
        return proposedxy, "ppd"
    elseif inenz(proposedxy)
        confirmxy = north!(proposedxy, initxy, dx, dy, "ppd", "enz")
        return confirmxy, "enz"
    elseif inwater(proposedxy)
        confirmxy = north!(proposedxy, initxy, dx, dy, "ppd", "water")
        return confirmxy, "water"
    else # collision
        return proposedxy, "ppd"
    end
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
function print_house(house, floor)
    for i in 1:3
        for j in 1:10
            print(" ", house[floor, i, j])
        end
        println()
    end
end

house = zeros(Int, 4, 3, 10)
n = parse(Int, readline())
for i in 1:n
    b, f, r, v = parse.(Int, split(readline()))
    house[b,f,r] += v
end

print_house(house, 1)
println("#"^20)
print_house(house, 2)
println("#"^20)
print_house(house, 3)
println("#"^20)
print_house(house, 4)


text:
export SoilWaterModel, PrescribedWaterModel, get_water_content

abstract type AbstractWaterModel <: AbstractSoilComponentModel end

"""
    PrescribedWaterModel{F1, F2} <: AbstractWaterModel

Model structure for a prescribed water content model.

The user supplies functions of space and time for both `ϑ_l` and
`θ_i`. No auxiliary or state variables are added, no PDE is solved.
The defaults are no moisture anywhere, for all time.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct PrescribedWaterModel{FN1, FN2} <: AbstractWaterModel
    "Augmented liquid fraction"
    ϑ_l::FN1
    "Volumetric fraction of ice"
    θ_i::FN2
end

"""
    function PrescribedWaterModel(
        ϑ_l::Function = (aux, t) -> eltype(aux)(0.0),
        θ_i::Function = (aux, t) -> eltype(aux)(0.0),
    )
        args = (ϑ_l, θ_i)
        return PrescribedWaterModel{typeof.(args)...}(args...)
    end
Outer constructor for the PrescribedWaterModel defining default values.

The functions supplied by the user are point-wise evaluated and are
evaluated in the Balance Law functions compute_gradient_argument,
nodal_update, etc. whenever the prescribed water content variables are
needed by the heat model.
"""
function PrescribedWaterModel(
    ϑ_l::Function = (aux, t) -> eltype(aux)(0.0),
    θ_i::Function = (aux, t) -> eltype(aux)(0.0),
)
    args = (ϑ_l, θ_i)
    return PrescribedWaterModel{typeof.(args)...}(args...)
end

"""
    SoilWaterModel{FT, IF, VF, MF, HM, Fiϑl, Fiθi, BCT} <: AbstractWaterModel

The necessary components for solving the equations for water (liquid or ice) in soil.

Without freeze/thaw source terms added (separately), this model reduces to
Richard's equation for liquid water. Note that the default for `θ_i` is zero.
Without freeze/thaw source terms added to both the liquid and ice equations,
the default should never be changed, because we do not enforce that the total
volumetric water fraction is less than or equal to porosity otherwise.

When freeze/thaw source terms are included, this model encompasses water in both
liquid and ice form, and water content is conserved upon phase change.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SoilWaterModel{FT, IF, VF, MF, HM, Fiϑl, Fiθi, BCT} <: AbstractWaterModel
    "Impedance Factor - will be 1 or ice dependent"
    impedance_factor::IF
    "Viscosity Factor - will be 1 or temperature dependent"
    viscosity_factor::VF
    "Moisture Factor - will be 1 or moisture dependent"
    moisture_factor::MF
    "Hydraulics Model - used in matric potential and moisture factor of hydraulic conductivity"
    hydraulics::HM
    "Initial condition: augmented liquid fraction"
    initialϑ_l::Fiϑl
    "Initial condition: volumetric ice fraction"
    initialθ_i::Fiθi
    "Boundary Condition Type"
    boundaries::BCT
end

"""
    SoilWaterModel(
        ::Type{FT};
        impedance_factor::AbstractImpedanceFactor{FT} = NoImpedance{FT}(),
        viscosity_factor::AbstractViscosityFactor{FT} = ConstantViscosity{FT}(),
        moisture_factor::AbstractMoistureFactor{FT} = MoistureIndependent{FT}(),
        hydraulics::AbstractHydraulicsModel{FT} = vanGenuchten{FT}(),
        initialϑ_l = (aux) -> FT(NaN),
        initialθ_i = (aux) -> FT(0.0),
        boundaries::AbstractBoundaryConditions,
    ) where {FT}

Constructor for the SoilWaterModel. Defaults imply a constant K = K_sat model.
"""
function SoilWaterModel(
    ::Type{FT};
    impedance_factor::AbstractImpedanceFactor{FT} = NoImpedance{FT}(),
    viscosity_factor::AbstractViscosityFactor{FT} = ConstantViscosity{FT}(),
    moisture_factor::AbstractMoistureFactor{FT} = MoistureIndependent{FT}(),
    hydraulics::AbstractHydraulicsModel{FT} = vanGenuchten{FT}(),
    initialϑ_l::Function = (aux) -> eltype(aux)(NaN),
    initialθ_i::Function = (aux) -> eltype(aux)(0.0),
    boundaries::AbstractBoundaryConditions,
) where {FT}
    args = (
        impedance_factor,
        viscosity_factor,
        moisture_factor,
        hydraulics,
        initialϑ_l,
        initialθ_i,
        boundaries,
    )
    return SoilWaterModel{FT, typeof.(args)...}(args...)
end

"""
    get_water_content(
        water::SoilWaterModel,
        aux::Vars,
        state::Vars,
        t::Real
    )

Return the moisture variables for the balance law soil water model.
"""
function get_water_content(
    water::SoilWaterModel,
    aux::Vars,
    state::Vars,
    t::Real,
)
    FT = eltype(state)
    return FT(state.soil.water.ϑ_l), FT(state.soil.water.θ_i)
end

"""
    get_water_content(
        water::PrescribedWaterModel,
        aux::Vars,
        state::Vars,
        t::Real
    )

Return the moisture variables for the prescribed soil water model.
"""
function get_water_content(
    water::PrescribedWaterModel,
    aux::Vars,
    state::Vars,
    t::Real,
)
    FT = eltype(aux)
    ϑ_l = water.ϑ_l(aux, t)
    θ_i = water.θ_i(aux, t)
    return FT(ϑ_l), FT(θ_i)
end

"""
    function get_diffusive_water_term(
        water::SoilWaterModel,
        diffusive::Vars
    )

Returns the diffusive water term from the `diffusive` vector.
"""
function get_diffusive_water_term(water::SoilWaterModel, diffusive::Vars)
    return diffusive.soil.water.K∇h
end

"""
    function get_diffusive_water_term(
        water::PrescribedWaterModel,
        diffusive::Vars
    )

Returns zero for the diffusive water term, under the PrescribedWaterModel.
"""
function get_diffusive_water_term(water::PrescribedWaterModel, diffusive::Vars)
    FT = eltype(diffusive)
    return SVector{3, FT}(0, 0, 0)
end

vars_state(water::SoilWaterModel, st::Prognostic, FT) = @vars(ϑ_l::FT, θ_i::FT)

vars_state(water::SoilWaterModel, st::Auxiliary, FT) = @vars(h::FT, K::FT)

vars_state(water::SoilWaterModel, st::Gradient, FT) = @vars(h::FT)

vars_state(::SoilWaterModel, st::GradientFlux, FT) = @vars(K∇h::SVector{3, FT})#really, the flux is - K∇h

function soil_init_aux!(
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    aux::Vars,
    geom::LocalGeometry,
)
    FT = eltype(aux)
    T = get_initial_temperature(land.soil.heat, aux, FT(0.0))
    S_l = effective_saturation(
        soil.param_functions.porosity,
        water.initialϑ_l(aux),
    )
    ψ = pressure_head(
        water.hydraulics,
        soil.param_functions.porosity,
        soil.param_functions.S_s,
        water.initialϑ_l(aux),
        water.initialθ_i(aux),
    )
    aux.soil.water.h = hydraulic_head(aux.z, ψ)
    aux.soil.water.K =
        soil.param_functions.Ksat * hydraulic_conductivity(
            water.impedance_factor,
            water.viscosity_factor,
            water.moisture_factor,
            water.hydraulics,
            water.initialθ_i(aux),
            soil.param_functions.porosity,
            T,
            S_l,
        )
end

function land_nodal_update_auxiliary_state!(
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    T = get_temperature(land.soil.heat, aux, t)
    S_l = effective_saturation(
        soil.param_functions.porosity,
        state.soil.water.ϑ_l,
    )
    ψ = pressure_head(
        water.hydraulics,
        soil.param_functions.porosity,
        soil.param_functions.S_s,
        state.soil.water.ϑ_l,
        state.soil.water.θ_i,
    )
    aux.soil.water.h = hydraulic_head(aux.z, ψ)
    aux.soil.water.K =
        soil.param_functions.Ksat * hydraulic_conductivity(
            water.impedance_factor,
            water.viscosity_factor,
            water.moisture_factor,
            water.hydraulics,
            state.soil.water.θ_i,
            soil.param_functions.porosity,
            T,
            S_l,
        )
end

function compute_gradient_argument!(
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)

    S_l = effective_saturation(
        soil.param_functions.porosity,
        state.soil.water.ϑ_l,
    )
    ψ = pressure_head(
        water.hydraulics,
        soil.param_functions.porosity,
        soil.param_functions.S_s,
        state.soil.water.ϑ_l,
        state.soil.water.θ_i,
    )
    transform.soil.water.h = hydraulic_head(aux.z, ψ)

end

function compute_gradient_flux!(
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.soil.water.K∇h = aux.soil.water.K * ∇transform.soil.water.h
end

function flux_second_order!(
    land::LandModel,
    soil::SoilModel,
    water::SoilWaterModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux.soil.water.ϑ_l -= diffusive.soil.water.K∇h
end

max_forks_repo_name: jm-c/CLIMA


text:
module JSServe

import Sockets
using UUIDs, Hyperscript, JSON3, Observables
import Sockets: send, TCPServer
using Hyperscript: Node, children, tag
using HTTP, Markdown
using HTTP: Response, Request
using HTTP.Streams: Stream
using WebSockets
using WebSockets: WebSocket

using Base64
using MsgPack
using WidgetsBase
using WidgetsBase: vertical, horizontal
using SHA
using Tables
using Colors
using LinearAlgebra
using CodecZlib

include("types.jl")
include("server.jl")
include("js_source.jl")
include("session.jl")
include("observables.jl")
include("dependencies.jl")
include("http.jl")
include("util.jl")
include("widgets.jl")
include("hyperscript_integration.jl")
include("display.jl")
include("markdown_integration.jl")
include("serialization.jl")
include("offline.jl")
include("browser_display.jl")

const JSSERVE_CONFIGURATION = (
    # The URL used to which the default server listens to
    listen_url = Ref("127.0.0.1"),
    # The Port to which the default server listens to
    listen_port = Ref(9284),
    # The url Javascript uses to connect to the websocket.
    # if empty, it will use:
    # `window.location.protocol + "//" + window.location.host`
    external_url = Ref(""),
    # The url prepended to assets when served!
    # if `""`, urls are inserted into HTML in relative form!
    content_delivery_url = Ref(""),
    # Verbosity for logging!
    verbose = Ref(false)
)

"""
    configure_server!(;
            listen_url::String=JSSERVE_CONFIGURATION.listen_url[],
            listen_port::Integer=JSSERVE_CONFIGURATION.listen_port[],
            forwarded_port::Integer=listen_port,
            external_url=nothing,
            content_delivery_url=nothing
        )

Configures the parameters for the automatically started server.

    Parameters:

    * listen_url=JSSERVE_CONFIGURATION.listen_url[]
        The address the server listens to.
        must be 0.0.0.0, 127.0.0.1 or localhost.
        If not set differently by an ENV variable, will default to 127.0.0.1

    * listen_port::Integer=JSSERVE_CONFIGURATION.listen_port[],
        The Port to which the default server listens to
        If not set differently by an ENV variable, will default to 9284

    * forwarded_port::Integer=listen_port,
        if port gets forwarded to some other port, set it here!

    * external_url=nothing
        The url from which the server is reachable.
        If served on "127.0.0.1", this will default to http://localhost:forwarded_port
        if listen_url is "0.0.0.0", this will default to http://\$(Sockets.getipaddr()):forwarded_port
        so that the server is reachable inside the local network.
        If the server should be reachable from some external dns server,
        this needs to be set here.

    * content_delivery_url=nothing
        You can server files from another server.
        Make sure any file referenced from Julia is reachable at
        content_delivery_url * "/the_file"
"""
function configure_server!(;
        listen_url=nothing,
        # The Port to which the default server listens to
        listen_port::Integer=JSSERVE_CONFIGURATION.listen_port[],
        # if port gets forwarded to some other port, set it here!
        forwarded_port::Integer=listen_port,
        # The url from which the server is reachable.
        # If served on "127.0.0.1", this will default to
        # if listen_url is "0.0.0.0", this will default to
        # Sockets.getipaddr() so that it's the ip address in the local network.
        # If the server should be reachable from some external dns server,
        # this needs to be set here
        external_url=nothing,
        # You can server files from another server.
        # Make sure any file referenced from Julia is reachable at
        # content_delivery_url * "/the_file"
        content_delivery_url=nothing
    )

    if isnothing(listen_url)
        if !isnothing(external_url)
            # if we serve to an external url, server must listen to 0.0.0.0
            listen_url = "0.0.0.0"
        else
            listen_url = JSSERVE_CONFIGURATION.listen_url[]
        end
    end

    if isnothing(external_url)
        if listen_url == "0.0.0.0"
            external_url = string(Sockets.getipaddr(), ":$forwarded_port")
        elseif listen_url in ("127.0.0.1", "localhost")
            external_url = "http://localhost:$forwarded_port"
        else
            error("Trying to listen to $(listen_url), while only \"127.0.0.1\", \"0.0.0.0\" and \"localhost\" are supported")
        end
    end
    # set the config!
    JSSERVE_CONFIGURATION.listen_url[] = listen_url
    JSSERVE_CONFIGURATION.external_url[] = external_url
    JSSERVE_CONFIGURATION.listen_port[] = listen_port
    if content_delivery_url === nothing
        JSSERVE_CONFIGURATION.content_delivery_url[] = external_url
    else
        JSSERVE_CONFIGURATION.content_delivery_url[] = content_delivery_url
    end
    return
end

function has_html_display()
    for display in Base.Multimedia.displays
        # Ugh, why would textdisplay say it supports HTML??
        display isa TextDisplay && continue
        displayable(display, MIME"text/html"()) && return true
    end
    return false
end

function __init__()
    url = if haskey(ENV, "JULIA_WEBIO_BASEURL")
        ENV["JULIA_WEBIO_BASEURL"]
    else
        ""
    end
    if endswith(url, "/")
        url = url[1:end-1]
    end
    JSSERVE_CONFIGURATION.listen_url[] = get(ENV, "JSSERVE_LISTEN_URL", "127.0.0.1")
    JSSERVE_CONFIGURATION.external_url[] = url
    JSSERVE_CONFIGURATION.content_delivery_url[] = url

    if haskey(ENV, "WEBIO_HTTP_PORT")
        JSSERVE_CONFIGURATION.listen_port[] = parse(Int, ENV["WEBIO_HTTP_PORT"])
    end

    # If there is no html inline display in the IDE that JSServe is running
    # we display things in the local browser
    if !has_html_display()
        browser_display()
    end
end

# Core functionality
export Page, Session, App, DOM, @js_str
export Slider, Button

end # module


text:
using Poly
export @poly_loop, @depends_on, sub_sym_in_ex
using MacroTools

"""
helpers for substituting values for symbols in an expression
    will sub an interpolated symbol for the value as well as a non-interpolated symbol
"""
function sub_sym_in_ex(ex::Expr, sub_sym::Symbol, val::Number)
    if ex.head == :$ && ex.args[1] == sub_sym
        return val
    end
    new_args = []
    for arg in ex.args
        push!(new_args, sub_sym_in_ex(arg, sub_sym, val))
    end
    return Expr(ex.head, new_args...)
end

function sub_sym_in_ex(sym::Symbol, sub_sym::Symbol, val::Number)
    if sym == sub_sym
        return val
    end
    return sym
end

function sub_sym_in_ex(any, sub_sym::Symbol, val::Number)
    return any
end

"""
get instructions out of block with cond
"""
function get_insts_in_block(args, cond::Expr, parent_doms::Vector{Domain}, dom::Domain)::Vector{Instruction}
    instructions = []
    for arg in args
        if typeof(arg) != LineNumberNode
            symbol = arg
            while typeof(symbol) != Symbol
                symbol = symbol.args[1]
            end
            iname = gensym(string(symbol))
            inst = Instruction(iname, arg, Set(), cond)
            for d in parent_doms
                push!(inst.dependencies, d.iname)
            end
            push!(inst.dependencies, dom.iname)
            push!(dom.instructions, inst)
            push!(instructions, inst)
        end
    end
    return instructions
end

"""
helper for macro
"""
function poly_loop(ex::Expr; parent_doms=Domain[])::LoopKernel
    # generate loop domain
    bounds = ex.args[1]
    iname = bounds.args[1]
    space = bounds.args[2]
    if typeof(space) != Expr || !(length(space.args) == 3 || length(space.args) == 4)
        error("expected a loop in the form bound:bound or bound:step:bound")
    end
    if length(space.args) == 3
        # bound : bound
        lowerbound = space.args[2]
        upperbound = space.args[3]
        step = 1
    else
        # bound : step : bound
        lowerbound = space.args[2]
        upperbound = space.args[4]
        step = space.args[3]
    end
    dom = Domain(iname, lowerbound, upperbound, step, Set(), [])
    domains = [dom]
    new_parent_doms = copy(parent_doms)
    append!(new_parent_doms, [dom])

    # generate loop body instructions
    body = ex.args[2]
    instructions = []
    for line in body.args
        if typeof(line) == Expr && line.head == :for
            # nested loops
            subkern = poly_loop(line, parent_doms=new_parent_doms)
            append!(instructions, subkern.instructions)
            append!(domains, subkern.domains)
            push!(dom.instructions, subkern.domains[1])
            for d in parent_doms
                push!(subkern.domains[1].dependencies, d.iname)
            end
            push!(subkern.domains[1].dependencies, dom.iname)
        elseif typeof(line) == Expr && line.head == :if
            # if block
            ifcond = line.args[1]
            # instructions in if block
            append!(instructions, get_insts_in_block(line.args[2].args, ifcond, parent_doms, dom))
            # elseif or else or both
            if length(line.args) == 3
                rest = line.args[3]
                if rest.head == :elseif
                    # elseif block
                    elifcond = rest.args[1].args[2] # since line number node
                    elifcond = :($elifcond && !($ifcond))
                    # instructions in elseif block
                    append!(instructions, get_insts_in_block(rest.args[2].args, elifcond, parent_doms, dom))
                    # else block also
                    if length(rest.args) == 3
                        # else block
                        elsecond = :(!($ifcond) && !($elifcond))
                        # instructions in else block
                        append!(instructions, get_insts_in_block(rest.args[3].args, elsecond, parent_doms, dom))
                    end
                else
                    # else block
                    elsecond = :(!($ifcond))
                    # instructions in else block
                    append!(instructions, get_insts_in_block(rest.args, elsecond, parent_doms, dom))
                end
            end
        elseif typeof(line) != LineNumberNode
            symbol = line
            while typeof(symbol) != Symbol
                symbol = symbol.args[1]
            end
            iname = gensym(string(symbol))
            inst = Instruction(iname, line, Set(), :())
            push!(instructions, inst)
            for d in parent_doms
                push!(inst.dependencies, d.iname)
            end
            push!(inst.dependencies, dom.iname)
            push!(dom.instructions, inst)
        end
    end

    kern = LoopKernel(instructions, domains, [], [], [])
    return kern
end

"""
helper to get the value in a value
"""
function __get(::Type{Val{N}}) where N
    return N
end

"""
Runs a Polyhedral model on the input for loop. May reorder instructions (including loop orderings) and vectorize code. Only the outermost loop should use the macro. By default, will also tile code. See below for more details.

# Example
```jldoctest
@poly_loop for i = 1:n
    for j = 1:m
        for k = 1:r
            out[i, j] += A[i, k] * B[k, j]
        end
        out[i, j] *= 2
    end
end
```
(runs a matrix multiplication and double)

!!! note "Loop Types"
    `@poly_loop` requires a normal for loop (i.e `i = lowerbound:upperbound` or `i = lowerbound:step:upperbound`)

## Tiling
Loops will be automatically tiled based on l1 cache sizes unless `tile=0` is set.
### Example No Tiling:
```jldoctest
@poly_loop tile=0 for i=1:n
    arr[i] = 1.0
end
```
Tiling is usually faster, but in cases where the loop bounds are small it will just introduce more loop overhead. Unless `tile=0` is set, the outermost band (grouping of permutable loops) will be tiled.

!!! warning "Invalid Tiling"
    Some non-uniform access patterns can result in invalid tiling. This is pretty rare. If this occurs, code will likely error during compilation.

## Threading
To enable multithreading of the outermost loop, pass thread=true (disabled by default):
### Example:
```jldoctest
@poly_loop thread=true for i=1:n
    arr[i] = 1.0
end
```

!!! note "Threading Condition"
    Julia must be run with `julia --threads=auto` (or a number of threads instead of `auto`) to enable threading at all.

## Debugging
Debugging and verbosity options (see run_polyhedral_model for details) can be passed like:
```jldoctest
julia> @poly_loop debug=true verbose=2 for i=1:4:n-1
    arr[i] *= arr[i + 1]
end
```
This will print out the following information, depending on the verbosity level set:
- 0: no printing (unless `debug=true`, where ISL errors are printed)
- 1: final Julia expression is printed
- 2: initial C code, loop orderings, new schedule, final C code, and final Julia expression printed
- 3: all of 2 plus domain, access relations, dependence relations, original schedule, and new schedule constraints printed

## Functions

!!! warning "Functions"
    Function calls are not currently supported in most cases

If a function modifies any inputs, then there is no way to know which inputs are modified or how those inputs are modified. If a function returns an output such as a matrix, the write to each individual index of the matrix (for example `A[i, j]`) cannot be inferrred by ISL. In addition, `@inbounds` is not neccessary (or allowed) as the aggresive transformation will add `@inbounds` automatically.

### Loop Bounds
This also applies to loop bounds. One easy workaround is something like:
```jldoctest
juila> n = size(out, 1)
julia> @poly_loop for i = 1:n
    out[i] += 1.0
end
```

### Exceptions
The only exceptions are the following functions (which can be used in loop bounds and in instructions):
`min()`:
```jldoctest
@poly_loop for i=1:min(n, m)
    arr[i, i] = 1.0
end
```
`max()`:
```jldoctest
@poly_loop for i=1:max(n, m)
    arr[i, i] -= 1.0
end
```
`floor()` -> DO NOT cast to `Int()` (will be added later):
```jldoctest
@poly_loop for i=1:floor(n/2)
    arr[i] += 1.0
end
```

!!! note "Loop Names"
    All loop iterators must have unique symbols, even if in different scopes. This is so an original schedule can be extracted from the code.

## Example
```jldoctest
@poly_loop for i=1:n
    for j=1:floor(n/2)
        arr[i, j] += 1.0
    end
    for jj=floor(n/2):n
        arr[i, j] *= 2.0
    end
end
```

## Interpolation

!!! warning "Striding"
    ISL does not work with non-numerical striding, i.e. `for i=1:TILE:n`.

As a workaround, this macro supports interpolating the values into the function so that ISL has access to the numerical values (and non-numerical strides MUST be interpolated). This delays the analysis to runtime, so it is advised that interpolated values are typically constant so that compilation does not occur over and over.

### Examples:
```jldoctest
@poly_loop for i=1:\$c:n
    arr[i] = 1
end
```
```jldoctest
@poly_loop for i=1:\$c:\$n
    arr[i] = 1
end
```
In addition, any constants can be interpolated, and sometimes this may improve the performance of ISL. Numerical bounds can help with tiling and inference on spaces. However, it is important that interpolated values are actually constants, or very rarely change, so that this delayed compilation does not occur regularly during runtime.

!!! note "Indexing"
    All array indices must be in terms of loop bounds or constants, not of variables.

### Example:
```jldoctest
julia> c = i + j
julia> A[c]
```
Must become
```jldoctest
julia> A[i + j]
```

!!! warning "Multiplicative Indexing"
    Multilpicative indexing (i.e. `A[i*t, j]`) is not supported. If needed, stride over the iterator.

#### Example:
```jldoctest
NUM_TILES = Int(N/TILE_SIZE)
@poly_loop for t=1:NUM_TILES
    A[TILE_SIZE*t] = 1
end
```
Would need to become this:
```jldoctest
@poly_loop for t=1:\$TILE_SIZE:N
    A[t] = 1
end
```

## If/elseif/else blocks:

!!! note "If Block Conditions"
    If blocks must contain conditions on the loop iterators and constants only (not on any elements of an array, for example)

### Example:
```jldoctest
@poly_loop for i=1:n
    if i%2 == 0
        A[i] = 1
    else
        A[i] = 2
    end
end
```

"""
macro poly_loop(ex0...)
    if length(ex0) > 4
        error("expected at most 3 options and a loop, got: ", length(ex0), "elements total")
    end
    debug = false
    verbose = 0
    tile = -1 # all tiles by default
    thread = false # don't thread by default
    # get set options
    for arg in ex0[1:end-1] # last is for loop
        if arg.args[1] == :debug
            debug = arg.args[2]
        elseif arg.args[1] == :verbose
            verbose = arg.args[2]
        elseif arg.args[1] == :tile
            tile = arg.args[2]
        elseif arg.args[1] == :thread
            thread = arg.args[2]
        end
    end
    if typeof(debug) != Bool
        error("expected boolean for debug, got: ", debug)
    end
    if typeof(verbose) != Int
        error("expected integer for verbose, got: ", verbose)
    end
    if typeof(tile) != Int
        error("expected integer for tile, got: ", tile)
    elseif tile < -1
        error("expected tile >= -1, got: ", tile)
    end
    if typeof(thread) != Bool
        error("expected boolean for thread, got: ", thread)
    end

    ex0 = ex0[end] # expression is last arg

    if ex0.head != :for
        error("expected a for loop, got: ", ex0.head)
    end

    liftedex = copy(ex0)
    letargs = Base._lift_one_interp!(liftedex)
    if letargs != []
        # have interpolation, need to evaluate at runtime to get values of symbols
        sub_symbols = []
        # gather all the symbols that need to be replaced
        for arg in letargs
            push!(sub_symbols, arg.args[2].args[1])
        end
        # compile the Val() expressions that will be passed to the generated function
        exprs = []
        for sym in sub_symbols
            push!(exprs, Expr(:call, :Val, Expr(:escape, sym)))
        end
        # get the function args for kernel
        func_args = collect(get_kernel_args(poly_loop(ex0), use_defined=true))
        # unique name
        name = gensym("delay_poly_loop")
        eval(quote
            # generated function will have the values of the symbols for replacing at runtime
            @generated function $name($(func_args...), vals...)
                expr = $(QuoteNode(ex0))
                syms = $(QuoteNode(sub_symbols))
                d = $(QuoteNode(debug))
                v = $(QuoteNode(verbose))
                t = $(QuoteNode(tile))
                th = $(QuoteNode(thread))
                for (i, val) in enumerate(vals)
                    N = __get(val)
                    expr = sub_sym_in_ex(expr, syms[i], N)
                end
                kernel = poly_loop(expr)
                expr = compile_expr(kernel, debug=d, verbose=v, tile=t, thread=th)
                expr
            end
        end)
        quote
            $name($(map(esc, func_args)...), $(exprs...))
        end
    else
        # make kernel
        kernel = poly_loop(ex0)
        expr = MacroTools.prettify(compile_expr(kernel, debug=debug, verbose=verbose, tile=tile, thread=thread))
        esc(quote
            $(expr)
        end)
    end
end


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The multi-commodity flow problem

# JuMP implementation of the multicommodity transportation model AMPL: A Modeling
# Language for Mathematical Programming, 2nd ed by Robert Fourer, David Gay, and
# Brian W. Kernighan 4-1.
#
# Originally contributed by Louis Luangkesorn, February 26, 2015.

using JuMP
import GLPK
import Test

function example_multi(; verbose = true)
    orig = ["GARY", "CLEV", "PITT"]
    dest = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"]
    prod = ["bands", "coils", "plate"]
    numorig = length(orig)
    numdest = length(dest)
    numprod = length(prod)
    ## supply(prod, orig) amounts available at origins
    supply = [
        400    700    800;
        800   1600   1800;
        200    300    300
    ]
    ## demand(prod, dest) amounts required at destinations
    demand = [
        300   300   100    75   650   225   250;
        500   750   400   250   950   850   500;
        100   100     0    50   200   100   250
    ]
    ## limit(orig, dest) of total units from any origin to destination
    limit = [625.0 for j in 1:numorig, i in 1:numdest]
    ## cost(dest, orig, prod) Shipment cost per unit
    cost = reshape([
        [
            [  30,   10,    8,   10,   11,   71,    6];
            [  22,    7,   10,    7,   21,   82,   13];
            [  19,   11,   12,   10,   25,   83,   15]
        ];
        [
            [  39,   14,   11,   14,   16,   82,    8];
            [  27,    9,   12,    9,   26,   95,   17];
            [  24,   14,   17,   13,   28,   99,   20]
        ];
        [
            [  41,   15,   12,   16,   17,   86,    8];
            [  29,    9,   13,    9,   28,   99,   18];
            [  26,   14,   17,   13,   31,  104,   20]
        ]
    ], 7, 3, 3)
    ## DECLARE MODEL
    multi = Model(GLPK.Optimizer)
    ## VARIABLES
    @variable(multi, trans[1:numorig, 1:numdest, 1:numprod] >= 0)
    ## OBJECTIVE
    @objective(
        multi,
        Max,
        sum(
            cost[j, i, p] * trans[i, j, p]
            for i in 1:numorig, j in 1:numdest, p in 1:numprod
        )
    )
    ## CONSTRAINTS
    ## Supply constraint
    @constraint(
        multi,
        supply_con[i in 1:numorig, p in 1:numprod],
        sum(trans[i, j, p] for j in 1:numdest) == supply[p, i]
    )
    ## Demand constraint
    @constraint(
        multi,
        demand_con[j in 1:numdest, p in 1:numprod],
        sum(trans[i, j, p] for i in 1:numorig) == demand[p, j]
    )
    ## Total shipment constraint
    @constraint(
        multi,
        total_con[i in 1:numorig, j in 1:numdest],
        sum(trans[i, j, p] for p in 1:numprod) - limit[i, j] <= 0
    )
    optimize!(multi)
    Test.@test termination_status(multi) == MOI.OPTIMAL
    Test.@test primal_status(multi) == MOI.FEASIBLE_POINT
    Test.@test objective_value(multi) == 225_700.0
    if verbose
        println("RESULTS:")
        for i in 1:length(orig)
            for j in 1:length(dest)
                for p in 1:length(prod)
                    print(" $(prod[p]) $(orig[i]) $(dest[j]) = $(value(trans[i, j, p]))\t")
                end
                println()
            end
        end
    end
    return
end

example_multi()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode