text:
# include("pascaltriangle.jl")
# include("in2b.jl")
# include("b2in.jl")
include("in2bEne.jl")

function cutMsizeEne(Msize0::Int64, Np::Int64, matp::Matrix{Int64}, Ene0minumhalf::Int64)

    Msize = Msize0*2

    # matp = zeros(Int64,Msize+1,Np+1)
    # pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    maxmatp = matp[Msize+1,Np+1] # the indices are m+1 and n+1 for N^m_ns

    indvec = zeros(Int64,maxmatp)
    mm = 0

    for nn = 1:maxmatp

        Enenn = in2bEne(nn,Msize,Np,matp)

        if Enenn <= Ene0minumhalf
           mm = mm + 1
           indvec[mm] = nn
        end

    end

    indvec1 = copy(indvec[1:mm])

    return indvec1

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
using Ipopt
using Base.Test

# First of all, test that hs071 example works
include("hs071_test.jl")

# Test setting some options
# String option
println("\nString option")
addOption(prob, "hessian_approximation", "exact")
@test_throws ErrorException addOption(prob, "hessian_approximation", "badoption")
println("\nInt option")
# Int option
addOption(prob, "file_print_level", 3) == nothing
@test_throws ErrorException addOption(prob, "file_print_level", -1)
# Double option
println("\nFloat option")
addOption(prob, "derivative_test_tol", 0.5)
@test_throws ErrorException addOption(prob, "derivative_test_tol", -1.0)

# Test opening an output file
openOutputFile(prob, "blah.txt", 5)

# Test MathProgBase stuff
#include(joinpath(Pkg.dir("MathProgBase"),"test","linprog.jl"))
#linprogtest(IpoptSolver())
    using MathProgBase
    solver = IpoptSolver()
    sol = linprog([-1,0],[2 1],'<',1.5,solver)
    @test sol.status == :Optimal
    @test abs(sol.objval - -0.75) <= 1e-6
    @test norm(sol.sol - [0.75,0.0]) <= 1e-6

    sol = linprog([-1,0],sparse([2 1]),'<',1.5,solver)
    @test sol.status == :Optimal
    @test abs(sol.objval - -0.75) <= 1e-6
    @test norm(sol.sol - [0.75,0.0]) <= 1e-6

    # test infeasible problem:
    # min x
    # s.t. 2x+y <= -1
    # x,y >= 0
    sol = linprog([1,0],[2 1],'<',-1,solver)
#    @test sol.status == :Infeasible

    # test unbounded problem:
    # min -x-y
    # s.t. -x+2y <= 0
    # x,y >= 0
    sol = linprog([-1,-1],[-1 2],'<',[0],solver)
#    @test sol.status == :Unbounded

include(joinpath(Pkg.dir("MathProgBase"),"test","nlp.jl"))
nlptest(IpoptSolver())
nlptest_nohessian(IpoptSolver())
convexnlptest(IpoptSolver())
rosenbrocktest(IpoptSolver())

include(joinpath(Pkg.dir("MathProgBase"),"test","quadprog.jl"))
quadprogtest(IpoptSolver())
qpdualtest(IpoptSolver())

# Test retoration only options
#
# Warm start with infeasible solution, force restoration on initial iteration.
# But limit to 0 iterations. Forces :UserLimit exit.
m = MathProgBase.NonlinearModel(IpoptSolver(start_with_resto="yes", resto_max_iter=0))
l = [1,1,1,1]
u = [5,5,5,5]
lb = [25, 40]
ub = [Inf, 40]
MathProgBase.loadproblem!(m, 4, 2, l, u, lb, ub, :Min, HS071())
MathProgBase.setwarmstart!(m,[0,15,15,11])
MathProgBase.optimize!(m)
@test MathProgBase.status(m) == :UserLimit

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ e6f1768c-d881-4ec8-acad-caa2bcd01aad
begin
	using Pkg
	Pkg.activate(".")
	using PlutoUI
	PlutoUI.TableOfContents(title = "Chapter 1",depth = 6)
end

# ╔═╡ 2633c520-5076-49e8-be7c-491c65da2665
using MLDatasets

# ╔═╡ 477e8407-a0aa-437a-9390-71977483c30a
md"""
# Machine Learning: A Probabilistic Perspective (MLPP)
## Chapter One: Introduction
Chapter one covers a broad sweep of machine learning. The chapter details a sample of the main points covered in the book. The first chapter ends in 3 exercises. Each exercise requires the use of Matlab code. We will convert these problems into Julia code.
"""

# ╔═╡ 4030ee45-be0e-4718-89a9-99bc687b6449
md"""
### Set up notebook
1. Change to root directory of project.
2. Activate the project.
3. Load packages
"""

# ╔═╡ 5faece8a-8d6d-11ec-3651-1d3f197a98e9
cd("/home/pushinglimits/Projects/MLPPProblemsInJulia")

# ╔═╡ 73d41226-84ff-45f9-bb5f-48762722ca67
md"""
## Details about chapter
"""

# ╔═╡ cbe518bd-b39e-4751-8c0b-853f01114fe4

# ╔═╡ 285b3f09-6ca7-4947-b4e4-00e137eaaa92
md"""
## Exercises
### Exercise 1.1

"""

# ╔═╡ cdc6563e-f717-49d2-92e1-ac7bd1e77e26

# ╔═╡ d1aa0ce8-0a41-48ce-8d93-6f4baa537231
md"""
### Exercise 1.2
"""

# ╔═╡ ac909932-d184-417c-b748-c41c937df477

# ╔═╡ 0baba8bf-64a3-4584-82d1-61703366511f
md"""
### Exercise 1.3
"""

# ╔═╡ 1435d6f2-98a0-437d-b41e-94126e59e46d

# ╔═╡ cf4c20b0-d70c-428c-bbda-a5ddeec613fc
#train_x, train_y = MNIST.traindata()

# ╔═╡ 51c6c59a-67fe-4b56-a0ce-f31e005dfebe
#test_x,  test_y  = MNIST.testdata()

# ╔═╡ 868dfc4c-2a7d-4fdc-901c-9f42a5c6106b

# ╔═╡ Cell order:
# ╟─477e8407-a0aa-437a-9390-71977483c30a
# ╠═4030ee45-be0e-4718-89a9-99bc687b6449
# ╠═5faece8a-8d6d-11ec-3651-1d3f197a98e9
# ╠═e6f1768c-d881-4ec8-acad-caa2bcd01aad
# ╠═2633c520-5076-49e8-be7c-491c65da2665
# ╟─73d41226-84ff-45f9-bb5f-48762722ca67
# ╠═cbe518bd-b39e-4751-8c0b-853f01114fe4
# ╟─285b3f09-6ca7-4947-b4e4-00e137eaaa92
# ╠═cdc6563e-f717-49d2-92e1-ac7bd1e77e26
# ╟─d1aa0ce8-0a41-48ce-8d93-6f4baa537231
# ╠═ac909932-d184-417c-b748-c41c937df477
# ╟─0baba8bf-64a3-4584-82d1-61703366511f
# ╠═1435d6f2-98a0-437d-b41e-94126e59e46d
# ╠═cf4c20b0-d70c-428c-bbda-a5ddeec613fc
# ╠═51c6c59a-67fe-4b56-a0ce-f31e005dfebe
# ╠═868dfc4c-2a7d-4fdc-901c-9f42a5c6106b

max_stars_repo_name: fieldofnodes/MLPPProblemsInJulia

max_issues_repo_name: fieldofnodes/MLPPProblemsInJulia

max_forks_repo_name: fieldofnodes/MLPPProblemsInJulia


text:
# * Derivatives

function diffop!(dest::AbstractMatrix,
                 L::BSplineOrRestricted, R::BSplineOrRestricted,
                 a, b, op=I)
    assert_compatible_bases(L,R)

    k = max(order(L),order(R))
    (bandwidth(dest,1) ≤ k-1) && (bandwidth(dest,2) ≤ k-1) ||
        throw(DimensionMismatch("Insufficient bandwidths of destination matrix"))

    # ∂' = -∂ ⇒ for weak Laplacians (in case of Dirichlet0 boundary
    # conditions), we negate the left basis.
    χ = (iseven(a) ? 1 : -1)*basis_functions(L, a)
    ξ = basis_functions(R, b)

    overlap_matrix!(dest, χ, op*ξ, weights(parent(R)))

    dest
end

diffop!(dest::BandedMatrix, B::BSplineOrRestricted, o) =
    diffop!(dest, B, B, o)

@materialize function *(Ac::AdjointBasisOrRestricted{<:BSpline},
                        D::Derivative,
                        B::BasisOrRestricted{<:BSpline})
    T -> begin
        Matrix(undef, parent(Ac), B, T)
    end
    dest::AbstractMatrix{T} -> begin
        diffop!(dest, parent(Ac), B, 0, 1)
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:BSpline},
                        Dc::QuasiAdjoint{<:Any,<:Derivative},
                        D::Derivative,
                        B::BasisOrRestricted{<:BSpline})
    T -> begin
        Matrix(undef, parent(Ac), B, T)
    end
    dest::AbstractMatrix{T} -> begin
        diffop!(dest, parent(Ac), B, 1, 1)
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:BSpline},
                        Dc::QuasiAdjoint{<:Any,<:Derivative},
                        O::QuasiDiagonal,
                        D::Derivative,
                        B::BasisOrRestricted{<:BSpline})
    T -> begin
        Matrix(undef, parent(Ac), B, T)
    end
    dest::AbstractMatrix{T} -> begin
        # Evaluate the quasi-diagonal operator O on the quadrature roots
        # of B.
        op = Diagonal(getindex.(Ref(O.diag), locs(parent(B))))
        diffop!(dest, parent(Ac), B, 1, 1, op)
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:BSpline},
                        Dc::QuasiAdjoint{<:Any,<:Derivative},
                        O::QuasiDiagonal,
                        P::QuasiDiagonal,
                        D::Derivative,
                        B::BasisOrRestricted{<:BSpline})
    T -> begin
        Matrix(undef, parent(Ac), B, T)
    end
    dest::AbstractMatrix{T} -> begin
        # Evaluate the quasi-diagonal operators O & P on the
        # quadrature roots of B.
        op = Diagonal(getindex.(Ref(O.diag), locs(parent(B))) .*
                      getindex.(Ref(P.diag), locs(parent(B))))
        diffop!(dest, parent(Ac), B, 1, 1, op)
    end
end


text:
module SwJ1644

#Example
export greet
greet() = print("Hello World!")

# input and output
using CSV
#using FileIO

export ReadData
function ReadData(filename::String)
    data = CSV.read(filename)
    #data = load(filename)
    return data
end

# To find peaks from the light curve
export peaks
function peaks(t,y)
    N = length(t)
    j = 0
    Peak = ones(N,2)
    for i=4:N-3
        if (max(y[i+1],y[i+2],y[i+3])<y[i])
            if (max(y[i-1],y[i-2],y[i-3])<y[i])
                j=j+1
                Peak[j,1]=t[i]
                Peak[j,2]=y[i]
            end
        end
#        if j>3
#            if (Peak[j,1]-Peak[j-1,1])<2*(Peak[j-1,1]-Peak[j-2,1])
#                j=j-1
#            end
#        end
    end
    Peak2 = ones(j,2)
    for i=1:j
        Peak2[i,1] = Peak[i,1]
        Peak2[i,2] = Peak[i,2]
    end
    return Peak2
end

# another method to find the peaks and also the dips
export peaks2
function peaks2(t,y)
    N = length(t)
    jPeak = 1
    jDip = 1
    ratio = 2.5
    flag = "Peak"
    Peak = ones(N,2) # 1 for t, 2 for y
    # Dip = Peak # Be careful, this makes these two are equal everywhere
    Dip = ones(N,2)
    PeakTmp = ones(N,2)
    DipTmp = ones(N,2)
    PeakTmp[1,1:2] = [t[1],y[1]]
    DipTmp[1,:] = [t[1],y[1]]
    for i in 2:N
        if y[i] > y[i-1]
            if y[i] > PeakTmp[jPeak,2]
                PeakTmp[jPeak,:] = [t[i],y[i]]
            end
            if flag == "Peak"
                if DipTmp[jDip,1] <= PeakTmp[jPeak,1]
                    DipTmp[jDip,:] = [t[i],y[i]]
                end
            #end
            elseif flag == "Dip"
                if PeakTmp[jPeak,2]/ratio > DipTmp[jDip,2]
                    Dip[jDip,1:2] = DipTmp[jDip,1:2]
                    jDip = jDip + 1
                    flag = "Peak"
                    DipTmp[jDip,:] = [t[i],y[i]]
                end
            end
        end
        if y[i] < y[i-1]
            if y[i] < DipTmp[jDip,2]
                DipTmp[jDip,:] = [t[i],y[i]]
            end
            if flag == "Dip"
                if PeakTmp[jPeak,1] <= DipTmp[jDip,1]
                    PeakTmp[jPeak,:] = [t[i],y[i]]
                end
            #end
            elseif flag == "Peak"
                #println(i, " Now in Peak 2, PeakTmp: ", PeakTmp[jPeak,2],", DipTmp:", DipTmp[jDip,2])
                if PeakTmp[jPeak,2]/ratio > DipTmp[jDip,2]
                    Peak[jPeak,1:2] = PeakTmp[jPeak,1:2]
                    jPeak = jPeak + 1
                    flag = "Dip"
                    PeakTmp[jPeak,:] = [t[i],y[i]]
                end
            end
        end
    end
    # cut the ends of the Peaks and Dips
    Peak2 = ones(jPeak-1,2)
    Dip2 = ones(jPeak-1,2)
    for i in 1:jPeak-1
        Peak2[i,:] = Peak[i,:]
        Dip2[i,:] = Dip[i,:]
    end
    return Peak2,Dip2
end

# plot
#using Plots
# export plotlc
function plotlc(x::Array,y::Array)
    plot(x,y)
end

end # module


text:
using Documenter, Luxor

makedocs(
    modules = [Luxor],
    sitename = "Luxor",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages    = Any[
        "Introduction to Luxor"   => "index.md",
        "A few examples"          => "examples.md",
        "Tutorial"                => "tutorial.md",
        "Basic concepts"          => "basics.md",
        "Simple shapes"           => "simplegraphics.md",
        "More graphics"           => "moregraphics.md",
        "Geometry tools"          => "geometrytools.md",
        "Tables and grids"        => "tables-grids.md",
        "Colors and styles"       => "colors-styles.md",
        "Polygons and paths"      => "polygons.md",
        "Text"                    => "text.md",
        "Transforms and matrices" => "transforms.md",
        "Clipping"                => "clipping.md",
        "Images"                  => "images.md",
        "Turtle graphics"         => "turtle.md",
        "Animation"               => "animation.md",
        "More examples"           => "moreexamples.md",
        "Index"                   => "functionindex.md"
        ]
    )

deploydocs(
    repo = "github.com/JuliaGraphics/Luxor.jl.git",
    target = "build"
)


text:
using JuMP
import MathOptInterface
const MOI = MathOptInterface
include("utils.jl")

function add_continuous_variables(m, p)
    @variables m begin
	    dvSize[p.Tech] >= 0     #X^{\sigma}_{t}: System Size of Technology t [kW]   (NEW)
    	dvSystemSizeSegment[p.Tech, p.Subdivision, p.Seg] >= 0   #X^{\sigma s}_{tks}: System size of technology t allocated to segmentation k, segment s [kW]  (NEW)
		dvGridPurchase[p.PricingTier, p.TimeStep] >= 0   # X^{g}_{uh}: Power from grid dispatched to meet electrical load in demand tier u during time step h [kW]  (NEW)
	    dvRatedProduction[p.Tech, p.TimeStep] >= 0   #X^{rp}_{th}: Rated production of technology t during time step h [kW]  (NEW)
	    dvProductionToGrid[p.Tech, p.SalesTiers, p.TimeStep] >= 0  # X^{ptg}_{tuh}: Exports from electrical production to the grid by technology t in demand tier u during time step h [kW]   (NEW)
	    dvStorageToGrid[p.StorageSalesTiers, p.TimeStep] >= 0  # X^{stg}_{uh}: Exports from electrical storage to the grid in demand tier u during time step h [kW]  (NEW)
		dvProductionToStorage[p.Storage, p.Tech, p.TimeStep] >= 0  # X^{ptg}_{bth}: Power from technology t used to charge storage system b during time step h [kW]  (NEW)
	    dvDischargeFromStorage[p.Storage, p.TimeStep] >= 0 # X^{pts}_{bh}: Power discharged from storage system b during time step h [kW]  (NEW)
	    dvGridToStorage[p.TimeStep] >= 0 # X^{gts}_{h}: Electrical power delivered to storage by the grid in time step h [kW]  (NEW)
	    dvStorageSOC[p.Storage, p.TimeStepBat] >= 0  # X^{se}_{bh}: State of charge of storage system b in time step h   (NEW)
	    dvStorageCapPower[p.Storage] >= 0   # X^{bkW}_b: Power capacity of storage system b [kW]  (NEW)
	    dvStorageCapEnergy[p.Storage] >= 0   # X^{bkWh}_b: Energy capacity of storage system b [kWh]  (NEW)
	    dvProdIncent[p.Tech] >= 0   # X^{pi}_{t}: Production incentive collected for technology [$]
		dvPeakDemandE[p.Ratchets, p.DemandBin] >= 0  # X^{de}_{re}:  Peak electrical power demand allocated to tier e during ratchet r [kW]
		dvPeakDemandEMonth[p.Month, p.DemandMonthsBin] >= 0  #  X^{dn}_{mn}: Peak electrical power demand allocated to tier n during month m [kW]
		dvPeakDemandELookback[p.Month] >= 0  # X^{lp}: Peak electric demand look back [kW]
        MinChargeAdder >= 0   #to be removed
		#UtilityMinChargeAdder[p.Month] >= 0   #X^{mc}_m: Annual utility minimum charge adder in month m [\$]
		#CHP and Fuel-burning variables
		dvFuelUsage[p.Tech, p.TimeStep] >= 0  # Fuel burned by technology t in time step h
		#dvFuelBurnYIntercept[p.Tech, p.TimeStep]  #X^{fb}_{th}: Y-intercept of fuel burned by technology t in time step h
		#dvThermalProduction[p.Tech, p.TimeStep]  #X^{tp}_{th}: Thermal production by technology t in time step h
		#dvAbsorptionChillerDemand[p.TimeStep]  #X^{ac}_h: Thermal power consumption by absorption chiller in time step h
		#dvElectricChillerDemand[p.TimeStep]  #X^{ec}_h: Electrical power consumption by electric chiller in time step h
    end
end

function add_integer_variables(m, p)
    @variables m begin
        binNMLorIL[p.NMILRegime], Bin    # Z^{nmil}_{v}: 1 If generation is in net metering interconnect limit regime v; 0 otherwise
        binProdIncent[p.Tech], Bin # Z^{pi}_{t}:  1 If production incentive is available for technology t; 0 otherwise
		binSegmentSelect[p.Tech, p.Subdivision, p.Seg], Bin # Z^{\sigma s}_{tks} 1 if technology t, segmentation k is in segment s; 0 ow. (NEW)
        binSingleBasicTech[p.Tech,p.TechClass], Bin   #  Z^\text{sbt}_{tc}: 1 If technology t is used for technology class c; 0 otherwise
        binTechIsOnInTS[p.Tech, p.TimeStep], Bin  # 1 Z^{to}_{th}: If technology t is operating in time step h; 0 otherwise
		binDemandTier[p.Ratchets, p.DemandBin], Bin  # 1 If tier e has allocated demand during ratchet r; 0 otherwise
        binDemandMonthsTier[p.Month, p.DemandMonthsBin], Bin # 1 If tier n has allocated demand during month m; 0 otherwise
		binEnergyTier[p.Month, p.PricingTier], Bin    #  Z^{ut}_{mu} 1 If demand tier $u$ is active in month m; 0 otherwise (NEW)
    end
end

function add_parameters(m, p)
    m[:r_tax_fraction_owner] = (1 - p.r_tax_owner)
    m[:r_tax_fraction_offtaker] = (1 - p.r_tax_offtaker)

    m[:PVTechs] = filter(t->startswith(t, "PV"), p.Tech)
	if !isempty(p.Tech)
		m[:WindTechs] = p.TechsInClass["WIND"]
		m[:GeneratorTechs] = p.TechsInClass["GENERATOR"]
	else
		m[:WindTechs] = []
		m[:GeneratorTechs] = []
	end
end

function add_cost_expressions(m, p)
	m[:TotalTechCapCosts] = @expression(m, p.two_party_factor * (
		sum( p.CapCostSlope[t,s] * m[:dvSystemSizeSegment][t,"CapCost",s] for t in p.Tech, s in 1:p.SegByTechSubdivision["CapCost",t] ) +
		sum( p.CapCostYInt[t,s] * m[:binSegmentSelect][t,"CapCost",s] for t in p.Tech, s in 1:p.SegByTechSubdivision["CapCost",t] )
	))
	m[:TotalStorageCapCosts] = @expression(m, p.two_party_factor *
		sum( p.StorageCostPerKW[b]*m[:dvStorageCapPower][b] + p.StorageCostPerKWH[b]*m[:dvStorageCapEnergy][b] for b in p.Storage )
	)
	m[:TotalPerUnitSizeOMCosts] = @expression(m, p.two_party_factor * p.pwf_om *
		sum( p.OMperUnitSize[t] * m[:dvSize][t] for t in p.Tech )
	)
    if !isempty(m[:GeneratorTechs])
		m[:TotalPerUnitProdOMCosts] = @expression(m, p.two_party_factor * p.pwf_om *
			sum( p.OMcostPerUnitProd[t] * m[:dvRatedProduction][t,ts] for t in p.FuelBurningTechs, ts in p.TimeStep )
		)
    else
        m[:TotalPerUnitProdOMCosts] = @expression(m, 0.0)
	end
end

function add_export_expressions(m, p)
	if !isempty(p.Tech)
		# NOTE: LevelizationFactor is baked into m[:dvProductionToGrid]
		m[:TotalExportBenefit] = @expression(m, p.pwf_e * p.TimeStepScaling * sum(
			sum(p.GridExportRates[u,ts] * m[:dvStorageToGrid][u,ts] for u in p.StorageSalesTiers)
			+ sum(p.GridExportRates[u,ts] * m[:dvProductionToGrid][t,u,ts]
				  for u in p.SalesTiers, t in p.TechsBySalesTier[u]
			) for ts in p.TimeStep )
		)
		m[:CurtailedElecWIND] = @expression(m,
			p.TimeStepScaling * sum(m[:dvProductionToGrid][t,u,ts]
				for t in m[:WindTechs], u in p.CurtailmentTiers, ts in p.TimeStep)
		)
		m[:ExportedElecWIND] = @expression(m,
			p.TimeStepScaling * sum(m[:dvProductionToGrid][t,u,ts]
				for t in m[:WindTechs], u in p.SalesTiersByTech[t], ts in p.TimeStep) - m[:CurtailedElecWIND]
		)
		m[:ExportedElecGEN] = @expression(m,
			p.TimeStepScaling * sum(m[:dvProductionToGrid][t,u,ts]
				for t in m[:GeneratorTechs], u in p.SalesTiersByTech[t], ts in p.TimeStep)
		)
		m[:ExportBenefitYr1] = @expression(m,
			p.TimeStepScaling * sum(
			sum( p.GridExportRates[u,ts] * m[:dvStorageToGrid][u,ts] for u in p.StorageSalesTiers)
			+ sum( p.GridExportRates[u,ts] * m[:dvProductionToGrid][t,u,ts]
				for u in p.SalesTiers, t in p.TechsBySalesTier[u])
			for ts in p.TimeStep )
		)
	else
		m[:TotalExportBenefit] = 0
		m[:CurtailedElecWIND] = 0
		m[:ExportedElecWIND] = 0
		m[:ExportedElecGEN] = 0
		m[:ExportBenefitYr1] = 0
	end
end

function add_bigM_adjustments(m, p)
	m[:NewMaxUsageInTier] = Array{Float64,2}(undef,12, p.PricingTierCount+1)
	m[:NewMaxDemandInTier] = Array{Float64,2}(undef, length(p.Ratchets), p.DemandBinCount)
	m[:NewMaxDemandMonthsInTier] = Array{Float64,2}(undef,12, p.DemandMonthsBinCount)
	m[:NewMaxSize] = Dict()

	NewMaxSizeByHour = Array{Float64,2}(undef,length(p.Tech),p.TimeStepCount)
	# m[:NewMaxDemandMonthsInTier] sets a new minimum if the new peak demand for the month, minus the size of all previous bins, is less than the existing bin size.
	if !isempty(p.ElecStorage)
		added_power = p.StorageMaxSizePower["Elec"]
		added_energy = p.StorageMaxSizeEnergy["Elec"]
	else
		added_power = 1.0e-3
		added_energy = 1.0e-3
	end

	for n in p.DemandMonthsBin
		for mth in p.Month
			if n > 1
				m[:NewMaxDemandMonthsInTier][mth,n] = minimum([p.MaxDemandMonthsInTier[n],
					added_power + maximum([p.ElecLoad[ts] #+ LoadProfileChillerElectric[ts]
					for ts in p.TimeStepRatchetsMonth[mth]])  -
					sum(m[:NewMaxDemandMonthsInTier][mth,np] for np in 1:(n-1)) ]
				)
			else
				m[:NewMaxDemandMonthsInTier][mth,n] = minimum([p.MaxDemandMonthsInTier[n],
					added_power + maximum([p.ElecLoad[ts] #+ LoadProfileChillerElectric[ts]
					for ts in p.TimeStepRatchetsMonth[mth]])   ])
			end
		end
	end

	# m[:NewMaxDemandInTier] sets a new minimum if the new peak demand for the ratchet, minus the size of all previous bins for the ratchet, is less than the existing bin size.
	for e in p.DemandBin
		for r in p.Ratchets
			if e > 1
				m[:NewMaxDemandInTier][r,e] = minimum([p.MaxDemandInTier[e],
				added_power + maximum([p.ElecLoad[ts] #+ p.LoadProfileChillerElectric[ts]
					for ts in p.TimeStep])  -
				sum(m[:NewMaxDemandInTier][r,ep] for ep in 1:(e-1))
				])
			else
				m[:NewMaxDemandInTier][r,e] = minimum([p.MaxDemandInTier[e],
				added_power + maximum([p.ElecLoad[ts] #+ p.LoadProfileChillerElectric[ts]
					for ts in p.TimeStep])
				])
			end
		end
	end

	# m[:NewMaxUsageInTier] sets a new minumum if the total demand for the month, minus the size of all previous bins, is less than the existing bin size.
	for u in p.PricingTier
		for mth in p.Month
			if u > 1
				m[:NewMaxUsageInTier][mth,u] = minimum([p.MaxUsageInTier[u],
					added_energy + sum(p.ElecLoad[ts] #+ p.LoadProfileChillerElectric[ts]
					for ts in p.TimeStepRatchetsMonth[mth]) - sum(m[:NewMaxUsageInTier][mth,up] for up in 1:(u-1))
				])
			else
				m[:NewMaxUsageInTier][mth,u] = minimum([p.MaxUsageInTier[u],
					added_energy + sum(p.ElecLoad[ts] #+ p.LoadProfileChillerElectric[ts]
					for ts in p.TimeStepRatchetsMonth[mth])
				])
			end
		end
	end

	# NewMaxSize generates a new maximum size that is equal to the largest monthly load of the year.  This is intended to be a reasonable upper bound on size that would never be exceeeded, but is sufficiently small to replace much larger big-M values placed as a default.
	TempHeatingTechs = [] #temporarily replace p.HeatingTechs which is undefined
	TempCoolingTechs = [] #temporarily replace p.CoolingTechs which is undefined

	for t in TempHeatingTechs
		m[:NewMaxSize][t] = maximum([sum(p.HeatingLoad[ts] for ts in p.TimeStepRatchetsMonth[mth]) for mth in p.Month])
		if (m[:NewMaxSize][t] > p.MaxSize[t])
			m[:NewMaxSize][t] = p.MaxSize[t]
		end
	end
	for t in TempCoolingTechs
		m[:NewMaxSize][t] = maximum([sum(p.CoolingLoad[ts] for ts in p.TimeStepRatchetsMonth[mth]) for mth in p.Month])
		if (m[:NewMaxSize][t] > p.MaxSize[t])
			m[:NewMaxSize][t] = p.MaxSize[t]
		end
	end
	for t in p.ElectricTechs
		m[:NewMaxSize][t] = maximum([sum(p.ElecLoad[ts] for ts in p.TimeStepRatchetsMonth[mth]) for mth in p.Month])
		if (m[:NewMaxSize][t] > p.MaxSize[t])
			m[:NewMaxSize][t] = p.MaxSize[t]
		end
	end

	# NewMaxSizeByHour is designed to scale the right-hand side of the constraint limiting rated production in each hour to the production factor; in most cases this is unaffected unless the production factor is zero, in which case the right-hand side is set to zero.
	#for t in p.ElectricTechs
	#	for ts in p.TimeStep
	#		NewMaxSizeByHour[t,ts] = minimum([m[:NewMaxSize][t],
	#			sum(p.ProdFactor[t,d,ts] for d in p.Load if p.LoadProfile[d,ts] > 0)  * m[:NewMaxSize][t],
	#			sum(p.LoadProfile[d,ts] for d in ["1R"], ts in p.TimeStep)
	#		])
	#	end
	#end
end

function add_no_grid_constraints(m, p)
	for ts in p.TimeStepsWithoutGrid
		fix(m[:dvGridToStorage][ts], 0.0, force=true)
		for u in p.PricingTier
			fix(m[:dvGridPurchase][u,ts], 0.0, force=true)
		end
		#for u in p.StorageSalesTiers  #don't allow curtailment of storage either
		#	fix(m[:dvStorageToGrid][u,ts], 0.0, force=true)
		#end
	end
end

function add_fuel_constraints(m, p)

	##Constraint (1a): Sum of fuel used must not exceed prespecified limits
	@constraint(m, TotalFuelConsumptionCon[f in p.FuelType],
		sum( m[:dvFuelUsage][t,ts] for t in p.TechsByFuelType[f], ts in p.TimeStep ) <=
		p.FuelLimit[f]
	)

	# Constraint (1b): Fuel burn for non-CHP Constraints
	@constraint(m, FuelBurnCon[t in p.FuelBurningTechs, ts in p.TimeStep],
		m[:dvFuelUsage][t,ts]  == (p.FuelBurnSlope[t] * p.ProductionFactor[t,ts] * m[:dvRatedProduction][t,ts]) +
			(p.FuelBurnYInt[t] * m[:binTechIsOnInTS][t,ts])
	)
	m[:TotalGenFuelCharges] = @expression(m, p.pwf_e * p.TimeStepScaling * sum( p.FuelCost[f] *
		sum(m[:dvFuelUsage][t,ts] for t in p.TechsByFuelType[f], ts in p.TimeStep)
		for f in p.FuelType)
	)
	#if !isempty(CHPTechs)
		#Constraint (1c): Total Fuel burn for CHP
		#@constraint(m, CHPFuelBurnCon[t in CHPTechs, ts in p.TimeStep],
		#			m[:dvFuelUsage][t,ts]  == p.FuelBurnAmbientFactor[t,ts] * (dvFuelBurnYIntercept[t,th] +
		#				p.ProductionFactor[t,ts] * p.FuelBurnRateM[t] * m[:dvRatedProduction][t,ts])
		#			)

		#Constraint (1d): Y-intercept fuel burn for CHP
		#@constraint(m, CHPFuelBurnYIntCon[t in CHPTechs, ts in p.TimeStep],
		#			p.FuelBurnYIntRate[t] * m[:dvSize][t] - m[:NewMaxSize][t] * (1-m[:binTechIsOnInTS][t,ts])  <= dvFuelBurnYIntercept[t,th]
		#			)
	#end

	#if !isempty(NonCHPHeatingTechs)
		#Constraint (1e): Total Fuel burn for Boiler
		#@constraint(m, BoilerFuelBurnCon[t in NonCHPHeatingTechs, ts in p.TimeStep],
		#			m[:dvFuelUsage][t,ts]  ==  dvThermalProduction[t,ts] / p.BoilerEfficiency
		#			)
	#end

	### Constraint set (2): CHP Thermal Production Constraints
	#if !isempty(CHPTechs)
		#Constraint (2a-1): Upper Bounds on Thermal Production Y-Intercept
		#@constraint(m, CHPYInt2a1Con[t in CHPTechs, ts in p.TimeStep],
		#			dvThermalProductionYIntercept[t,ts] <= CHPThermalProdIntercept[t] * m[:dvSize][t]
		#			)
		# Constraint (2a-2): Upper Bounds on Thermal Production Y-Intercept
		#@constraint(m, CHPYInt2a1Con[t in CHPTechs, ts in p.TimeStep],
		#			dvThermalProductionYIntercept[t,ts] <= CHPThermalProdIntercept[t] * m[:NewMaxSize][t] * m[:binTechIsOnInTS][t,ts]
		#			)
		# Constraint (2b): Thermal Production of CHP
		#@constraint(m, CHPThermalProductionCpn[t in CHPTechs, ts in p.TimeStep],
		#			dvThermalProduction[t,ts] <=  HotWaterAmbientFactor[t,ts] * HotWaterThermalFactor[t,ts] * (
		#			CHPThermalProdSlope[t] * ProductionFactor[t,ts] * m[:dvRatedProduction][t,ts] + dvThermalProductionYIntercept[t,ts]
		#				)
		#			)
	#end
end

function add_binTechIsOnInTS_constraints(m, p)
	### Section 3: Switch Constraints
	#Constraint (3a): Technology must be on for nonnegative output (fuel-burning only)
	@constraint(m, ProduceIfOnCon[t in p.FuelBurningTechs, ts in p.TimeStep],
		m[:dvRatedProduction][t,ts] <= m[:NewMaxSize][t] * m[:binTechIsOnInTS][t,ts]
	)
	#Constraint (3b): Technologies that are turned on must not be turned down
	@constraint(m, MinTurndownCon[t in p.FuelBurningTechs, ts in p.TimeStep],
		p.MinTurndown[t] * m[:dvSize][t] - m[:dvRatedProduction][t,ts] <= m[:NewMaxSize][t] * (1-m[:binTechIsOnInTS][t,ts])
	)
end

function add_storage_size_constraints(m, p)
	# Constraint (4a): Reconcile initial state of charge for storage systems
	@constraint(m, InitStorageCon[b in p.Storage], m[:dvStorageSOC][b,0] == p.StorageInitSOC[b] * m[:dvStorageCapEnergy][b])
	# Constraint (4b)-1: Lower bound on Storage Energy Capacity
	@constraint(m, StorageEnergyLBCon[b in p.Storage], m[:dvStorageCapEnergy][b] >= p.StorageMinSizeEnergy[b])
	# Constraint (4b)-2: Upper bound on Storage Energy Capacity
	@constraint(m, StorageEnergyUBCon[b in p.Storage], m[:dvStorageCapEnergy][b] <= p.StorageMaxSizeEnergy[b])
	# Constraint (4c)-1: Lower bound on Storage Power Capacity
	@constraint(m, StoragePowerLBCon[b in p.Storage], m[:dvStorageCapPower][b] >= p.StorageMinSizePower[b])
	# Constraint (4c)-2: Upper bound on Storage Power Capacity
	@constraint(m, StoragePowerUBCon[b in p.Storage], m[:dvStorageCapPower][b] <= p.StorageMaxSizePower[b])
end

function add_storage_op_constraints(m, p)
	### Battery Operations
	# Constraint (4d): Electrical production sent to storage or grid must be less than technology's rated production
	@constraint(m, ElecTechProductionFlowCon[b in p.ElecStorage, t in p.ElectricTechs, ts in p.TimeStepsWithGrid],
		m[:dvProductionToStorage][b,t,ts] + sum(m[:dvProductionToGrid][t,u,ts] for u in p.SalesTiersByTech[t]) <=
		p.ProductionFactor[t,ts] * p.LevelizationFactor[t] * m[:dvRatedProduction][t,ts]
	)
	# Constraint (4e): Electrical production sent to storage or grid must be less than technology's rated production - no grid
	@constraint(m, ElecTechProductionFlowNoGridCon[b in p.ElecStorage, t in p.ElectricTechs, ts in p.TimeStepsWithoutGrid],
		m[:dvProductionToStorage][b,t,ts] + sum(m[:dvProductionToGrid][t,u,ts] for u in p.CurtailmentTiers)  <=
		p.ProductionFactor[t,ts] * p.LevelizationFactor[t] * m[:dvRatedProduction][t,ts]
	)
	# Constraint (4f)-1: (Hot) Thermal production sent to storage or grid must be less than technology's rated production
	#@constraint(m, HeatingTechProductionFlowCon[b in p.HotTES, t in p.HeatingTechs, ts in p.TimeStep],
    #	        m[:dvProductionToStorage][b,t,ts]  <=
	#			p.ProductionFactor[t,ts] * dvThermalProduction[t,ts]
	#			)
	# Constraint (4f)-2: (Cold) Thermal production sent to storage or grid must be less than technology's rated production
	#@constraint(m, CoolingTechProductionFlowCon[b in p.ColdTES, t in p.CoolingTechs, ts in p.TimeStep],
    #	        m[:dvProductionToStorage][b,t,ts]  <=
	#			p.ProductionFactor[t,ts] * dvThermalProduction[t,ts]
	#			)
	# Constraint (4g): Reconcile state-of-charge for electrical storage - with grid
	@constraint(m, ElecStorageInventoryCon[b in p.ElecStorage, ts in p.TimeStepsWithGrid],
		m[:dvStorageSOC][b,ts] == m[:dvStorageSOC][b,ts-1] + p.TimeStepScaling * (
			sum(p.ChargeEfficiency[t,b] * m[:dvProductionToStorage][b,t,ts] for t in p.ElectricTechs) +
			p.GridChargeEfficiency*m[:dvGridToStorage][ts] - m[:dvDischargeFromStorage][b,ts]/p.DischargeEfficiency[b]
		)
	)

	# Constraint (4h): Reconcile state-of-charge for electrical storage - no grid
	@constraint(m, ElecStorageInventoryConNoGrid[b in p.ElecStorage, ts in p.TimeStepsWithoutGrid],
		m[:dvStorageSOC][b,ts] == m[:dvStorageSOC][b,ts-1] + p.TimeStepScaling * (
			sum(p.ChargeEfficiency[t,b] * m[:dvProductionToStorage][b,t,ts] for t in p.ElectricTechs) - m[:dvDischargeFromStorage][b,ts]/p.DischargeEfficiency[b]
		)
	)

	# Constraint (4i)-1: Reconcile state-of-charge for (hot) thermal storage
	#@constraint(m, HotTESInventoryCon[b in p.HotTES, ts in p.TimeStep],
    #	        m[:dvStorageSOC][b,ts] == m[:dvStorageSOC][b,ts-1] + p.TimeStepScaling * (
	#				sum(p.ChargeEfficiency[t,b] * m[:dvProductionToStorage][b,t,ts] for t in p.HeatingTechs) -
	#				m[:dvDischargeFromStorage][b,ts]/p.DischargeEfficiency[b]
	#				)
	#			)

	# Constraint (4i)-2: Reconcile state-of-charge for (cold) thermal storage
	#@constraint(m, ColdTESInventoryCon[b in p.ColdTES, ts in p.TimeStep],
    #	        m[:dvStorageSOC][b,ts] == m[:dvStorageSOC][b,ts-1] + p.TimeStepScaling * (
	#				sum(p.ChargeEfficiency[t,b] * m[:dvProductionToStorage][b,t,ts] for t in p.CoolingTechs) -
	#				m[:dvDischargeFromStorage][b,ts]/p.DischargeEfficiency[b]
	#				)
	#			)

	# Constraint (4j): Minimum state of charge
	@constraint(m, MinStorageLevelCon[b in p.Storage, ts in p.TimeStep],
		m[:dvStorageSOC][b,ts] >= p.StorageMinSOC[b] * m[:dvStorageCapEnergy][b]
	)

	#Constraint (4i)-1: Dispatch to electrical storage is no greater than power capacity
	@constraint(m, ElecChargeLEQCapCon[b in p.ElecStorage, ts in p.TimeStep],
		m[:dvStorageCapPower][b] >= (
			sum(m[:dvProductionToStorage][b,t,ts] for t in p.ElectricTechs) + m[:dvGridToStorage][ts]
		)
	)

	#Constraint (4i)-2: Dispatch to hot storage is no greater than power capacity
	#@constraint(m, HotTESChargeLEQCapCon[b in p.HotTES, ts in p.TimeStep],
    #	        m[:dvStorageCapPower][b] >= (
	#				sum(m[:dvProductionToStorage][b,t,ts] for t in p.HeatingTechs)
	#				)
	#			)

	#Constraint (4i)-3: Dispatch to cold storage is no greater than power capacity
	#@constraint(m, ColdTESChargeLEQCapCon[b in p.ColdTES, ts in p.TimeStep],
    #	        m[:dvStorageCapPower][b] >= (
	#				sum(m[:dvProductionToStorage][b,t,ts] for t in p.CoolingTechs)
	#				)
	#			)

	#Constraint (4j): Dispatch from storage is no greater than power capacity
	@constraint(m, DischargeLEQCapCon[b in p.Storage, ts in p.TimeStep],
		m[:dvStorageCapPower][b] >= m[:dvDischargeFromStorage][b,ts]
	)

	#Constraint (4k)-alt: Dispatch to and from electrical storage is no greater than power capacity
	@constraint(m, ElecChargeLEQCapConAlt[b in p.ElecStorage, ts in p.TimeStepsWithGrid],
		m[:dvStorageCapPower][b] >=   m[:dvDischargeFromStorage][b,ts] +
			sum(m[:dvProductionToStorage][b,t,ts] for t in p.ElectricTechs) + m[:dvGridToStorage][ts]
	)
	#Constraint (4l)-alt: Dispatch from electrical storage is no greater than power capacity
	@constraint(m, DischargeLEQCapConNoGridAlt[b in p.ElecStorage, ts in p.TimeStepsWithoutGrid],
		m[:dvStorageCapPower][b] >= m[:dvDischargeFromStorage][b,ts] +
			sum(m[:dvProductionToStorage][b,t,ts] for t in p.ElectricTechs)
	)

	#Constraint (4m)-1: Dispatch from thermal storage is no greater than power capacity
	#@constraint(m, DischargeLEQCapCon[b in p.HotTES, ts in p.TimeStep],
    #	        m[:dvStorageCapPower][b] >= sum(m[:dvProductionToStorage][b,t,ts] for t in p.HeatingTechs)
	#			)
	#Constraint (4m)-2: Dispatch from thermal storage is no greater than power capacity
	#@constraint(m, DischargeLEQCapCon[b in p.ColdTES, ts in p.TimeStep],
    #	        m[:dvStorageCapPower][b] >= sum(m[:dvProductionToStorage][b,t,ts] for t in p.CoolingTechs)
	#			)

	#Constraint (4n): State of charge upper bound is storage system size
	@constraint(m, StorageEnergyMaxCapCon[b in p.Storage, ts in p.TimeStep],
		m[:dvStorageSOC][b,ts] <= m[:dvStorageCapEnergy][b]
	)

	if !p.StorageCanGridCharge
		for ts in p.TimeStepsWithGrid
			fix(m[:dvGridToStorage][ts], 0.0, force=true)
		end
	end
end

function add_prod_incent_constraints(m, p)
	##Constraint (6a)-1: Production Incentive Upper Bound (unchanged)
	@constraint(m, ProdIncentUBCon[t in p.Tech],
		m[:dvProdIncent][t] <= m[:binProdIncent][t] * p.MaxProdIncent[t] * p.pwf_prod_incent[t] * p.two_party_factor)
	##Constraint (6a)-2: Production Incentive According to Production (updated)
	@constraint(m, IncentByProductionCon[t in p.Tech],
		m[:dvProdIncent][t] <= p.TimeStepScaling * p.ProductionIncentiveRate[t] * p.pwf_prod_incent[t] * p.two_party_factor *
			sum(p.ProductionFactor[t, ts] * m[:dvRatedProduction][t,ts] for ts in p.TimeStep)
	)
	##Constraint (6b): System size max to achieve production incentive
	@constraint(m, IncentBySystemSizeCon[t in p.Tech],
		m[:dvSize][t]  <= p.MaxSizeForProdIncent[t] + m[:NewMaxSize][t] * (1 - m[:binProdIncent][t]))

	if !isempty(p.Tech)
		m[:TotalProductionIncentive] = @expression(m, sum(m[:dvProdIncent][t] for t in p.Tech))
	else
		m[:TotalProductionIncentive] = 0
	end
end

function add_tech_size_constraints(m, p)
	# PV techs can be constrained by space available based on location at site (roof, ground, both)
	@constraint(m, TechMaxSizeByLocCon[loc in p.Location],
		sum( m[:dvSize][t] * p.TechToLocation[t, loc] for t in p.Tech) <= p.MaxSizesLocation[loc]
	)

	#Constraint (7a): Single Basic Technology Constraints
	@constraint(m, TechMaxSizeByClassCon[c in p.TechClass, t in p.TechsInClass[c]],
		m[:dvSize][t] <= m[:NewMaxSize][t] * m[:binSingleBasicTech][t,c]
		)
	##Constraint (7b): At most one Single Basic Technology per Class
	@constraint(m, TechClassMinSelectCon[c in p.TechClass],
		sum( m[:binSingleBasicTech][t,c] for t in p.TechsInClass[c] ) <= 1
		)
	##Constraint (7c): Minimum size for each tech class
	@constraint(m, TechClassMinSizeCon[c in p.TechClass],
				sum( m[:dvSize][t] for t in p.TechsInClass[c] ) >= p.TechClassMinSize[c]
			)

	## Constraint (7d): Non-turndown technologies are always at rated production
	@constraint(m, RenewableRatedProductionCon[t in p.TechsNoTurndown, ts in p.TimeStep],
		m[:dvRatedProduction][t,ts] == m[:dvSize][t]
	)

	##Constraint (7e): Derate factor limits production variable (separate from ProductionFactor)
	@constraint(m, TurbineRatedProductionCon[t in p.Tech, ts in p.TimeStep; !(t in p.TechsNoTurndown)],
		m[:dvRatedProduction][t,ts]  <= p.ElectricDerate[t,ts] * m[:dvSize][t]
	)

	##Constraint (7f)-1: Minimum segment size
	@constraint(m, SegmentSizeMinCon[t in p.Tech, k in p.Subdivision, s in 1:p.SegByTechSubdivision[k,t]],
		m[:dvSystemSizeSegment][t,k,s]  >= p.SegmentMinSize[t,k,s] * m[:binSegmentSelect][t,k,s]
	)

	##Constraint (7f)-2: Maximum segment size
	@constraint(m, SegmentSizeMaxCon[t in p.Tech, k in p.Subdivision, s in 1:p.SegByTechSubdivision[k,t]],
		m[:dvSystemSizeSegment][t,k,s]  <= p.SegmentMaxSize[t,k,s] * m[:binSegmentSelect][t,k,s]
	)

	##Constraint (7g):  Segments add up to system size
	@constraint(m, SegmentSizeAddCon[t in p.Tech, k in p.Subdivision],
		sum(m[:dvSystemSizeSegment][t,k,s] for s in 1:p.SegByTechSubdivision[k,t]) == m[:dvSize][t]
	)

	##Constraint (7h): At most one segment allowed
	@constraint(m, SegmentSelectCon[c in p.TechClass, t in p.TechsInClass[c], k in p.Subdivision],
		sum(m[:binSegmentSelect][t,k,s] for s in 1:p.SegByTechSubdivision[k,t]) <= m[:binSingleBasicTech][t,c]
	)
end

function add_load_balance_constraints(m, p)
	@constraint(m, ElecLoadBalanceCon[ts in p.TimeStepsWithGrid],
		sum(p.ProductionFactor[t,ts] * p.LevelizationFactor[t] * m[:dvRatedProduction][t,ts] for t in p.ElectricTechs) +
		sum( m[:dvDischargeFromStorage][b,ts] for b in p.ElecStorage ) +
		sum( m[:dvGridPurchase][u,ts] for u in p.PricingTier ) ==
		sum( sum(m[:dvProductionToStorage][b,t,ts] for b in p.ElecStorage) +
			sum(m[:dvProductionToGrid][t,u,ts] for u in p.SalesTiersByTech[t]) for t in p.ElectricTechs) +
		sum(m[:dvStorageToGrid][u,ts] for u in p.StorageSalesTiers) + m[:dvGridToStorage][ts] +
		## sum(dvThermalProduction[t,ts] for t in p.CoolingTechs )/ p.ElectricChillerEfficiency +
		p.ElecLoad[ts]
	)

	##Constraint (8b): Electrical Load Balancing without Grid
	@constraint(m, ElecLoadBalanceNoGridCon[ts in p.TimeStepsWithoutGrid],
		sum(p.ProductionFactor[t,ts] * p.LevelizationFactor[t] * m[:dvRatedProduction][t,ts] for t in p.ElectricTechs) +
		sum( m[:dvDischargeFromStorage][b,ts] for b in p.ElecStorage )  ==
		sum( sum(m[:dvProductionToStorage][b,t,ts] for b in p.ElecStorage) +
			sum(m[:dvProductionToGrid][t,u,ts] for u in p.CurtailmentTiers) for t in p.ElectricTechs) +
		## sum(dvThermalProduction[t,ts] for t in p.CoolingTechs )/ p.ElectricChillerEfficiency +
		p.ElecLoad[ts]
	)
end

function add_storage_grid_constraints(m, p)
	##Constraint (8c): Grid-to-storage no greater than grid purchases
	@constraint(m, GridToStorageCon[ts in p.TimeStepsWithGrid],
		sum( m[:dvGridPurchase][u,ts] for u in p.PricingTier)  >= m[:dvGridToStorage][ts]
	)

	##Constraint (8d): Storage-to-grid no greater than discharge from Storage
	@constraint(m, StorageToGridCon[ts in p.TimeStepsWithGrid],
		sum( m[:dvDischargeFromStorage][b,ts] for b in p.ElecStorage)  >= sum(m[:dvStorageToGrid][u,ts] for u in p.StorageSalesTiers)
	)
end

function add_prod_grid_constraints(m, p)
	##Constraint (8e): Production-to-grid no greater than production
	@constraint(m, ProductionToGridCon[t in p.Tech, ts in p.TimeStepsWithGrid],
	 p.ProductionFactor[t,ts] * p.LevelizationFactor[t] * m[:dvRatedProduction][t,ts] >= sum(m[:dvProductionToGrid][t,u,ts] for u in p.SalesTiersByTech[t])
	)

	##Constraint (8f): Total sales to grid no greater than annual allocation - storage tiers
	@constraint(m,  AnnualGridSalesLimitCon,
	 p.TimeStepScaling * (
		sum( m[:dvStorageToGrid][u,ts] for u in p.StorageSalesTiers, ts in p.TimeStepsWithGrid if !(u in p.CurtailmentTiers)) +  sum(m[:dvProductionToGrid][t,u,ts] for u in p.SalesTiers, t in p.TechsBySalesTier[u], ts in p.TimeStepsWithGrid if !(u in p.CurtailmentTiers))) <= p.MaxGridSales[1]
	)
end

function add_NMIL_constraints(m, p)
	#Constraint (9a): exactly one net-metering regime must be selected
	@constraint(m, sum(m[:binNMLorIL][n] for n in p.NMILRegime) == 1)

	##Constraint (9b): Maximum system sizes for each net-metering regime
	@constraint(m, NetMeterSizeLimit[n in p.NMILRegime],
		sum(p.TurbineDerate[t] * m[:dvSize][t]
		for t in p.TechsByNMILRegime[n]) <= p.NMILLimits[n] * m[:binNMLorIL][n]
	)
end

function add_nem_constraint(m, p)
	# NEM is SalesTier 1
	# dvStorageToGrid is always fixed at 0.0, remove it?
	@constraint(m, GridSalesLimit,
		p.TimeStepScaling * sum(m[:dvProductionToGrid][t,1,ts] for t in p.TechsBySalesTier[1], ts in p.TimeStep)  +
		sum(m[:dvStorageToGrid][u,ts] for u in p.StorageSalesTiers, ts in p.TimeStep)
		<= p.TimeStepScaling * sum(m[:dvGridPurchase][u,ts] for u in p.PricingTier, ts in p.TimeStep)
	)
end

function add_no_grid_export_constraint(m, p)
	for ts in p.TimeStepsWithoutGrid
		for u in p.SalesTiers
			if !(u in p.CurtailmentTiers)
				for t in p.TechsBySalesTier[u]
					fix(m[:dvProductionToGrid][t, u, ts], 0.0, force=true)
				end
			end
		end
	end
end

function add_energy_price_constraints(m, p)
	##Constraint (10a): Usage limits by pricing tier, by month
	@constraint(m, [u in p.PricingTier, mth in p.Month],
		p.TimeStepScaling * sum( m[:dvGridPurchase][u, ts] for ts in p.TimeStepRatchetsMonth[mth] ) <= m[:binEnergyTier][mth, u] * m[:NewMaxUsageInTier][mth,u])
	##Constraint (10b): Ordering of pricing tiers
	@constraint(m, [u in 2:p.FuelBinCount, mth in p.Month],   #Need to fix, update purchase vs. sales pricing tiers
		m[:binEnergyTier][mth, u] - m[:binEnergyTier][mth, u-1] <= 0)
	## Constraint (10c): One tier must be full before any usage in next tier
	@constraint(m, [u in 2:p.FuelBinCount, mth in p.Month],
		m[:binEnergyTier][mth, u] * m[:NewMaxUsageInTier][mth,u-1] - sum( m[:dvGridPurchase][u-1, ts] for ts in p.TimeStepRatchetsMonth[mth] ) <= 0
	)
	m[:TotalEnergyChargesUtil] = @expression(m, p.pwf_e * p.TimeStepScaling *
		sum( p.ElecRate[u,ts] * m[:dvGridPurchase][u,ts] for ts in p.TimeStep, u in p.PricingTier)
	)
end

function add_monthly_demand_charge_constraints(m, p)
	## Constraint (11a): Upper bound on peak electrical power demand by tier, by month, if tier is selected (0 o.w.)
	@constraint(m, [n in p.DemandMonthsBin, mth in p.Month],
		m[:dvPeakDemandEMonth][mth,n] <= m[:NewMaxDemandMonthsInTier][mth,n] * m[:binDemandMonthsTier][mth,n])

	## Constraint (11b): Monthly peak electrical power demand tier ordering
	@constraint(m, [mth in p.Month, n in 2:p.DemandMonthsBinCount],
		m[:binDemandMonthsTier][mth, n] <= m[:binDemandMonthsTier][mth, n-1])

	## Constraint (11c): One monthly peak electrical power demand tier must be full before next one is active
	@constraint(m, [mth in p.Month, n in 2:p.DemandMonthsBinCount],
		m[:binDemandMonthsTier][mth, n] * m[:NewMaxDemandMonthsInTier][mth,n-1] <= m[:dvPeakDemandEMonth][mth, n-1])

	## Constraint (11d): Monthly peak demand is >= demand at each hour in the month`
	@constraint(m, [mth in p.Month, ts in p.TimeStepRatchetsMonth[mth]],
		sum( m[:dvPeakDemandEMonth][mth, n] for n in p.DemandMonthsBin ) >=
		sum( m[:dvGridPurchase][u, ts] for u in p.PricingTier )
	)
	if !isempty(p.DemandRatesMonth)
		m[:DemandFlatCharges] = @expression(m, p.pwf_e * sum( p.DemandRatesMonth[mth,n] * m[:dvPeakDemandEMonth][mth,n] for mth in p.Month, n in p.DemandMonthsBin) )
	else
		m[:DemandFlatCharges] = 0
	end

end

function add_tou_demand_charge_constraints(m, p)
	## Constraint (12a): Upper bound on peak electrical power demand by tier, by ratchet, if tier is selected (0 o.w.)
	@constraint(m, [r in p.Ratchets, e in p.DemandBin],
		m[:dvPeakDemandE][r, e] <= m[:NewMaxDemandInTier][r,e] * m[:binDemandTier][r, e])

	## Constraint (12b): Ratchet peak electrical power ratchet tier ordering
	@constraint(m, [r in p.Ratchets, e in 2:p.DemandBinCount],
		m[:binDemandTier][r, e] <= m[:binDemandTier][r, e-1])

	## Constraint (12c): One ratchet peak electrical power demand tier must be full before next one is active
	@constraint(m, [r in p.Ratchets, e in 2:p.DemandBinCount],
		m[:binDemandTier][r, e] * m[:NewMaxDemandInTier][r,e-1] <= m[:dvPeakDemandE][r, e-1])

	## Constraint (12d): Ratchet peak demand is >= demand at each hour in the ratchet`
	@constraint(m, [r in p.Ratchets, ts in p.TimeStepRatchets[r]],
		sum( m[:dvPeakDemandE][r, e] for e in p.DemandBin ) >=
		sum( m[:dvGridPurchase][u, ts] for u in p.PricingTier )
	)

	if p.DemandLookbackRange != 0  # then the dvPeakDemandELookback varies by month

		##Constraint (12e): dvPeakDemandELookback is the highest peak demand in DemandLookbackMonths
		for mth in p.Month
			if mth > p.DemandLookbackRange
				@constraint(m, [lm in 1:p.DemandLookbackRange, ts in p.TimeStepRatchetsMonth[mth - lm]],
					m[:dvPeakDemandELookback][mth]
					≥ sum( m[:dvGridPurchase][u, ts] for u in p.PricingTier )
				)
			else  # need to handle rollover months
				for lm in 1:p.DemandLookbackRange
					lkbkmonth = mth - lm
					if lkbkmonth ≤ 0
						lkbkmonth += 12
					end
					@constraint(m, [ts in p.TimeStepRatchetsMonth[lkbkmonth]],
						m[:dvPeakDemandELookback][mth]
						≥ sum( m[:dvGridPurchase][u, ts] for u in p.PricingTier )
					)
				end
			end
		end

		##Constraint (12f): Ratchet peak demand charge is bounded below by lookback
		@constraint(m, [mth in p.Month],
			sum( m[:dvPeakDemandEMonth][mth, n] for n in p.DemandMonthsBin ) >=
			p.DemandLookbackPercent * m[:dvPeakDemandELookback][mth]
		)

	else  # dvPeakDemandELookback does not vary by month

		##Constraint (12e): dvPeakDemandELookback is the highest peak demand in DemandLookbackMonths
		@constraint(m, [lm in p.DemandLookbackMonths],
			m[:dvPeakDemandELookback][1] >= sum(m[:dvPeakDemandEMonth][lm, n] for n in p.DemandMonthsBin)
		)

		##Constraint (12f): Ratchet peak demand charge is bounded below by lookback
		@constraint(m, [mth in p.Month],
			sum( m[:dvPeakDemandEMonth][mth, n] for n in p.DemandMonthsBin ) >=
			p.DemandLookbackPercent * m[:dvPeakDemandELookback][1]
		)
	end

	if !isempty(p.DemandRates)
		m[:DemandTOUCharges] = @expression(m, p.pwf_e * sum( p.DemandRates[r,e] * m[:dvPeakDemandE][r,e] for r in p.Ratchets, e in p.DemandBin) )
	end

end

function add_util_fixed_and_min_charges(m, p)
    m[:TotalFixedCharges] = p.pwf_e * p.FixedMonthlyCharge * 12

	### Constraint (13): Annual minimum charge adder
	if p.AnnualMinCharge > 12 * p.MonthlyMinCharge
        m[:TotalMinCharge] = p.AnnualMinCharge
    else
        m[:TotalMinCharge] = 12 * p.MonthlyMinCharge
    end

	if m[:TotalMinCharge] >= 1e-2
        @constraint(m, MinChargeAddCon, m[:MinChargeAdder] >= m[:TotalMinCharge] - (
			m[:TotalEnergyChargesUtil] + m[:TotalDemandCharges] + m[:TotalExportBenefit] + m[:TotalFixedCharges])
		)
	else
		@constraint(m, MinChargeAddCon, m[:MinChargeAdder] == 0)
	end
end

function add_cost_function(m, p)
	m[:REcosts] = @expression(m,
		# Capital Costs
		m[:TotalTechCapCosts] + m[:TotalStorageCapCosts] +

		## Fixed O&M, tax deductible for owner
		m[:TotalPerUnitSizeOMCosts] * m[:r_tax_fraction_owner] +

        ## Variable O&M, tax deductible for owner
		m[:TotalPerUnitProdOMCosts] * m[:r_tax_fraction_owner] +

		# Utility Bill, tax deductible for offtaker
		(m[:TotalEnergyChargesUtil] + m[:TotalDemandCharges] + m[:TotalExportBenefit] + m[:TotalFixedCharges] + 0.999*m[:MinChargeAdder]) * m[:r_tax_fraction_offtaker] +

        ## Total Generator Fuel Costs, tax deductible for offtaker
        m[:TotalGenFuelCharges] * m[:r_tax_fraction_offtaker] -

        # Subtract Incentives, which are taxable
		m[:TotalProductionIncentive] * m[:r_tax_fraction_owner]
	)
    #= Note: 0.9999*m[:MinChargeAdder] in Obj b/c when m[:TotalMinCharge] > (TotalEnergyCharges + m[:TotalDemandCharges] + TotalExportBenefit + m[:TotalFixedCharges])
		it is arbitrary where the min charge ends up (eg. could be in m[:TotalDemandCharges] or m[:MinChargeAdder]).
		0.0001*m[:MinChargeAdder] is added back into LCC when writing to results.  =#
end

function add_yearone_expressions(m, p)
    m[:Year1UtilityEnergy] = @expression(m,  p.TimeStepScaling * sum(
		m[:dvGridPurchase][u,ts] for ts in p.TimeStep, u in p.PricingTier)
	)
    m[:Year1EnergyCost] = m[:TotalEnergyChargesUtil] / p.pwf_e
    m[:Year1DemandCost] = m[:TotalDemandCharges] / p.pwf_e
    m[:Year1DemandTOUCost] = m[:DemandTOUCharges] / p.pwf_e
    m[:Year1DemandFlatCost] = m[:DemandFlatCharges] / p.pwf_e
    m[:Year1FixedCharges] = m[:TotalFixedCharges] / p.pwf_e
    m[:Year1MinCharges] = m[:MinChargeAdder] / p.pwf_e
    m[:Year1Bill] = m[:Year1EnergyCost] + m[:Year1DemandCost] + m[:Year1FixedCharges] + m[:Year1MinCharges]
end

function reopt(reo_model, model_inputs::Dict)

	t_start = time()
    p = Parameter(model_inputs)
	t = time() - t_start

	results = reopt_run(reo_model, p)
	results["julia_input_construction_seconds"] = t
	return results
end

function reopt_run(m, p::Parameter)

	t_start = time()
	results = Dict{String, Any}()
    Obj = 1  # 1 for minimize LCC, 2 for min LCC AND high mean SOC

	## Big-M adjustments; these need not be replaced in the parameter object.
	add_bigM_adjustments(m, p)
	results["julia_reopt_preamble_seconds"] = time() - t_start
	t_start = time()

	add_continuous_variables(m, p)
	add_integer_variables(m, p)

	results["julia_reopt_variables_seconds"] = time() - t_start
	t_start = time()
    ##############################################################################
	#############  		Constraints									 #############
	##############################################################################

	## Temporary workaround for outages TimeStepsWithoutGrid
	if !isempty(p.TimeStepsWithoutGrid)
		add_no_grid_constraints(m, p)
		add_no_grid_export_constraint(m, p)
	end

	#don't allow curtailment or sales of stroage
	for ts in p.TimeStep
		for u in p.StorageSalesTiers
			fix(m[:dvStorageToGrid][u,ts], 0.0, force=true)
		end
	end

	add_fuel_constraints(m, p)
	add_binTechIsOnInTS_constraints(m, p)

    ### Section 4: Storage System Constraints
	add_storage_size_constraints(m, p)
	add_storage_op_constraints(m, p)
	### Constraint set (5) - hot and cold thermal loads - reserved for later

	### Constraint set (6): Production Incentive Cap
	add_prod_incent_constraints(m, p)
    ### System Size and Production Constraints
	### Constraint set (7): System Size is zero unless single basic tech is selected for class
	if !isempty(p.Tech)
		add_tech_size_constraints(m, p)
	end

	### Constraint set (8): Electrical Load Balancing and Grid Sales
	##Constraint (8a): Electrical Load Balancing with Grid
	add_load_balance_constraints(m, p)

	add_storage_grid_constraints(m, p)
	if !isempty(p.SalesTiers)
		add_prod_grid_constraints(m, p)
	end
	## End constraint (8)

    ### Constraint set (9): Net Meter Module (copied directly from legacy model)
	if !isempty(p.Tech)
		add_NMIL_constraints(m, p)
	end
	##Constraint (9c): Net metering only -- can't sell more than you purchase
	if !isempty(p.SalesTiers)
		add_nem_constraint(m, p)
	end
	###End constraint set (9)

	### Constraint set (10): Electrical Energy Pricing tiers
	add_energy_price_constraints(m, p)

	### Constraint set (11): Peak Electrical Power Demand Charges: binDemandMonthsTier
	add_monthly_demand_charge_constraints(m, p)
	### Constraint set (12): Peak Electrical Power Demand Charges: Ratchets
	if !isempty(p.TimeStepRatchets)
		add_tou_demand_charge_constraints(m, p)
	else
		m[:DemandTOUCharges] = 0
	end
    m[:TotalDemandCharges] = @expression(m, m[:DemandTOUCharges] + m[:DemandFlatCharges])

	add_parameters(m, p)
	add_cost_expressions(m, p)
	add_export_expressions(m, p)
	add_util_fixed_and_min_charges(m, p)

	add_cost_function(m, p)

    if Obj == 1
		@objective(m, Min, m[:REcosts])
	elseif Obj == 2  # Keep SOC high
		@objective(m, Min, m[:REcosts] - sum(m[:dvStorageSOC]["Elec",ts] for ts in p.TimeStep)/8760.)
	end

	results["julia_reopt_constriants_seconds"] = time() - t_start
	t_start = time()

	optimize!(m)

	results["julia_reopt_optimize_seconds"] = time() - t_start
	t_start = time()

	if termination_status(m) == MOI.TIME_LIMIT
		results["status"] = "timed-out"
    elseif termination_status(m) == MOI.OPTIMAL
        results["status"] = "optimal"
    else
		results["status"] = "not optimal"
    end

	##############################################################################
    #############  		Outputs    									 #############
    ##############################################################################
	try
		results["lcc"] = round(JuMP.objective_value(m)+ 0.0001*value(m[:MinChargeAdder]))
	catch
		# not optimal, empty objective_value
		return results
	end

	add_yearone_expressions(m, p)

	results = reopt_results(m, p, results)
	results["julia_reopt_postprocess_seconds"] = time() - t_start
	return results
end

function reopt_results(m, p, r::Dict)
	add_storage_results(m, p, r)
	add_pv_results(m, p, r)
	if !isempty(m[:GeneratorTechs])
		add_generator_results(m, p, r)
    else
		r["gen_net_fixed_om_costs"] = 0
		r["gen_net_variable_om_costs"] = 0
		r["gen_total_fuel_cost"] = 0
		r["gen_year_one_fuel_cost"] = 0
		r["gen_year_one_variable_om_costs"] = 0
    	r["GENERATORtoBatt"] = []
		r["GENERATORtoGrid"] = []
		r["GENERATORtoLoad"] = []
	end
	if !isempty(m[:WindTechs])
		add_wind_results(m, p, r)
	else
		r["WINDtoLoad"] = []
    	r["WINDtoGrid"] = []
	end
	add_util_results(m, p, r)
	return r
end

function add_storage_results(m, p, r::Dict)
    r["batt_kwh"] = value(m[:dvStorageCapEnergy]["Elec"])
    r["batt_kw"] = value(m[:dvStorageCapPower]["Elec"])

    if r["batt_kwh"] != 0
    	@expression(m, soc[ts in p.TimeStep], m[:dvStorageSOC]["Elec",ts] / r["batt_kwh"])
        r["year_one_soc_series_pct"] = value.(soc)
    else
        r["year_one_soc_series_pct"] = []
    end
    @expression(m, GridToBatt[ts in p.TimeStep], m[:dvGridToStorage][ts])
	r["GridToBatt"] = round.(value.(GridToBatt), digits=3)

	@expression(m, ElecFromBatt[ts in p.TimeStep],
		sum(m[:dvDischargeFromStorage][b,ts] for b in p.ElecStorage))
	r["ElecFromBatt"] = round.(value.(ElecFromBatt), digits=3)
	@expression(m, ElecFromBattExport[ts in p.TimeStep],
		sum(m[:dvStorageToGrid][u,ts] for u in p.StorageSalesTiers))
	r["ElecFromBattExport"] = round.(value.(ElecFromBattExport), digits=3)
	nothing
end

function add_generator_results(m, p, r::Dict)
	m[:GenPerUnitSizeOMCosts] = @expression(m, p.two_party_factor *
		sum(p.OMperUnitSize[t] * p.pwf_om * m[:dvSize][t] for t in m[:GeneratorTechs])
	)
	m[:GenPerUnitProdOMCosts] = @expression(m, p.two_party_factor *
		sum(m[:dvRatedProduction][t,ts] * p.TimeStepScaling * p.ProductionFactor[t,ts] * p.OMcostPerUnitProd[t] * p.pwf_om
			for t in m[:GeneratorTechs], ts in p.TimeStep)
	)
	if value(sum(m[:dvSize][t] for t in m[:GeneratorTechs])) > 0
		r["generator_kw"] = value(sum(m[:dvSize][t] for t in m[:GeneratorTechs]))
		r["gen_net_fixed_om_costs"] = round(value(m[:GenPerUnitSizeOMCosts]) * m[:r_tax_fraction_owner], digits=0)
		r["gen_net_variable_om_costs"] = round(value(m[:GenPerUnitProdOMCosts]) * m[:r_tax_fraction_owner], digits=0)
		r["gen_total_fuel_cost"] = round(value(m[:TotalGenFuelCharges]) * m[:r_tax_fraction_offtaker], digits=2)
		r["gen_year_one_fuel_cost"] = round(value(m[:TotalGenFuelCharges]) / p.pwf_e, digits=2)
		r["gen_year_one_variable_om_costs"] = round(value(m[:GenPerUnitProdOMCosts]) / (p.pwf_om * p.two_party_factor), digits=0)
		r["gen_year_one_fixed_om_costs"] = round(value(m[:GenPerUnitSizeOMCosts]) / (p.pwf_om * p.two_party_factor), digits=0)
	end
	@expression(m, GENERATORtoBatt[ts in p.TimeStep],
				sum(m[:dvProductionToStorage]["Elec",t,ts] for t in m[:GeneratorTechs]))
	r["GENERATORtoBatt"] = round.(value.(GENERATORtoBatt), digits=3)

	@expression(m, GENERATORtoGrid[ts in p.TimeStep],
				sum(m[:dvProductionToGrid][t,u,ts] for t in m[:GeneratorTechs], u in p.SalesTiersByTech[t]))
	r["GENERATORtoGrid"] = round.(value.(GENERATORtoGrid), digits=3)

	@expression(m, GENERATORtoLoad[ts in p.TimeStep],
				sum(m[:dvRatedProduction][t, ts] * p.ProductionFactor[t, ts] * p.LevelizationFactor[t]
					for t in m[:GeneratorTechs]) -
					GENERATORtoBatt[ts] - GENERATORtoGrid[ts]
					)
	r["GENERATORtoLoad"] = round.(value.(GENERATORtoLoad), digits=3)

    @expression(m, GeneratorFuelUsed, sum(m[:dvFuelUsage][t, ts] for t in m[:GeneratorTechs], ts in p.TimeStep))
	r["fuel_used_gal"] = round(value(GeneratorFuelUsed), digits=2)

	m[:Year1GenProd] = @expression(m,
		p.TimeStepScaling * sum(m[:dvRatedProduction][t,ts] * p.ProductionFactor[t, ts]
			for t in m[:GeneratorTechs], ts in p.TimeStep)
	)
	r["year_one_gen_energy_produced"] = round(value(m[:Year1GenProd]), digits=0)
	m[:AverageGenProd] = @expression(m,
		p.TimeStepScaling * sum(m[:dvRatedProduction][t,ts] * p.ProductionFactor[t, ts] * p.LevelizationFactor[t]
			for t in m[:GeneratorTechs], ts in p.TimeStep)
	)
	r["average_yearly_gen_energy_produced"] = round(value(m[:AverageGenProd]), digits=0)

	nothing
end

function add_wind_results(m, p, r::Dict)
	r["wind_kw"] = round(value(sum(m[:dvSize][t] for t in m[:WindTechs])), digits=4)
	@expression(m, WINDtoBatt[ts in p.TimeStep],
	            sum(sum(m[:dvProductionToStorage][b, t, ts] for t in m[:WindTechs]) for b in p.ElecStorage))
	r["WINDtoBatt"] = round.(value.(WINDtoBatt), digits=3)

	@expression(m, WINDtoCurtail[ts in p.TimeStep],
				sum(m[:dvProductionToGrid][t,u,ts] for t in m[:WindTechs], u in p.CurtailmentTiers))

	r["WINDtoCurtail"] = round.(value.(WINDtoCurtail), digits=3)

	@expression(m, WINDtoGrid[ts in p.TimeStep],
				sum(m[:dvProductionToGrid][t,u,ts] for t in m[:WindTechs], u in p.SalesTiersByTech[t]) - WINDtoCurtail[ts])

	r["WINDtoGrid"] = round.(value.(WINDtoGrid), digits=3)

	@expression(m, WINDtoLoad[ts in p.TimeStep],
				sum(m[:dvRatedProduction][t, ts] * p.ProductionFactor[t, ts] * p.LevelizationFactor[t]
					for t in m[:WindTechs]) - WINDtoGrid[ts] - WINDtoBatt[ts] - WINDtoCurtail[ts] )
	r["WINDtoLoad"] = round.(value.(WINDtoLoad), digits=3)
	m[:Year1WindProd] = @expression(m,
		p.TimeStepScaling * sum(m[:dvRatedProduction][t,ts] * p.ProductionFactor[t, ts]
			for t in m[:WindTechs], ts in p.TimeStep)
	)
	r["year_one_wind_energy_produced"] = round(value(m[:Year1WindProd]), digits=0)
	m[:AverageWindProd] = @expression(m,
		p.TimeStepScaling * sum(m[:dvRatedProduction][t,ts] * p.ProductionFactor[t, ts] * p.LevelizationFactor[t]
			for t in m[:WindTechs], ts in p.TimeStep)
	)
	r["average_wind_energy_produced"] = round(value(m[:AverageWindProd]), digits=0)
	nothing
end

function add_pv_results(m, p, r::Dict)
	PVclasses = filter(tc->startswith(tc, "PV"), p.TechClass)
    for PVclass in PVclasses
		PVtechs_in_class = filter(t->startswith(t, PVclass), m[:PVTechs])

		if !isempty(PVtechs_in_class)

			r[string(PVclass, "_kw")] = round(value(sum(m[:dvSize][t] for t in PVtechs_in_class)), digits=4)

			# NOTE: must use anonymous expressions in this loop to overwrite values for cases with multiple PV
            if !isempty(p.ElecStorage)
				PVtoBatt = @expression(m, [ts in p.TimeStep],
					sum(m[:dvProductionToStorage][b, t, ts] for t in PVtechs_in_class, b in p.ElecStorage))
			else
				PVtoBatt = @expression(m, [ts in p.TimeStep], 0.0)
            end
			r[string(PVclass, "toBatt")] = round.(value.(PVtoBatt), digits=3)

			PVtoCurtail = @expression(m, [ts in p.TimeStep],
					sum(m[:dvProductionToGrid][t,u,ts] for t in PVtechs_in_class, u in p.CurtailmentTiers))
    	    r[string(PVclass, "toCurtail")] = round.(value.(PVtoCurtail), digits=3)

			PVtoGrid = @expression(m, [ts in p.TimeStep],
					sum(m[:dvProductionToGrid][t,u,ts] for t in PVtechs_in_class, u in p.SalesTiersByTech[t]) - PVtoCurtail[ts])
    	    r[string(PVclass, "toGrid")] = round.(value.(PVtoGrid), digits=3)

			PVtoLoad = @expression(m, [ts in p.TimeStep],
				sum(m[:dvRatedProduction][t, ts] * p.ProductionFactor[t, ts] * p.LevelizationFactor[t] for t in PVtechs_in_class)
				- PVtoGrid[ts] - PVtoBatt[ts] - PVtoCurtail[ts]
				)
            r[string(PVclass, "toLoad")] = round.(value.(PVtoLoad), digits=3)

			Year1PvProd = @expression(m, sum(m[:dvRatedProduction][t,ts] * p.ProductionFactor[t, ts]
				for t in PVtechs_in_class, ts in p.TimeStep) * p.TimeStepScaling)
			r[string("year_one_energy_produced_", PVclass)] = round(value(Year1PvProd), digits=0)

			AveragePvProd = @expression(m, sum(m[:dvRatedProduction][t,ts] * p.ProductionFactor[t, ts] * p.LevelizationFactor[t]
			    for t in PVtechs_in_class, ts in p.TimeStep) * p.TimeStepScaling)
            r[string("average_yearly_energy_produced_", PVclass)] = round(value(AveragePvProd), digits=0)

			ExportedElecPV = @expression(m, sum(m[:dvProductionToGrid][t,u,ts]
				for t in PVtechs_in_class, u in p.SalesTiersByTech[t], ts in p.TimeStep) * p.TimeStepScaling)
            r[string("average_annual_energy_exported_", PVclass)] = round(value(ExportedElecPV), digits=0)

            PVPerUnitSizeOMCosts = @expression(m, sum(p.OMperUnitSize[t] * p.pwf_om * m[:dvSize][t] for t in PVtechs_in_class))
            r[string(PVclass, "_net_fixed_om_costs")] = round(value(PVPerUnitSizeOMCosts) * m[:r_tax_fraction_owner], digits=0)
        end
	end
end

function add_util_results(m, p, r::Dict)
    net_capital_costs_plus_om = value(m[:TotalTechCapCosts] + m[:TotalStorageCapCosts]) +
                                value(m[:TotalPerUnitSizeOMCosts] + m[:TotalPerUnitProdOMCosts]) * m[:r_tax_fraction_owner] +
                                value(m[:TotalGenFuelCharges]) * m[:r_tax_fraction_offtaker]

    push!(r, Dict("year_one_utility_kwh" => round(value(m[:Year1UtilityEnergy]), digits=2),
						 "year_one_energy_cost" => round(value(m[:Year1EnergyCost]), digits=2),
						 "year_one_demand_cost" => round(value(m[:Year1DemandCost]), digits=2),
						 "year_one_demand_tou_cost" => round(value(m[:Year1DemandTOUCost]), digits=2),
						 "year_one_demand_flat_cost" => round(value(m[:Year1DemandFlatCost]), digits=2),
						 "year_one_export_benefit" => round(value(m[:ExportBenefitYr1]), digits=0),
						 "year_one_fixed_cost" => round(m[:Year1FixedCharges], digits=0),
						 "year_one_min_charge_adder" => round(value(m[:Year1MinCharges]), digits=2),
						 "year_one_bill" => round(value(m[:Year1Bill]), digits=2),
						 "year_one_payments_to_third_party_owner" => round(value(m[:TotalDemandCharges]) / p.pwf_e, digits=0),
						 "total_energy_cost" => round(value(m[:TotalEnergyChargesUtil]) * m[:r_tax_fraction_offtaker], digits=2),
						 "total_demand_cost" => round(value(m[:TotalDemandCharges]) * m[:r_tax_fraction_offtaker], digits=2),
						 "total_fixed_cost" => round(m[:TotalFixedCharges] * m[:r_tax_fraction_offtaker], digits=2),
						 "total_export_benefit" => round(value(m[:TotalExportBenefit]) * m[:r_tax_fraction_offtaker], digits=2),
						 "total_min_charge_adder" => round(value(m[:MinChargeAdder]) * m[:r_tax_fraction_offtaker], digits=2),
						 "net_capital_costs_plus_om" => round(net_capital_costs_plus_om, digits=0),
						 "average_annual_energy_exported_wind" => round(value(m[:ExportedElecWIND]), digits=0),
						 "average_annual_energy_curtailed_wind" => round(value(m[:CurtailedElecWIND]), digits=0),
                         "average_annual_energy_exported_gen" => round(value(m[:ExportedElecGEN]), digits=0),
						 "net_capital_costs" => round(value(m[:TotalTechCapCosts] + m[:TotalStorageCapCosts]), digits=2))...)

    @expression(m, GridToLoad[ts in p.TimeStep],
                sum(m[:dvGridPurchase][u,ts] for u in p.PricingTier) - m[:dvGridToStorage][ts] )
    r["GridToLoad"] = round.(value.(GridToLoad), digits=3)
end


text:
using Plots

function f(r)

    x = sqrt(r/M)
    x0 = sqrt(R_isco/M)
    x1 = 2*cos((1/3)*(acos(a_star))-(pi/3))
    x2 = 2*cos((1/3)*(acos(a_star))+(pi/3))
    x3 = -2*cos((1/3)*(acos(a_star)))

    flux = (3/(2*M))*(1/(x^(2)*(x^(3)-(3*x)+(2*a_star))))*(x-x0-((3/2)*a_star*log(x/x0))
            - (3*(x1-a_star)^2)/(x1*(x1-x2)*(x1-x3))*log((x-x1)/(x0-x1))
            - (3*(x2-a_star)^2)/(x2*(x2-x1)*(x2-x3))*log((x-x2)/(x0-x2))
            - (3*(x3-a_star)^2)/(x3*(x3-x1)*(x3-x2))*log((x-x3)/(x0-x3)))

    return flux

end

function D(R, f)
    D = mdot*f/(4*pi*R)
    return D
end

function (T(D))
    T = (D/sigma_SB)^(1/4)
    return T
end

function r_isco(a_star)

    a = a_star
    r_g = 2*G*M/(c^(2))
    z1 = 1+(1-a^2)^(1/3)*((1+a)^(1/3)+(1-a)^(1/3))
    z2 = (3*a^(2)+z1^(2))^(1/2)
    if a >= 0
        r_isco = (3+z2-(3-z1)*(3+z1+2*z2)^(1/2))*r_g
    elseif a < 0
        r_isco = (3+z2+(3-z1)*(3+z1+2*z2)^(1/2))*r_g
    end

    return r_isco
end

G = 6.67e-11
c = 3e8
L_sun = 3.8e26
M_sun = 1.99e30
sigma_SB = 5.67e-8

eta = 0.1
M = 10e6*M_sun
a_star = 0.5
R_isco = r_isco(a_star)
L_edd = 3e4*L_sun*(M/M_sun)
Mdot = -L_edd/(c^2*eta)
mdot = 0.1*Mdot

r_vals = LinRange(R_isco, 10*R_isco, 10000)
D_vals = []
T_vals = []

for r in r_vals
    push!(D_vals, D(r, f(r)))
    push!(T_vals, T(D(r, f(r))))
end

plot(r_vals, D_vals)
plot(r_vals, T_vals)

max_stars_repo_path: src/Rel_PT_Emission Profile.jl

max_issues_repo_path: src/Rel_PT_Emission Profile.jl

max_forks_repo_path: src/Rel_PT_Emission Profile.jl


text:
"""

# `dag_vars`

$(SIGNATURES)

Part of the API, exported
"""
function dag_vars(d::OrderedDict)
  vars = Symbol[]
  for var in keys(d)
    if isnothing(var)
      @warn "LHS can't be an empty set: $var."
    elseif typeof(var) == Symbol
      append!(vars, [var])
      handle_rhs!(vars, d[var])
    elseif typeof(var) == Vector{Symbol}
      append!(vars, var)
      handle_rhs!(vars, d[var])
    end
  end
  unique(vars)
end

function handle_rhs!(vars::Vector{Symbol}, rhs::SymbolList)
  if isnothing(var)
    append!(vars, [])
  elseif typeof(rhs) == Symbol
    append!(vars, [rhs])
  elseif typeof(rhs) == Vector{Symbol}
    append!(vars, rhs)
  end
end

"""

# `edge_matrix`

$(SIGNATURES)

Part of the API, exported
"""
function edge_matrix(d::OrderedDict)
  vars = dag_vars(d)
  l = length(vars)
  a = zeros(Int, l, l)

  for key in keys(d)
    if typeof(key) == Symbol
      for (ind, var) in enumerate(vars)
        # Is this Symbol the LHS
        if var == key
          # What's in the RHS
          for (i, vr) in enumerate(vars)
            if typeof(d[key]) == Symbol
              a[i, ind] += vr == d[key] ? 1 : 0
            elseif typeof(d[key]) == Vector{Symbol}
              if vr in d[key]
                a[i, ind] += 1
              end
            end
          end
        end
      end
    elseif typeof(key) == Vector{Symbol}
      for (ind, var) in enumerate(vars)
        # Are these Symbols (key could hold a single Symbol) the LHS
        if var in key
          for (i, vr) in enumerate(vars)
            if typeof(d[key]) == Symbol
              a[i, ind] += vr == d[key] ? 1 : 0
            elseif typeof(d[key]) == Vector{Symbol}
              if vr in d[key]
                a[i, ind] += 1
              end
            end
          end
        end
      end
    end
  end
  NamedArray(Int.(a), (vars, vars), ("Rows", "Cols"))
end

"""

# `edge_matrix`

$(SIGNATURES)

Part of the API, exported
"""
function edge_matrix(a::NamedArray, inv=false)
  a = sign.(a)
  if inv
    ord = topological_order(a)
    ord = reverse(ord)
    a = a[ord, ord]
  end
  transpose(a) + I(size(a, 1))
end

"""

# `adjacency_matrix`

$(SIGNATURES)

Part of the API, exported
"""
function adjacency_matrix(d::OrderedDict)
  transpose(edge_matrix(d))
end

"""

# `adjacency_matrix`

$(SIGNATURES)

Part of the API, exported
"""
function adjacency_matrix(e::NamedArray)
  a = transpose(e)
  a - I(size(a, 1))
end

"""

# `adjacency_matrix_to_dict`

$(SIGNATURES)

Part of the API, exported.
"""
function adjacency_matrix_to_dict(a::NamedArray)
  vars = names(a, 1)
  dct = OrderedDict()
  for (ind, r) in enumerate(eachrow(a))
    rhs = vars[findall(x -> x ==1, r)]
    if length(rhs) == 1
      dct[vars[ind]] = vars[findall(x -> x ==1, r)][1]
    elseif length(rhs) > 1
      dct[vars[ind]] = vars[findall(x -> x ==1, r)]
    end
  end
  dct
end

"""

# `undirected_matrix`

$(SIGNATURES)

Internal
"""
function undirected_matrix(d::DAG)
  u = Array(d.a) + Array(d.e)
  vars = names(d.a)[1]
  NamedArray(Int.(u), (vars, vars), ("Rows", "Cols"))
end

"""

# `topological_order`

$(SIGNATURES)

Part of the API, exported
"""
function topological_order(a::NamedArray)
  #@assert is_dag(a)
  n = size(a, 1)
  nod = 1:n
  indeg = zeros(Int, n)
  if sum(tril(a)) == sum(tril(ones(Int, n, n)))
    return nod
  end
  zero_indeg = Int[]
  for i in nod
    indeg[i] = sum(a[:, i])
    if indeg[i] == 0
      append!(zero_indeg, [i])
    end
  end
  s = 1
  ord = zeros(Int, n)
  while length(zero_indeg) > 0
    v = pop!(zero_indeg)
    ord[s] = v
    s += 1
    cs = filter(i -> a[v, i] !== 0, nod)
    length(cs) == 0 && continue
    for j in 1:length(cs)
      k = cs[j]
      indeg[k] = indeg[k] - 1
      indeg[k] == 0 && push!(zero_indeg, k)
    end
  end
  reverse(ord)
end

"""

# `topological_sort`

$(SIGNATURES)

Part of the API, exported
"""
function topological_sort(a::NamedArray)
  ord = topological_order(a)
  a[ord, ord]
end

"""

# `topological_sort`

$(SIGNATURES)

Part of the API, exported
"""
function topological_sort(dag::DAG)
  new_dag = deepcopy(dag)
  ord = topological_order(new_dag.a)
  new_dag.a = new_dag.a[ord, ord]
  new_dag.e = new_dag.a[ord, ord]
  new_dag.vars = new_dag.vars[ord]
  new_dag
end

"""

# `topological_sort!`

$(SIGNATURES)

Part of the API, exported
"""
function topological_sort!(dag::DAG)
  ord = topological_order(dag.a)
  dag.a = dag.a[ord, ord]
  dag.e = dag.a[ord, ord]
  dag.vars = dag.vars[ord]
end

export
  dag_vars,
  adjacency_matrix,
  edge_matrix,
  topological_sort,
  topological_order,
  adjacency_matrix_to_dict


text:
# """
#    abstract type Source has two subtype: SingleSource and SimultaneousSource
# """
# abstract type Source end

"""
Defines an explosive source for time domain finite difference modelling.

# Fields

- `isz::Ti`: vertical index location
- `isx::Ti`: horizontal index location
- `src2spt::Ti` : index mapping source location to snapshot
- `src2wfd::Ti` : index mapping source location to wavefield
- `dt::Tv`      : time sampling
- `it_min::Ti`  : source lower time index
- `it_max::Ti`  : source upper time index
- `p::Vector{Tv}` :  Vector of Source wavelets
- `hello`
"""
struct Source{Ti<:Int64, Tv<:AbstractFloat}
    isz     :: Ti  # vertical index location
    isx     :: Ti  # horizontal index location
    src2spt :: Ti  # index mapping to snapshot
    src2wfd :: Ti  # index mapping to wavefield
    dt      :: Tv
    it_min  :: Ti  # the index lower bound in time axis
    it_max  :: Ti  # the index upper bound in time axis
    p       :: Vector{Tv} # source wavelet
end

"""
   create a ricker wavelet
"""
function ricker(fdom, dt)

 	  nw = 2.2 / fdom / dt         # decide number of samples
	  nw = 2*floor(Int64, nw/2)+1  # gaurantee the number of samples is odd

	  nc = floor(Int64, nw/2)
	  w  = zeros(nw)               # allocate space for wavelet
	  k  = collect(1:nw)

    for i = 1 : nw
        alpha = (nc-i+1) * fdom * dt * pi
        beta  = alpha^2
        w[i]  = (1 - beta * 2) * exp(-beta)
    end
	  return w
end

"""
  convert wavelet to its corresponding minimum phase version
"""
function convert2miniphase(w::Vector{Tv}) where {Tv<:AbstractFloat}

    nw = length(w)

    # number of samples
    nf = 8*nextpow(2, nw)

    W = fft(vcat(w,zeros(nf-nw)))
    A = log.(abs.(W) .+ 1.e-8)
    a = 2*ifft(A)
    n2 = floor(Int, nf/2)
    a[n2+2:nf] .= 0.0
    a[1] = a[1]/2
    A = exp.(fft(a))
    a = real(ifft(A))
    wmin = real(a[1:nw])

    return wmin
end

"""
   fc is the high cut frequency of low-pass filter, hl is half length of filter
and dt is sampling rate
"""
function tapered_sinc(fc, hl, dt) where {Tv <: AbstractFloat}

    # length of filter
    l = 2*hl
    filter = zeros(l+1)
    ff = fc / (1.0/dt)

    for i = 0 : l
        idx = i - hl
        if idx == 0
           filter[i+1] = 2*pi*ff
        else
           filter[i+1] = sin(2*pi*ff*idx)/idx
        end
        taper = 0.42 - 0.5 * cos(2*pi*i/l) + 0.08*cos(4*pi*i/l)
        filter[i+1] = filter[i+1] * taper
    end

    # normalize filter
    filter = filter / sum(filter)
    return filter
end

"""
   scaled sinc function without tapering, so the decaying of sinc function is slow
"""
function scaled_sinc(fc, hl, dt::Tv) where {Tv<:AbstractFloat}

    l = 2*hl
    filter = zeros(l+1)
    # high cut frequency fraction
    ff = fc / (1.0/dt)

    for i = 0 : l
        idx = i - hl
        if idx == 0
           filter[i+1] = 2*pi*ff
        else
           filter[i+1] = sin(2*pi*ff*idx)/idx
        end
    end

    filter = filter / sum(filter)
    return filter
end

"""
   S=Source(sz, sx, params; <keyword arguments>)

Constructs an explosive Source structure.

Please see the extended help with `?? Source`

# Arguments

-`sz::`: source coordinate
-`sx::`: soirce cordinate
-`params::TdParams`: Parameters for time domain modelling

# Keyword Arguments

-`location_flag`: Source(s) positioning by `"index"` or `"distance"`
-`ot`: Source(s) firing time
-`amp`: Source(s) amplitude
-`fdom`: Dominant frequency of wavelet
-`type_flag`: Type of source wavelet.
-`hl`:Window lenght for sinc interpolation (for type_flag="sinc")
-`p`: Either a user provided source wavelet or a vector of source wavelets.

# Extended help

   S=Source(sz, sx, params; location_flag="index", ot=0.0, amp=1.0, fdom=20.0, hl=128, type_flag="ricker", p::Vector{Float32}=p)

The structure defines the source location by either its position or relative index in the model through the `location_flag`.

If the source wavelet (or a vector of source wavelets) `p` is not provieded, the user can specify it with "type_flag".

-Default Ricker wavelet: `type_flag="ricker"`
-Sinc wavelet: `type_flag="sinc"`
-Minimum phase wavelet: `type_flag=miniphase`
"""
function Source(sz, sx, params::TdParams; location_flag="index", ot=0.0, amp=1.0, fdom=20.0, hl=128,
                type_flag="ricker", p=Vector{Float32}(undef,0)) where {Tv <: Real}

    # source location given index
    if location_flag == "index"
       isz = round(Int64, sz)
       isx = round(Int64, sx)

    # source location given as distance
    elseif location_flag == "distance"
       isz = round(Int64, sz/params.dz) + 1
       isx = round(Int64, sx/params.dx) + 1
    else
       error("wrong specification of source location")
    end

    # error checking
    if isz > params.nz || isx > params.nx || isz < 1 || isx < 1
       error("source located outside of modeling area")
    end

    # can't put on the free-surface
    if params.free_surface && isz == 1
       error("can't inject source on the surface")
    end

    # the auxillary index mapping the location of source to snapshot or wavefield
    src2spt = (isx+params.npml-1) * params.Nz + isz + params.ntop
    src2wfd = (isx            -1) * params.nz + isz

    # non-user provided source wavelet
    if length(p) == 0
       if type_flag == "ricker"
          p = ricker(fdom, params.dt)
          p = amp * p
       elseif type_flag == "miniphase"
          p = ricker(fdom, params.dt)
          p = convert2miniphase(p)
          p = amp * p
       elseif type_flag == "sinc"
          fc = params.data_format(fdom)
          hl = round(Int64, hl)
          p  = tapered_sinc(fc, hl, params.dt)
          p  = amp * p
       end
    end

    # make the data type consistent
    p = convert(Vector{params.data_format}, p)
    nt= length(p)

    # the index time range of source wavelet
    it_min = floor(Int64, ot/params.dt) + 1
    it_max = it_min + nt - 1

    # call the default constructor
    return Source(isz, isx, src2spt, src2wfd, params.dt, it_min, it_max, p)
end

"""
   source_array= get_multi_sources(sz::Vector, sx::Vector, params::TdParams; <keyword arguments>)

Get an array of sources with different positions.

Every element of the source array contains an initialized Source structure.

# Arguments:

-`sz::Vector`: source z coordinates.
-`sx::Vector`: source x coordinates.
-`params::TdParams`: Time domain parameters for modeling.

# Extended help

This function extends the function constructor `Source` to a vector of source.

-`location_flag::String`: Define the type of location for each source.
-`ot::Vector`: Vector of firing times for each source.
-`amp::Float64`: Amplitude for the array of sources.
-`fdom::Float64`: Dominant frequency the array of sources.
-`hl::`:
-`typef_flag::String`: Type of source (default: "ricker")
-`p::Vector{Float32}`: Vector of wavelets for the array of soruces.
"""
function get_multi_sources(sz::Vector, sx::Vector, params::TdParams; location_flag="index", ot=0.0, amp=1.0,
                          fdom=20.0, hl=128, type_flag="ricker", p=Vector{Float32}(undef,0))

    # number of source
    ns = length(sz)
    if length(sx) != ns
       error("length of source coordinate mismatch")
    end

    # source location given index
    isz = zeros(Int64, ns)
    isx = zeros(Int64, ns)
    if location_flag == "index"
       for i = 1 : ns
           isz[i] = round(Int64, sz[i])
           isx[i] = round(Int64, sx[i])
       end

    # source location given as distance
    elseif location_flag == "distance"
       isz[i] = round(Int64, sz[i]/params.dz) + 1
       isx[i] = round(Int64, sx[i]/params.dx) + 1
    else
       error("wrong specification of source location")
    end

    # error checking
    for i = 1 : ns
        if isz[i] > params.nz || isx[i] > params.nx
           error("source located outside of modeling area")
        end
    end

    # starting time
    if typeof(ot) <: Real  # all sources starting at same time
       ot = ot * ones(ns)
    end
    if length(ot) != ns
       error("length of starting time ot mismatch")
    end

    # amplitude for each source
    if typeof(amp) <: Real
       amp = amp * ones(ns)
    end
    if length(amp) != ns
       error("length of amplitude is wrong")
    end

    # dominant frequency for source wavelet
    if typeof(fdom) <: Real
       fdom = fdom * ones(ns)
    end
    if length(fdom) != ns
       error("length of dominant frequency is wrong")
    end

    wavelet = Vector{Vector{params.data_format}}(undef, ns)
    if length(p) == 0  # without user provided wavelet
       if typeof(type_flag) == String               # all sources share same wavelet
          wavelet_type = Vector{String}(undef, ns)
          wavelet_type[:] .= type_flag
          type_flag    = wavelet_type
       elseif typeof(type_flag) == Vector{String}   # use different wavelet
          if length(type_flag) != ns
             error("length of type_flag mismatch")
          end
       end

       # generate the specific source wavelet
       for i = 1 : ns
           if type_flag[i] == "ricker"
              tmp = amp[i] * ricker(fdom[i], params.dt)
              wavelet[i] = convert(Vector{params.data_format}, tmp)
           elseif type_flag[i] == "miniphase"
              tmp = ricker(fdom[i], params.dt)
              tmp = amp[i] * convert2miniphase(tmp)
              wavelet[i] = convert(Vector{params.data_format}, tmp)
           elseif type_flag[i] == "sinc"
              hl   = round(Int64, hl)
              tmp  = amp[i] * tapered_sinc(fdom[i], hl, params.dt)
              wavelet[i] = convert(Vector{params.data_format}, tmp)
           end
       end

    elseif eltype(p) <: AbstractFloat   # all sources share same user-given wavelet
       for i = 1 : ns
           wavelet[i] = convert(Vector{params.data_format}, p)
       end

    else                                # sources have different user-given wavelet
       for i = 1 : ns
           wavelet[i] = convert(Vector{params.data_format}, p[i])
       end
    end

    # allocate memory for vector of sources
    srcs = Vector{Source}(undef, ns)
    for i = 1 : ns
        srcs[i] = Source(isz[i], isx[i], params; p=wavelet[i], ot=ot[i])
    end
    return srcs
end

"""
   Add source to snapshot (we add source to pressure terms, so the wavelet is caled by dt)
"""
function add_source!(spt::Snapshot, src::Source, it::Int64)

    # add source
    if src.it_min <= it <= src.it_max

       idx_t = it - src.it_min + 1
       idx_p = src.src2spt

       # equvialent to adding to pressure field
       spt.px[idx_p] = spt.px[idx_p] + src.p[idx_t] * src.dt
    end

    return nothing
end

"""
   add source to wavefield, used for forward wavefield reconstruction
"""
function add_source!(wfd::Wavefield, src::Source, it::Int64)

    if src.it_min <= it <= src.it_max

       idx_t = it - src.it_min + 1
       idx_p = src.src2wfd

       wfd.p[idx_p] = wfd.p[idx_p] + src.p[idx_t] * src.dt
    end
    return nothing
end

"""
   add multi-sources simultaneously to snapshot
"""
function add_source!(spt::Snapshot, srcs::Vector{Source}, it::Int64)

    # number of sources
    ns = length(srcs)

    # loop over each source
    for i = 1 : ns
        add_source!(spt, srcs[i], it)
    end
    return nothing
end

"""
   add multi-sources simultaneously to snapshot
"""
function add_source!(wfd::Wavefield, srcs::Vector{Source}, it::Int64)

    # number of sources
    ns = length(srcs)

    # loop over each source
    for i = 1 : ns
        add_source!(wfd, srcs[i], it)
    end
    return nothing
end

"""
   subtract source from wavefield, used for backward wavefield reconstruction
"""
function subtract_source!(wfd::Wavefield, src::Source, it::Int64)

    if src.it_min <= it <= src.it_max
       idx_t = it - src.it_min + 1
       idx_p = src.src2wfd

       wfd.p[idx_p] = wfd.p[idx_p] - src.p[idx_t] * src.dt
    end
    return nothing
end

"""
   subtract multiple sources from wavefield, used for backward wavefield reconstruction
"""
function subtract_source!(wfd::Wavefield, srcs::Vector{Source}, it::Int64)

    # number of sources
    ns = length(srcs)

    # loop over sources
    for i = 1 : ns
        if srcs[i].it_min <= it <= srcs[i].it_max
           idx_t = it - srcs[i].it_min + 1
           idx_p = srcs[i].src2wfd

           wfd.p[idx_p] = wfd.p[idx_p] - srcs[i].p[idx_t] * srcs[i].dt
        end
    end

    return nothing
end

"""
   find the time range of multiple sources
"""
function time_range_source(srcs::Vector{Source})

    tmax = 0.0
    tmin = 0.0
    ns   = length(srcs)

    for i = 1 : ns
        tmp1 = (srcs[i].it_min-1) * srcs[i].dt
        tmp2 = (srcs[i].it_max-1) * srcs[i].dt

        if tmp1 < tmin
           tmin = tmp1
        end
        if tmp2 > tmax
           tmax = tmp2
        end
    end

    return tmin, tmax
end

function write_source(path_src::String, src::Source)

    fid = open(path_src, "w")

    write(fid, src.isz)
    write(fid, src.isx)
    write(fid, src.src2spt)
    write(fid, src.src2wfd)
    write(fid, src.it_min)
    write(fid, src.it_max)

    typeof(src.dt) == Float32 && write(fid, 1)
    typeof(src.dt) == Float64 && write(fid, 2)

    write(fid, src.dt)
    write(fid, src.p )

    close(fid)
    return nothing
end

function read_source(path_src::String)

    fid     = open(path_src, "r")
    isz     = read(fid, Int64)
    isx     = read(fid, Int64)
    src2spt = read(fid, Int64)
    src2wfd = read(fid, Int64)
    it_min  = read(fid, Int64)
    it_max  = read(fid, Int64)
    nt      = it_max - it_min + 1
    format_code = read(fid, Int64)

    if format_code == 1
       dt = read(fid, Float32)
       p  = zeros(Float32, nt)
       read!(fid, p)

    elseif format_code == 2
       dt = read(fid, Float64)
       p  = zeros(Float64, nt)
       read!(fid, p)

    else
       error("wrong data format")
    end

    return Source(isz, isx, src2spt, src2wfd, dt, it_min, it_max, p)
end

function source_isequal(src1::Source, src2::Source)

    src1.isz     != src2.isz && (return false)
    src1.isx     != src2.isx && (return false)
    src1.src2spt != src2.src2spt && (return false)
    src1.src2wfd != src2.src2wfd && (return false)
    src1.it_min  != src2.it_min  && (return false)
    src1.it_max  != src2.it_max  && (return false)
    src1.dt      != src2.dt      && (return false)
    norm(src1.p-src2.p)/norm(src1.p) > 10*eps(eltype(src1.p)) && (return false)

    return true
end

# ========================need to be updated ===================================
# """
#    randomly source starting time shift and make the starting time is the integer
# times the time interval.
# """
# function random_start_time(ns::Ti, dt::Tv, s::Tv) where {Ti<:Int64, Tv<:AbstractFloat}
#
#     ot = s * rand(ns)
#     indt = round(Int64, ot/dt)
#     ot[:] = indt*dt
#
#     return convert(Vector{typeof(dt)}, ot)
# end

# """
#    generate a snapshot via multiple sources
# """
# function srcs2spt(srcs::Array{Source,1}, it::Int64)
#
#     nz = srcs[1].nz ; nx = srcs[1].nx;
#     ext= srcs[1].ext; iflag=srcs[1].iflag;
#     dt = srcs[1].dt
#     spt = InitSnapShot(nz, nx, ext, iflag, dt, it)
#
#     AddMultiSources!(spt, srcs)
#     return spt
# end
#
# """
#    generate wavefield via multiple sources and write the snapshots to hard drive
# """
# function WriteSrcs2SnapShots(path::String , srcs::Array{Source,1})
#     nz = srcs[1].nz
#     (tl, tu) = SrcRange(srcs)
#     nt  = round(Int64, (tu-tl)/dt) + 1
#     fid = open(path, "w")
#     write(fid, nz, nx, ext, iflag, dt)
#     for it = 1 : nt
#         it1 = round(Int64, (tl+(it-1)*dt)/dt+1)
#         spt = srcs2spt(srcs, it1)
#         writeSnapShot(fid, spt)
#     end
#     close(fid)
#     return nothing
# end
#
# """
#    extract a wavelet from the stress field
# """
# function multiStress2wlet(path::String , iz::Int64, ix::Int64)
#     fid = open(path, "r")
#     (nz, nx, ext, iflag, dt, nt) = InfoStress(path)
#     w = zeros(nt)
#     for it = 1 : nt
#         sts= readStress(path, it)
#         w[it] = sts.p[iz,ix]
#     end
#     return w
# end


text:
function compare(c1::Channel,
  c2::Channel;
  skip::Vector{Symbol}=[] )
#
TP = true
TP = TP && c1.state == c2.state
TP = TP && c1.sz_max == c2.sz_max
TP = TP && c1.data |> length == c2.data |> length
# exit early if tests already failed
!TP && (return false)
# now check contents of data
for i in 1:length(c1.data)
TP = TP && c1.data[i] == c2.data[i]
end
return TP
end

compare(a::Int, b::Int) = a == b
compare(a::Bool, b::Bool) = a == b
compare(a::Dict, b::Dict) = a == b


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

### Parsing utilities

_directives(::Type{DateFormat{S,T}}) where {S,T} = T.parameters

character_codes(df::Type{DateFormat{S,T}}) where {S,T} = character_codes(_directives(df))
function character_codes(directives::SimpleVector)
    letters = sizehint!(Char[], length(directives))
    for (i, directive) in enumerate(directives)
        if directive <: DatePart
            letter = first(directive.parameters)
            push!(letters, letter)
        end
    end
    return letters
end

genvar(t::DataType) = Symbol(lowercase(string(nameof(t))))

"""
    tryparsenext_core(str::AbstractString, pos::Int, len::Int, df::DateFormat, raise=false)

Parse the string according to the directives within the `DateFormat`. Parsing will start at
character index `pos` and will stop when all directives are used or we have parsed up to
the end of the string, `len`. When a directive cannot be parsed the returned value tuple
will be `nothing` if `raise` is false otherwise an exception will be thrown.

Return a 3-element tuple `(values, pos, num_parsed)`:
* `values::Union{Tuple, Nothing}`: Either `nothing`, or a tuple which contains a value
  for each `DatePart` within the `DateFormat` in the order
  in which they occur. If the string ends before we finish parsing all the directives
  the missing values will be filled in with default values.
* `pos::Int`: The character index at which parsing stopped.
* `num_parsed::Int`: The number of values which were parsed and stored within `values`.
  Useful for distinguishing parsed values from default values.
"""
@generated function tryparsenext_core(str::AbstractString, pos::Int, len::Int,
                                      df::DateFormat, raise::Bool=false)
    directives = _directives(df)
    letters = character_codes(directives)

    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]
    value_names = Symbol[genvar(t) for t in tokens]
    value_defaults = Tuple(CONVERSION_DEFAULTS[t] for t in tokens)

    # Pre-assign variables to defaults. Allows us to use `@goto done` without worrying about
    # unassigned variables.
    assign_defaults = Expr[]
    for (name, default) in zip(value_names, value_defaults)
        push!(assign_defaults, quote
            $name = $default
        end)
    end

    vi = 1
    parsers = Expr[]
    for i = 1:length(directives)
        if directives[i] <: DatePart
            name = value_names[vi]
            val = Symbol(:val, name)
            vi += 1
            push!(parsers, quote
                pos > len && @goto done
                $val, next_pos = tryparsenext(directives[$i], str, pos, len, locale)
                $val === nothing && @goto error
                $name = $val
                pos = next_pos
                num_parsed += 1
                directive_index += 1
            end)
        else
            push!(parsers, quote
                pos > len && @goto done
                delim, next_pos = tryparsenext(directives[$i], str, pos, len, locale)
                delim === nothing && @goto error
                pos = next_pos
                directive_index += 1
            end)
        end
    end

    quote
        directives = df.tokens
        locale::DateLocale = df.locale

        num_parsed = 0
        directive_index = 1

        $(assign_defaults...)
        $(parsers...)

        pos > len || @goto error

        @label done
        return $(Expr(:tuple, value_names...)), pos, num_parsed

        @label error
        if raise
            if directive_index > length(directives)
                throw(ArgumentError("Found extra characters at the end of date time string"))
            else
                d = directives[directive_index]
                throw(ArgumentError("Unable to parse date time. Expected directive $d at char $pos"))
            end
        end
        return nothing, pos, 0
    end
end

"""
    tryparsenext_internal(::Type{<:TimeType}, str, pos, len, df::DateFormat, raise=false)

Parse the string according to the directives within the `DateFormat`. The specified `TimeType`
type determines the type of and order of tokens returned. If the given `DateFormat` or string
does not provide a required token a default value will be used. When the string cannot be
parsed the returned value tuple will be `nothing` if `raise` is false otherwise an exception will
be thrown.

Return a 2-element tuple `(values, pos)`:
* `values::Union{Tuple, Nothing}`: Either `nothing`, or a tuple which contains a value
  for each token as specified by the passed in type.
* `pos::Int`: The character index at which parsing stopped.
"""
@generated function tryparsenext_internal(::Type{T}, str::AbstractString, pos::Int, len::Int,
                                          df::DateFormat, raise::Bool=false) where T<:TimeType
    letters = character_codes(df)

    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]
    value_names = Symbol[genvar(t) for t in tokens]

    output_tokens = CONVERSION_TRANSLATIONS[T]
    output_names = Symbol[genvar(t) for t in output_tokens]
    output_defaults = Tuple(CONVERSION_DEFAULTS[t] for t in output_tokens)

    # Pre-assign output variables to defaults. Ensures that all output variables are
    # assigned as the value tuple returned from `tryparsenext_core` may not include all
    # of the required variables.
    assign_defaults = Expr[
        quote
            $name = $default
        end
        for (name, default) in zip(output_names, output_defaults)
    ]

    # Unpacks the value tuple returned by `tryparsenext_core` into separate variables.
    value_tuple = Expr(:tuple, value_names...)

    quote
        values, pos, num_parsed = tryparsenext_core(str, pos, len, df, raise)
        values === nothing && return nothing, pos
        $(assign_defaults...)
        $value_tuple = values
        return $(Expr(:tuple, output_names...)), pos
    end
end

@inline function tryparsenext_base10(str::AbstractString, i::Int, len::Int, min_width::Int=1, max_width::Int=0)
    i > len && (return nothing, i)
    min_pos = min_width <= 0 ? i : i + min_width - 1
    max_pos = max_width <= 0 ? len : min(i + max_width - 1, len)
    d::Int64 = 0
    @inbounds while i <= max_pos
        c, ii = next(str, i)
        if '0' <= c <= '9'
            d = d * 10 + (c - '0')
        else
            break
        end
        i = ii
    end
    if i <= min_pos
        return nothing, i
    else
        return d, i
    end
end

@inline function tryparsenext_word(str::AbstractString, i, len, locale, maxchars=0)
    word_start, word_end = i, 0
    max_pos = maxchars <= 0 ? len : min(len, nextind(str, i, maxchars-1))
    @inbounds while i <= max_pos
        c, ii = next(str, i)
        if isalpha(c)
            word_end = i
        else
            break
        end
        i = ii
    end
    if word_end == 0
        return nothing, i
    else
        return SubString(str, word_start, word_end), i
    end
end

function Base.parse(::Type{DateTime}, s::AbstractString, df::typeof(ISODateTimeFormat))
    i, end_pos = firstindex(s), lastindex(s)

    dm = dd = Int64(1)
    th = tm = ts = tms = Int64(0)

    val, i = tryparsenext_base10(s, i, end_pos, 1)
    dy = val === nothing ? (@goto error) : val
    i > end_pos && @goto error

    c, i = next(s, i)
    c != '-' && @goto error
    i > end_pos && @goto done

    val, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    dm = val === nothing ? (@goto error) : val
    i > end_pos && @goto done

    c, i = next(s, i)
    c != '-' && @goto error
    i > end_pos && @goto done

    val, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    dd = val === nothing ? (@goto error) : val
    i > end_pos && @goto done

    c, i = next(s, i)
    c != 'T' && @goto error
    i > end_pos && @goto done

    val, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    th = val === nothing ? (@goto error) : val
    i > end_pos && @goto done

    c, i = next(s, i)
    c != ':' && @goto error
    i > end_pos && @goto done

    val, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    tm = val === nothing ? (@goto error) : val
    i > end_pos && @goto done

    c, i = next(s, i)
    c != ':' && @goto error
    i > end_pos && @goto done

    val, i = tryparsenext_base10(s, i, end_pos, 1, 2)
    ts = val === nothing ? (@goto error) : val
    i > end_pos && @goto done

    c, i = next(s, i)
    c != '.' && @goto error
    i > end_pos && @goto done

    val, j = tryparsenext_base10(s, i, end_pos, 1, 3)
    tms = val === nothing ? (@goto error) : val
    tms *= 10 ^ (3 - (j - i))

    j > end_pos || @goto error

    @label done
    return DateTime(dy, dm, dd, th, tm, ts, tms)

    @label error
    throw(ArgumentError("Invalid DateTime string"))
end

function Base.parse(::Type{T}, str::AbstractString, df::DateFormat=default_format(T)) where T<:TimeType
    pos, len = firstindex(str), lastindex(str)
    values, pos = tryparsenext_internal(T, str, pos, len, df, true)
    T(values...)
end

function Base.tryparse(::Type{T}, str::AbstractString, df::DateFormat=default_format(T)) where T<:TimeType
    pos, len = firstindex(str), lastindex(str)
    values, pos = tryparsenext_internal(T, str, pos, len, df, false)
    if values === nothing
        nothing
    elseif validargs(T, values...) === nothing
        # TODO: validargs gets called twice, since it's called again in the T constructor
        T(values...)
    else
        nothing
    end
end

"""
    parse_components(str::AbstractString, df::DateFormat) -> Array{Any}

Parse the string into its components according to the directives in the `DateFormat`.
Each component will be a distinct type, typically a subtype of Period. The order of the
components will match the order of the `DatePart` directives within the `DateFormat`. The
number of components may be less than the total number of `DatePart`.
"""
@generated function parse_components(str::AbstractString, df::DateFormat)
    letters = character_codes(df)
    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]

    quote
        pos, len = firstindex(str), lastindex(str)
        values, pos, num_parsed = tryparsenext_core(str, pos, len, df, true)
        t = values
        types = $(Expr(:tuple, tokens...))
        result = Vector{Any}(uninitialized, num_parsed)
        for (i, typ) in enumerate(types)
            i > num_parsed && break
            result[i] = typ(t[i])  # Constructing types takes most of the time
        end
        return result
    end
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
# This file holds the definition of the Recruitment type. This type is used to
#   describe a recruitment scheme for a specific part of the organisation and
#   defines which attributes need to be initialised and what the distributions
#   of those attributes are.

# The Recruitment type requires no additional types.
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export Recruitment
type Recruitment

    # The name of this recruitment type.
    name::String

    # The time between two recruitment drives.
    recruitFreq::Float64

    # The offset of the recruitment cycle.
    recruitOffset::Float64

    # Vary recruitment according to needs.
    isAdaptive::Bool

    # The min/max number of personnel members to recruit.
    minRecruit::Int
    maxRecruit::Int

    # The state to recruit into.
    recState::String

    # The type of the recruitment distribution.
    recDistType::Symbol

    # The nodes of the recruitment distribution, if applicable.
    recDistNodes::Dict{Float64, Float64}

    # This is the type of the age distribution.
    ageDistType::Symbol

    # The distribution function of the number of recruits in a single cycle.
    recDist::Function

    # These are the nodes (with probability weights) of the age distribution.
    ageDistNodes::Dict{Float64, Float64}

    # The distribution function of the age of a single recruit.
    ageDist::Function

    function Recruitment( name::String, freq::T1, offset::T2 = 0.0 ) where T1 <: Real where T2 <: Real

        if freq <= 0.0
            error( "Recruitment frequency must be > 0.0." )
        end  # if freq <= 0.0

        newRec = new()
        newRec.name = name
        newRec.recruitFreq = freq
        newRec.recruitOffset = offset % freq + ( offset < 0.0 ? freq : 0.0 )
        newRec.isAdaptive = true
        newRec.minRecruit = 0
        newRec.maxRecruit = typemax( Int )
        newRec.recState = ""
        newRec.recDistType = :disc
        newRec.recDistNodes = Dict{Float64, Float64}()
        newRec.recDist = function() return 0 end
        newRec.ageDistType = :disc
        newRec.ageDistNodes = Dict( 0.0 => 1.0 )
        newRec.ageDist = function() return 0.0 end
        return newRec

    end  # Recruitment( freq, offset )

end  # type Recruitment