text:
function update_changepoint_model!(m, σ, d, boundaries, np)
	@inbounds for i=1:np+1
		r = boundaries[i]:boundaries[i+1]
		for col = 1:size(d,2)
			m[r,col] .= nanmean(view(d,r,col))
			σ[r,col] .= nanstd(view(d,r,col))
		end
	end
end

function update_changepoint_mu!(m, d, boundaries, np)
	@inbounds for i=1:np+1
		r = boundaries[i]:boundaries[i+1]
		for col = 1:size(d,2)
			m[r,col] .= nanmean(view(d,r,col))
		end
	end
end

# # Not currently used
# function update_changepoint_sigma!(σ, d, boundaries, np)
# 	@inbounds for i=1:np+1
# 		r = boundaries[i]:boundaries[i+1]
# 		for col = 1:size(d,2)
# 			σ[r,col] .= nanstd(view(d,r,col))
# 		end
# 	end
# end

"""
```julia
changepoint(data, [sigma], nsteps; np, npmin, npmax)
```
Given an ordered array of `data` points, optionally with uncertainties `sigma`,
use a Markov chain Monte Carlo approach based on that of Gallagher et al., 2010
(10.1016/j.epsl.2011.09.015) to estimate the position (by index) and optionally
number of changepoints that best explain the `data`. Will return the results for
`nsteps` steps of the Markov chain.

Optionally, you may also specify as keyword arguments either `np` or
`npmin` and `npmax` to constrain the number of allowed changepoints.

Currently prints results to the terminal (stdout), one line per step of the Markov
chain, but a more efficent output format is probably desirable in a future version
of this function.
"""
function changepoint(data::AbstractArray, nsteps::Integer; np::Integer=0, npmin::Integer=0, npmax::Integer=min(size(data,1) ÷ 2, 11))

    MOVE = 0.70
    BIRTH = 0.15
    DEATH = 0.15

	DEBUG = false
	FORMATTED = true

	T = float(eltype(data))
	nrows = size(data,1)
	ncolumns = size(data,2)
	m = similar(data, T)
	σ = similar(data, T) #Array{T}(undef, ncolumns)
	σₚ = similar(data, T) #Array{T}(undef, ncolumns)

	# Number of possible changepoint locations
	K = nrows-1

	# Parse provided options
	if 0 < np <= K
		# If np is specified, use that
		npmin = npmax = np
	else
		# Otherwise, ensure all provided values are plausible and go with that
		npmax > K && (npmax = K)
		npmin < 0 && (npmin = 0)
		np = min(max(npmin, 2), npmax)
	end

	# Allocate output array of changepoints
	result = fill(0, nsteps, npmax)

	# Create and fill initial boundary point array
	boundaries = Array{Int}(undef, K+2)
	boundaries[1] = 1
	boundaries[np+2] = nrows
	boundaries[2:np+1] .= rand(2:nrows-1, np)
	boundariesₚ = similar(boundaries)
	boundary_sigma = nrows/np

	np = count_unique!(view(boundaries,1:np+2)) - 2

	# Calculate initial proposal and log likelihood
	update_changepoint_model!(m, σ, data, boundaries, np)
	ll = normpdf_ll(m, σ, data)

	# The actual loop
	@inbounds for i = 1:nsteps

		# Randomly choose a type of modification to the model
		r = rand()
		u = rand()

		# Update the model with the chosen modification
		if r < MOVE && np>0
			# Move a changepoint
			copyto!(boundariesₚ,1,boundaries,1,np+2)
			# Pick which changepoint to move
			pick = rand(2:np+1)
			# Move the changepoint
			boundary_adj = randn()*boundary_sigma
			boundary_prop = boundariesₚ[pick] + round(Int, boundary_adj)

			# Treat ends of array as periodic boundary conditions
			boundariesₚ[pick] = mod(boundary_prop - 1, nrows-1) + 1

			# Check if this has caused any redundancies
			npₚ = count_unique!(view(boundariesₚ,1:np+2)) - 2

			# Update the model
			update_changepoint_mu!(m, data, boundariesₚ, npₚ)

			# Calculate log likelihood for proposal
			if (1 < boundary_prop < nrows)
				llₚ = normpdf_ll(m, σ, data)
			else
				llₚ = -Inf
			end

			DEBUG && println("Move: llₚ-ll = $llₚ - $ll")
			if log(u) < llₚ-ll
				DEBUG && println("Accepted!")
				ll = llₚ
				boundary_sigma = abs(boundary_adj)*2.9
				# println("sigma: $boundary_sigma")
				copyto!(boundaries,1,boundariesₚ,1,np+2)
				# for n=1:np
				# 	print("$(boundariesₚ[n+1]),")
 				# end
				# FORMATTED && print("\n")
			end

		elseif r < MOVE+BIRTH
			# Add a changepoint
			if np < npmax
				copyto!(boundariesₚ,1,boundaries,1,np+2)

				# Propose a new changepoint
				boundariesₚ[np+3] = rand(2:nrows-1)
				npₚ = count_unique!(view(boundariesₚ,1:np+3)) - 2

				# Update the model
				# update_changepoint_model!(m, σ, data, boundariesₚ, npₚ)
				update_changepoint_mu!(m, data, boundariesₚ, npₚ)

				# Calculate log likelihood for proposal
				lqz = sum(1 ./ (2*σ.*σ))
				llₚ = normpdf_ll(m, σ, data)
				DEBUG && println("Birth: -lqz+llₚ-ll = $(-lqz) + $llₚ - $ll")
				if log(u) < llₚ-lqz-ll
					DEBUG && println("Accepted!")
					ll = llₚ
					np = npₚ
					copyto!(boundaries,1,boundariesₚ,1,np+2)
					# for n=1:np
					# 	print("$(boundariesₚ[n+1]),")
					# end
					# FORMATTED && print("\n")
				end
			end
		elseif r < MOVE+BIRTH+DEATH
			# Delete a changepoint
			if np > npmin
				copyto!(boundariesₚ,1,boundaries,1,np+2)

				# Pick which changepoint to delete
				pick = rand(2:np+1)
				boundariesₚ[pick]=boundariesₚ[pick+1]
				npₚ = count_unique!(view(boundariesₚ,1:np+2)) - 2

				# Update the model
				# update_changepoint_model!(m, σ, data, boundariesₚ, npₚ)
				update_changepoint_mu!(m, data, boundariesₚ, npₚ)

				# Calculate log likelihood for proposal
				llₚ = normpdf_ll(m, σ, data)

				lqz = sum(1 ./ (2*σ.*σ))
				DEBUG && println("Death: lqz+llₚ-ll = $lqz + $llₚ - $ll")
				if log(u) < llₚ+lqz-ll
					DEBUG && println("Accepted!")
					ll = llₚ
					np = npₚ
					copyto!(boundaries,1,boundariesₚ,1,np+2)
					# for n=1:np
					# 	print("$(boundariesₚ[n+1]),")
					# end
					# FORMATTED && print("\n")
				end
			end
		end

		# Record results
		result[i, 1:np] .= boundaries[2:(np+1)]
	end
	return result
end
function changepoint(data::AbstractArray, sigma::AbstractArray, nsteps::Integer; np::Integer=0, npmin::Integer=0, npmax::Integer=min(size(data,1) ÷ 2, 11))

    MOVE = 0.70
    BIRTH = 0.15
    DEATH = 0.15

	DEBUG = false
	FORMATTED = true

	T = float(eltype(data))
	nrows = size(data,1)
	ncolumns = size(data,2)
	m = similar(data, T)

	# Number of possible changepoint locations
	K = nrows-1

	# Parse provided options
	if 0 < np <= K
		# If np is specified, use that
		npmin = npmax = np
	else
		# Otherwise, ensure all provided values are plausible and go with that
		npmax > K && (npmax = K)
		npmin < 0 && (npmin = 0)
		np = min(max(npmin, 2), npmax)
	end

	# Allocate output array of changepoints
	result = fill(0, nsteps, npmax)

	# Create and fill initial boundary point array
	boundaries = Array{Int}(undef, K+2)
	boundaries[1] = 1
	boundaries[np+2] = nrows
	boundaries[2:np+1] .= rand(2:nrows-1, np)
	boundariesₚ = similar(boundaries)
	boundary_sigma = nrows/np

	np = count_unique!(view(boundaries,1:np+2)) - 2

	# Calculate initial proposal and log likelihood
	update_changepoint_mu!(m, data, boundaries, np)
	ll = normpdf_ll(m, sigma, data)

	# The actual loop
	@inbounds for i = 1:nsteps

		# Randomly choose a type of modification to the model
		r = rand()
		u = rand()

		# Update the model with the chosen modification
		if r < MOVE && np>0
			# Move a changepoint
			copyto!(boundariesₚ,1,boundaries,1,np+2)
			# Pick which changepoint to move
			pick = rand(2:np+1)
			# Move the changepoint
			boundary_adj = randn()*boundary_sigma
			boundary_prop = boundariesₚ[pick] + round(Int, boundary_adj)

			# Treat ends of array as periodic boundary conditions
			boundariesₚ[pick] = mod(boundary_prop - 1, nrows-1) + 1

			# Check if this has caused any redundancies
			npₚ = count_unique!(view(boundariesₚ,1:np+2)) - 2

			# Update the model
			update_changepoint_mu!(m, data, boundariesₚ, npₚ)

			# Calculate log likelihood for proposal
			if (1 < boundary_prop < nrows)
				llₚ = normpdf_ll(m, sigma, data)
			else
				llₚ = -Inf
			end

			DEBUG && println("Move: llₚ-ll = $llₚ - $ll")
			if log(u) < llₚ-ll
				DEBUG && println("Accepted!")
				ll = llₚ
				boundary_sigma = abs(boundary_adj)*2.9
				# println("sigma: $boundary_sigma")
				copyto!(boundaries,1,boundariesₚ,1,np+2)
				# for n=1:np
				# 	print("$(boundariesₚ[n+1]),")
				# end
				# FORMATTED && print("\n")
			end
		elseif r < MOVE+BIRTH
			# Add a changepoint
			if np < npmax
				copyto!(boundariesₚ,1,boundaries,1,np+2)

				# Propose a new changepoint
				boundariesₚ[np+3] = rand(2:nrows-1)
				npₚ = count_unique!(view(boundariesₚ,1:np+3)) - 2

				# Update the model
				update_changepoint_mu!(m, data, boundariesₚ, npₚ)

				# Calculate log likelihood for proposal
				lqz = sum(1 ./ (2*sigma.*sigma))
				llₚ = normpdf_ll(m, sigma, data)
				DEBUG && println("Birth: -lqz+llₚ-ll = $(-lqz) + $llₚ - $ll")
				if log(u) < llₚ-lqz-ll
					DEBUG && println("Accepted!")
					ll = llₚ
					np = npₚ
					copyto!(boundaries,1,boundariesₚ,1,np+2)
					# for n=1:np
					# 	print("$(boundariesₚ[n+1]),")
					# end
					# FORMATTED && print("\n")
				end
			end
		elseif r < MOVE+BIRTH+DEATH
			# Delete a changepoint
			if np > npmin
				copyto!(boundariesₚ,1,boundaries,1,np+2)

				# Pick which changepoint to delete
				pick = rand(2:np+1)
				boundariesₚ[pick]=boundariesₚ[pick+1]
				npₚ = count_unique!(view(boundariesₚ,1:np+2)) - 2

				# Update the model
				update_changepoint_mu!(m, data, boundariesₚ, npₚ)

				# Calculate log likelihood for proposal
				llₚ = normpdf_ll(m, sigma, data)

				lqz = sum(1 ./ (2*sigma.*sigma))
				DEBUG && println("Death: lqz+llₚ-ll = $lqz + $llₚ - $ll")
				if log(u) < llₚ+lqz-ll
					DEBUG && println("Accepted!")
					ll = llₚ
					np = npₚ
					copyto!(boundaries,1,boundariesₚ,1,np+2)
					# for n=1:np
					# 	print("$(boundariesₚ[n+1]),")
					# end
					# FORMATTED && print("\n")
				end
			end
		end

		# Record results
		result[i, 1:np] .= boundaries[2:(np+1)]
	end
	return result
end
export changepoint


text:
include("./load_csv.jl")

# Specify model parameters and/or initial values to optimize
function get_search_index()::Tuple{Array{Int64,1},Array{Int64,1}}
    # parameters
    search_idx_params::Vector{Int} = [
        C.kf1,
        C.kr1,
        C.kf2,
        C.kr2,
        C.kf3,
        C.kr3,
        C.kf4,
        C.kr4,
        C.kf5,
        C.kr5,
        C.kf6,
        C.kr6,
        C.kf7,
        C.kr7,
        C.kf8,
        C.kr8,
        C.kf9,
        C.kr9,
        C.kf10,
        C.kr10,
        C.kf11,
        C.kr11,
        C.kf12,
        C.kr12,
        C.kf13,
        C.kr13,
        C.kf14,
        C.kr14,
        C.kf15,
        C.kr15,
        C.kf16,
        C.kr16,
        C.kf17,
        C.kr17,
        C.kf18,
        C.kr18,
        C.kf19,
        C.kr19,
        C.kf20,
        C.kr20,
        C.kf21,
        C.kr21,
        C.kf22,
        C.kr22,
        C.kf23,
        C.kr23,
        C.kf24,
        C.kr24,
        C.kf25,
        C.kr25,
        C.kf26,
        C.kr26,
        C.kf27,
        C.kr27,
        C.kf28,
        C.kr28,
        C.kf29,
        C.kr29,
        C.kf30,
        C.kr30,
        C.kf31,
        C.kr31,
        C.kf32,
        C.kr32,
        C.kf33,
        C.kr33,
        C.kf34,
        C.kr34,
        C.kf35,
        C.kr35,
        C.kf36,
        C.kr36,
        C.kf37,
        C.kr37,
        C.kf38,
        C.kr38,
        # C.kf39,
        # C.kr39,
        # C.kf40,
        # C.kr40,
        # C.kf41,
        # C.kr41,
        # C.kf42,
        # C.kr42,
        # C.kf43,
        # C.kr43,
        C.kf44,
        # C.kr44,
        # C.kf45,
        # C.kr45,
        # C.kf46,
        # C.kr46,
        # C.kf47,
        # C.kr47,
        # C.kf48,
        # C.kr48,
        # C.kf49,
        # C.kr49,
        C.kf50,
        C.kr50,
        # C.kf51,
        # C.kr51,
        # C.kf52,
        # C.kr52,
        # C.kf53,
        # C.kr53,
        # C.kf54,
        # C.kr54,
        # C.kf55,
        # C.kr55,
        C.kf56,
        # C.kr56,
        # C.kf57,
        # C.kr57,
        # C.kf58,
        # C.kr58,
        # C.kf59,
        # C.kr59,
        # C.kf60,
        # C.kr60,
        # C.kf61,
        # C.kr61,
        C.kf62,
        C.kr62,
        # C.kf63,
        # C.kr63,
        # C.kf64,
        # C.kr64,
        # C.kf65,
        # C.kr65,
        # C.kf66,
        # C.kr66,
        # C.kf67,
        # C.kr67,
        # C.kf68,
        # C.kr68,
        # C.kf69,
        # C.kr69,
        # C.kf70,
        # C.kr70,
        # C.kf71,
        # C.kr71,
        # C.kf72,
        # C.kr72,
        # C.kf73,
        # C.kr73,
        C.kf74,
        C.kr74,
        # C.kf75,
        # C.kr75,
        # C.kf76,
        # C.kr76,
        # C.kf77,
        # C.kr77,
        # C.kf78,
        # C.kr78,
        # C.kf79,
        # C.kr79,
        # C.kf80,
        # C.kr80,
        # C.kf81,
        # C.kr81,
        # C.kf82,
        # C.kr82,
        # C.kf83,
        # C.kr83,
        # C.kf84,
        # C.kr84,
        # C.kf85,
        # C.kr85,
        C.kf86,
        C.kr86,
        # C.kf87,
        # C.kr87,
        # C.kf88,
        # C.kr88,
        # C.kf89,
        # C.kr89,
        # C.kf90,
        # C.kr90,
        # C.kf91,
        # C.kr91,
        # C.kf92,
        # C.kr92,
        # C.kf93,
        # C.kr93,
        # C.kf94,
        # C.kr94,
        # C.kf95,
        # C.kr95,
        # C.kf96,
        # C.kr96,
        # C.kf97,
        # C.kr97,
        C.kf98,
        # C.kr98,
        # C.kf99,
        # C.kr99,
        # C.kf100,
        # C.kr100,
        # C.kf101,
        # C.kr101,
        # C.kf102,
        # C.kr102,
        # C.kf103,
        # C.kr103,
        # C.kf104,
        # C.kr104,
        # C.kf105,
        # C.kr105,
        # C.kf106,
        # C.kr106,
        # C.kf107,
        # C.kr107,
        # C.kf108,
        # C.kr108,
        # C.kf109,
        # C.kr109,
        C.V110,
        C.K110,
        C.V111,
        C.K111,
        C.V112,
        C.K112,
        C.V113,
        C.K113,
        C.V114,
        C.K114,
        C.V115,
        C.K115,
        C.V116,
        C.K116,
        C.V117,
        C.K117,
        C.V118,
        C.K118,
        C.V119,
        C.K119,
        C.V120,
        C.K120,
        C.V121,
        C.K121,
        C.kf122,
        # C.kr122,
        # C.kf123,
        # C.kr123,
        # C.kf124,
        # C.kr124,
        # C.kf125,
        # C.kr125,
        # C.kf126,
        # C.kr126,
        # C.kf127,
        # C.kr127,
        # C.kf128,
        # C.kr128,
        # C.kf129,
        # C.kr129,
        # C.kf130,
        # C.kr130,
        # C.kf131,
        # C.kr131,
        # C.kf132,
        # C.kr132,
        # C.kf133,
        # C.kr133,
        # C.kf134,
        # C.kr134,
        C.kf135,
        C.kr135,
        # C.kf136,
        # C.kr136,
        # C.kf137,
        # C.kr137,
        # C.kf138,
        # C.kr138,
        # C.kf139,
        # C.kr139,
        # C.kf140,
        # C.kr140,
        # C.kf141,
        # C.kr141,
        # C.kf142,
        # C.kr142,
        # C.kf143,
        # C.kr143,
        # C.kf144,
        # C.kr144,
        # C.kf145,
        # C.kr145,
        # C.kf146,
        # C.kr146,
        # C.kf147,
        # C.kr147,
        C.kf148,
        # C.kr148,
        # C.kf149,
        # C.kr149,
        # C.kf150,
        # C.kr150,
        # C.kf151,
        # C.kr151,
        # C.kf152,
        # C.kr152,
        # C.kf153,
        # C.kr153,
        # C.kf154,
        # C.kr154,
        # C.kf155,
        # C.kr155,
        # C.kf156,
        # C.kr156,
        # C.kf157,
        # C.kr157,
        # C.kf158,
        # C.kr158,
        # C.kf159,
        # C.kr159,
        # C.kf160,
        # C.kr160,
        C.kf161,
        C.kr161,
        # C.kf162,
        # C.kr162,
        # C.kf163,
        # C.kr163,
        # C.kf164,
        # C.kr164,
        # C.kf165,
        # C.kr165,
        # C.kf166,
        # C.kr166,
        # C.kf167,
        # C.kr167,
        # C.kf168,
        # C.kr168,
        # C.kf169,
        # C.kr169,
        # C.kf170,
        # C.kr170,
        # C.kf171,
        # C.kr171,
        # C.kf172,
        # C.kr172,
        # C.kf173,
        # C.kr173,
        C.kf174,
        C.kr174,
        # C.kf175,
        # C.kr175,
        # C.kf176,
        # C.kr176,
        # C.kf177,
        # C.kr177,
        # C.kf178,
        # C.kr178,
        # C.kf179,
        # C.kr179,
        # C.kf180,
        # C.kr180,
        # C.kf181,
        # C.kr181,
        # C.kf182,
        # C.kr182,
        # C.kf183,
        # C.kr183,
        # C.kf184,
        # C.kr184,
        # C.kf185,
        # C.kr185,
        # C.kf186,
        # C.kr186,
        # C.kf187,
        # C.kr187,
        # C.kf188,
        # C.kr188,
        # C.kf189,
        # C.kr189,
        # C.kf190,
        # C.kr190,
        # C.kf191,
        # C.kr191,
        # C.kf192,
        # C.kr192,
        C.kf193,
        # C.kr193,
        # C.kf194,
        # C.kr194,
        # C.kf195,
        # C.kr195,
        # C.kf196,
        # C.kr196,
        # C.kf197,
        # C.kr197,
        # C.kf198,
        # C.kr198,
        # C.kf199,
        # C.kr199,
        # C.kf200,
        # C.kr200,
        # C.kf201,
        # C.kr201,
        # C.kf202,
        # C.kr202,
        # C.kf203,
        # C.kr203,
        # C.kf204,
        # C.kr204,
        # C.kf205,
        # C.kr205,
        # C.kf206,
        # C.kr206,
        # C.kf207,
        # C.kr207,
        # C.kf208,
        # C.kr208,
        # C.kf209,
        # C.kr209,
        # C.kf210,
        # C.kr210,
        # C.kf211,
        # C.kr211,
        C.kf212,
        C.kr212,
        # C.kf213,
        # C.kr213,
        # C.kf214,
        # C.kr214,
        # C.kf215,
        # C.kr215,
        # C.kf216,
        # C.kr216,
        # C.kf217,
        # C.kr217,
        # C.kf218,
        # C.kr218,
        # C.kf219,
        # C.kr219,
        # C.kf220,
        # C.kr220,
        # C.kf221,
        # C.kr221,
        # C.kf222,
        # C.kr222,
        # C.kf223,
        # C.kr223,
        # C.kf224,
        # C.kr224,
        C.kf225,
        C.kr225,
        # C.kf226,
        # C.kr226,
        # C.kf227,
        # C.kr227,
        # C.kf228,
        # C.kr228,
        # C.kf229,
        # C.kr229,
        # C.kf230,
        # C.kr230,
        # C.kf231,
        # C.kr231,
        # C.kf232,
        # C.kr232,
        # C.kf233,
        # C.kr233,
        # C.kf234,
        # C.kr234,
        # C.kf235,
        # C.kr235,
        # C.kf236,
        # C.kr236,
        # C.kf237,
        # C.kr237,
        # C.kf238,
        # C.kr238,
        # C.kf239,
        # C.kr239,
        # C.kf240,
        # C.kr240,
        # C.kf241,
        # C.kr241,
        C.kf242,
        # C.kr242,
        # C.kf243,
        # C.kr243,
        # C.kf244,
        # C.kr244,
        # C.kf245,
        # C.kr245,
        # C.kf246,
        # C.kr246,
        # C.kf247,
        # C.kr247,
        # C.kf248,
        # C.kr248,
        # C.kf249,
        # C.kr249,
        # C.kf250,
        # C.kr250,
        # C.kf251,
        # C.kr251,
        # C.kf252,
        # C.kr252,
        # C.kf253,
        # C.kr253,
        # C.kf254,
        # C.kr254,
        # C.kf255,
        # C.kr255,
        # C.kf256,
        # C.kr256,
        # C.kf257,
        # C.kr257,
        # C.kf258,
        # C.kr258,
        C.kf259,
        C.kr259,
        C.V260,
        C.K260,
        C.V261,
        C.K261,
        C.V262,
        C.K262,
        C.V263,
        C.K263,
        C.V264,
        C.K264,
        C.V265,
        C.K265,
        C.V266,
        C.K266,
        C.kf267,
        C.kf268,
        # C.kf269,
        # C.kf270,
        # C.kf271,
        # C.kf272,
        # C.kf273,
        # C.kf274,
        # C.kf275,
        # C.kf276,
        # C.kf277,
        # C.kf278,
        # C.kf279,
        # C.kf280,
        # C.kf281,
        # C.kf282,
        # C.kf283,
        # C.kf284,
        # C.kf285,
        # C.kf286,
        # C.kf287,
        # C.kf288,
        # C.kf289,
        # C.kf290,
        C.V291,
        C.K291,
        # C.n291,
        C.kf292,
        # C.kr292,
        C.kf293,
        C.kf294,
        C.V295,
        C.K295,
        C.kf296,
        C.kf297,
        C.kf298,
        C.kr298,
        C.kf299,
        C.kr299,
        C.kf300,
        C.kr300,
        C.kf301,
        C.kr301,
        C.kf302,
        C.kr302,
        # C.kf303,
        # C.kr303,
        # C.kf304,
        # C.kr304,
        # C.kf305,
        # C.kr305,
        # C.kf306,
        # C.kr306,
        # C.kf307,
        # C.kr307,
        C.V308,
        C.K308,
        C.V309,
        C.K309,
        C.V310,
        C.K310,
        # C.n310,
        C.KF310,
        # C.nF310,
        C.kf311,
        # C.kr311,
        C.kf312,
        C.kf313,
        C.V314,
        C.K314,
        C.V315,
        C.K315,
        C.V316,
        C.K316,
        C.V317,
        C.K317,
        C.kf318,
        C.kf319,
        #
        C.w_EGFR,
        C.w_ERBB2,
        C.w_ERBB3,
        C.w_ERBB4,
        #
        C.w_GRB2,
        #
        C.w_SHC1,
        C.w_SHC2,
        C.w_SHC3,
        C.w_SHC4,
        #
        C.w_RASA1,
        C.w_RASA2,
        C.w_RASA3,
        #
        C.w_PIK3CA,
        C.w_PIK3CB,
        C.w_PIK3CD,
        C.w_PIK3CG,
        #
        C.w_PTEN,
        #
        C.w_SOS1,
        C.w_SOS2,
        #
        C.w_GAB1,
        #
        C.w_HRAS,
        C.w_KRAS,
        C.w_NRAS,
        #
        C.w_ARAF,
        C.w_BRAF,
        C.w_RAF1,
        #
        C.w_MAP2K1,
        C.w_MAP2K2,
        #
        C.w_MAPK1,
        C.w_MAPK3,
        #
        C.w_AKT1,
        C.w_AKT2,
        #
        C.w_PTPN1,
        #
        C.w_GSK3B,
        #
        C.w_DUSP5,
        C.w_DUSP6,
        C.w_DUSP7,
        #
        C.w_MYC,
    ]

    # initial values
    search_idx_initials::Vector{Int} = [
        V.PIP2
    ]

    return search_idx_params, search_idx_initials
end

function get_search_region()::Matrix{Float64}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()
    search_param::Vector{Float64} = init_search_param(search_idx, p, u0)

    search_rgn::Matrix{Float64} = zeros(2, length(p) + length(u0))

    # Default: 0.1 ~ 10x
    for (i, j) in enumerate(search_idx[1])
        search_rgn[1,j] = 1e-4  # lower bound
        search_rgn[2,j] = 1e+5  # upper bound
    end

    # Default: 0.5 ~ 2x
    for (i, j) in enumerate(search_idx[2])
        search_rgn[1,j + length(p)] = search_param[i + length(search_idx[1])] * 1  # lower bound
        search_rgn[2,j + length(p)] = search_param[i + length(search_idx[1])] * 1000  # upper bound
    end

    # search_rgn[:, C.param_name] = [lower_bound, upper_bound]
    # search_rgn[:, V.var_name+length(p)] = [lower_bound, upper_bound]

    for input_layer in [
        C.kf1,
        C.kr1,
        C.kf2,
        C.kr2,
        C.kf3,
        C.kr3,
        C.kf4,
        C.kr4,
        C.kf5,
        C.kr5,
        C.kf6,
        C.kr6,
        C.kf7,
        C.kr7,
        C.kf8,
        C.kr8,
        C.kf9,
        C.kr9,
    ]
        search_rgn[:, input_layer] = [1e-2, 1e+2]
    end

    for receptor_phosphorylation in [
        C.kf10,
        C.kr10,
        C.kf11,
        C.kr11,
        C.kf12,
        C.kr12,
        C.kf13,
        C.kr13,
        C.kf14,
        C.kr14,
        C.kf15,
        C.kr15,
    ]
        search_rgn[:, receptor_phosphorylation] = [1e-2, 1e+2]
    end

    for recepotor_adaptor in [
        C.kf16,
        C.kr16,
        C.kf17,
        C.kr17,
        C.kf18,
        C.kr18,
        C.kf19,
        C.kr19,
        C.kf20,
        C.kr20,
        C.kf21,
        C.kr21,
        C.kf22,
        C.kr22,
        C.kf23,
        C.kr23,
        C.kf24,
        C.kr24,
        C.kf25,
        C.kr25,
        C.kf26,
        C.kr26,
        C.kf27,
        C.kr27,
        C.kf28,
        C.kr28,
        C.kf29,
        C.kr29,
        C.kf30,
        C.kr30,
        C.kf31,
        C.kr31,
        C.kf32,
        C.kr32,
        C.kf33,
        C.kr33,
        C.kf34,
        C.kr34,
        C.kf35,
        C.kr35,
        C.kf36,
        C.kr36,
        C.kf37,
        C.kr37,
    ]
        search_rgn[:, recepotor_adaptor] = [1e-4, 1e+2]
    end

    for degradation in [
        C.kf267,
        C.kf268,
        C.kf294,
        C.kf296,
        C.kf297,
        C.kf312,
        C.kf318,
        C.kf319,
    ]
        search_rgn[:, degradation] = [1e-4, 1e+0]
    end

    for EC50 in [
        C.K110,
        C.K111,
        C.K112,
        C.K113,
        C.K114,
        C.K115,
        C.K116,
        C.K117,
        C.K118,
        C.K119,
        C.K120,
        C.K121,
        C.K260,
        C.K261,
        C.K262,
        C.K263,
        C.K264,
        C.K265,
        C.K266,
        C.K291,
        C.K295,
        C.K308,
        C.K309,
        C.K310,
        C.K314,
        C.K315,
        C.K316,
        C.K317,
    ]
        search_rgn[:, EC50] = [1e+0, 1e+4]
    end

    search_rgn[:, C.w_EGFR] = [0.01, 100]
    search_rgn[:, C.w_ERBB2] = [0.01, 100]
    search_rgn[:, C.w_ERBB3] = [0.01, 100]
    search_rgn[:, C.w_ERBB4] = [0.01, 100]
    #
    search_rgn[:, C.w_GRB2] = [0.01, 100]
    #
    search_rgn[:, C.w_SHC1] = [0.01, 100]
    search_rgn[:, C.w_SHC2] = [0.01, 100]
    search_rgn[:, C.w_SHC3] = [0.01, 100]
    search_rgn[:, C.w_SHC4] = [0.01, 100]
    #
    search_rgn[:, C.w_RASA1] = [0.01, 100]
    search_rgn[:, C.w_RASA2] = [0.01, 100]
    search_rgn[:, C.w_RASA3] = [0.01, 100]
    #
    search_rgn[:, C.w_PIK3CA] = [0.01, 100]
    search_rgn[:, C.w_PIK3CB] = [0.01, 100]
    search_rgn[:, C.w_PIK3CD] = [0.01, 100]
    search_rgn[:, C.w_PIK3CG] = [0.01, 100]
    #
    search_rgn[:, C.w_PTEN] = [0.01, 100]
    #
    search_rgn[:, C.w_SOS1] = [0.01, 100]
    search_rgn[:, C.w_SOS2] = [0.01, 100]
    #
    search_rgn[:, C.w_GAB1] = [0.01, 100]
    #
    search_rgn[:, C.w_HRAS] = [0.01, 100]
    search_rgn[:, C.w_KRAS] = [0.01, 100]
    search_rgn[:, C.w_NRAS] = [0.01, 100]
    #
    search_rgn[:, C.w_ARAF] = [0.01, 100]
    search_rgn[:, C.w_BRAF] = [0.01, 100]
    search_rgn[:, C.w_RAF1] = [0.01, 100]
    #
    search_rgn[:, C.w_MAP2K1] = [0.01, 100]
    search_rgn[:, C.w_MAP2K2] = [0.01, 100]
    #
    search_rgn[:, C.w_MAPK1] = [0.01, 100]
    search_rgn[:, C.w_MAPK3] = [0.01, 100]
    #
    search_rgn[:, C.w_AKT1] = [0.01, 100]
    search_rgn[:, C.w_AKT2] = [0.01, 100]
    #
    # search_rgn[:, C.w_RPS6KA1] = [0.01, 100]
    # search_rgn[:, C.w_RPS6KA2] = [0.01, 100]
    # search_rgn[:, C.w_RPS6KA3] = [0.01, 100]
    # search_rgn[:, C.w_RPS6KA6] = [0.01, 100]
    #
    # search_rgn[:, C.w_PPP2CA] = [0.01, 100]
    # search_rgn[:, C.w_PPP2CB] = [0.01, 100]
    #
    search_rgn[:, C.w_PTPN1] = [0.01, 100]
    #
    # search_rgn[:, C.w_CREB1] = [0.01, 100]
    #
    # search_rgn[:, C.w_ELK1] = [0.01, 100]
    #
    search_rgn[:, C.w_GSK3B] = [0.01, 100]
    #
    search_rgn[:, C.w_DUSP5] = [0.01, 100]
    search_rgn[:, C.w_DUSP6] = [0.01, 100]
    search_rgn[:, C.w_DUSP7] = [0.01, 100]
    #
    # search_rgn[:, C.w_FOS] = [0.01, 100]
    #
    search_rgn[:, C.w_MYC] = [0.01, 100]

    search_rgn = conv_lin2log!(search_rgn, search_idx)

    return search_rgn
end

function update_param(indiv::Vector{Float64})::Tuple{Array{Float64,1},Array{Float64,1}}
    p::Vector{Float64} = param_values()
    u0::Vector{Float64} = initial_values()

    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()

    for (i, j) in enumerate(search_idx[1])
        @inbounds p[j] = indiv[i]
    end
    for (i, j) in enumerate(search_idx[2])
        @inbounds u0[j] = indiv[i + length(search_idx[1])]
    end

    (p, u0) = individualize!(p, u0, "MCF7_BREAST")

    # constraints --------------------------------------------------------------
    p[C.kf39] = p[C.kf38]
    p[C.kr39] = p[C.kr38]
    p[C.kf40] = p[C.kf38]
    p[C.kr40] = p[C.kr38]
    p[C.kf41] = p[C.kf38]
    p[C.kr41] = p[C.kr38]
    p[C.kf42] = p[C.kf38]
    p[C.kr42] = p[C.kr38]
    p[C.kf43] = p[C.kf38]
    p[C.kr43] = p[C.kr38]
    p[C.kf45] = p[C.kf44]
    p[C.kr45] = p[C.kr44]
    p[C.kf46] = p[C.kf44]
    p[C.kr46] = p[C.kr44]
    p[C.kf47] = p[C.kf44]
    p[C.kr47] = p[C.kr44]
    p[C.kf48] = p[C.kf44]
    p[C.kr48] = p[C.kr44]
    p[C.kf49] = p[C.kf44]
    p[C.kr49] = p[C.kr44]
    p[C.kf51] = p[C.kf50]
    p[C.kr51] = p[C.kr50]
    p[C.kf52] = p[C.kf50]
    p[C.kr52] = p[C.kr50]
    p[C.kf53] = p[C.kf50]
    p[C.kr53] = p[C.kr50]
    p[C.kf54] = p[C.kf50]
    p[C.kr54] = p[C.kr50]
    p[C.kf55] = p[C.kf50]
    p[C.kr55] = p[C.kr50]
    p[C.kf57] = p[C.kf56]
    p[C.kr57] = p[C.kr56]
    p[C.kf58] = p[C.kf56]
    p[C.kr58] = p[C.kr56]
    p[C.kf59] = p[C.kf56]
    p[C.kr59] = p[C.kr56]
    p[C.kf60] = p[C.kf56]
    p[C.kr60] = p[C.kr56]
    p[C.kf61] = p[C.kf56]
    p[C.kr61] = p[C.kr56]
    p[C.kf63] = p[C.kf62]
    p[C.kr63] = p[C.kr62]
    p[C.kf64] = p[C.kf62]
    p[C.kr64] = p[C.kr62]
    p[C.kf65] = p[C.kf62]
    p[C.kr65] = p[C.kr62]
    p[C.kf66] = p[C.kf62]
    p[C.kr66] = p[C.kr62]
    p[C.kf67] = p[C.kf62]
    p[C.kr67] = p[C.kr62]
    p[C.kf68] = p[C.kf62]
    p[C.kr68] = p[C.kr62]
    p[C.kf69] = p[C.kf62]
    p[C.kr69] = p[C.kr62]
    p[C.kf70] = p[C.kf62]
    p[C.kr70] = p[C.kr62]
    p[C.kf71] = p[C.kf62]
    p[C.kr71] = p[C.kr62]
    p[C.kf72] = p[C.kf62]
    p[C.kr72] = p[C.kr62]
    p[C.kf73] = p[C.kf62]
    p[C.kr73] = p[C.kr62]
    p[C.kf75] = p[C.kf74]
    p[C.kr75] = p[C.kr74]
    p[C.kf76] = p[C.kf74]
    p[C.kr76] = p[C.kr74]
    p[C.kf77] = p[C.kf74]
    p[C.kr77] = p[C.kr74]
    p[C.kf78] = p[C.kf74]
    p[C.kr78] = p[C.kr74]
    p[C.kf79] = p[C.kf74]
    p[C.kr79] = p[C.kr74]
    p[C.kf80] = p[C.kf74]
    p[C.kr80] = p[C.kr74]
    p[C.kf81] = p[C.kf74]
    p[C.kr81] = p[C.kr74]
    p[C.kf82] = p[C.kf74]
    p[C.kr82] = p[C.kr74]
    p[C.kf83] = p[C.kf74]
    p[C.kr83] = p[C.kr74]
    p[C.kf84] = p[C.kf74]
    p[C.kr84] = p[C.kr74]
    p[C.kf85] = p[C.kf74]
    p[C.kr85] = p[C.kr74]
    p[C.kf87] = p[C.kf86]
    p[C.kr87] = p[C.kr86]
    p[C.kf88] = p[C.kf86]
    p[C.kr88] = p[C.kr86]
    p[C.kf89] = p[C.kf86]
    p[C.kr89] = p[C.kr86]
    p[C.kf90] = p[C.kf86]
    p[C.kr90] = p[C.kr86]
    p[C.kf91] = p[C.kf86]
    p[C.kr91] = p[C.kr86]
    p[C.kf92] = p[C.kf86]
    p[C.kr92] = p[C.kr86]
    p[C.kf93] = p[C.kf86]
    p[C.kr93] = p[C.kr86]
    p[C.kf94] = p[C.kf86]
    p[C.kr94] = p[C.kr86]
    p[C.kf95] = p[C.kf86]
    p[C.kr95] = p[C.kr86]
    p[C.kf96] = p[C.kf86]
    p[C.kr96] = p[C.kr86]
    p[C.kf97] = p[C.kf86]
    p[C.kr97] = p[C.kr86]
    p[C.kf99] = p[C.kf98]
    p[C.kr99] = p[C.kr98]
    p[C.kf100] = p[C.kf98]
    p[C.kr100] = p[C.kr98]
    p[C.kf101] = p[C.kf98]
    p[C.kr101] = p[C.kr98]
    p[C.kf102] = p[C.kf98]
    p[C.kr102] = p[C.kr98]
    p[C.kf103] = p[C.kf98]
    p[C.kr103] = p[C.kr98]
    p[C.kf104] = p[C.kf98]
    p[C.kr104] = p[C.kr98]
    p[C.kf105] = p[C.kf98]
    p[C.kr105] = p[C.kr98]
    p[C.kf106] = p[C.kf98]
    p[C.kr106] = p[C.kr98]
    p[C.kf107] = p[C.kf98]
    p[C.kr107] = p[C.kr98]
    p[C.kf108] = p[C.kf98]
    p[C.kr108] = p[C.kr98]
    p[C.kf109] = p[C.kf98]
    p[C.kr109] = p[C.kr98]
    p[C.kf123] = p[C.kf122]
    p[C.kr123] = p[C.kr122]
    p[C.kf124] = p[C.kf122]
    p[C.kr124] = p[C.kr122]
    p[C.kf125] = p[C.kf122]
    p[C.kr125] = p[C.kr122]
    p[C.kf126] = p[C.kf122]
    p[C.kr126] = p[C.kr122]
    p[C.kf127] = p[C.kf122]
    p[C.kr127] = p[C.kr122]
    p[C.kf128] = p[C.kf122]
    p[C.kr128] = p[C.kr122]
    p[C.kf129] = p[C.kf122]
    p[C.kr129] = p[C.kr122]
    p[C.kf130] = p[C.kf122]
    p[C.kr130] = p[C.kr122]
    p[C.kf131] = p[C.kf122]
    p[C.kr131] = p[C.kr122]
    p[C.kf132] = p[C.kf122]
    p[C.kr132] = p[C.kr122]
    p[C.kf133] = p[C.kf122]
    p[C.kr133] = p[C.kr122]
    p[C.kf134] = p[C.kf122]
    p[C.kr134] = p[C.kr122]
    p[C.kf136] = p[C.kf135]
    p[C.kr136] = p[C.kr135]
    p[C.kf137] = p[C.kf135]
    p[C.kr137] = p[C.kr135]
    p[C.kf138] = p[C.kf135]
    p[C.kr138] = p[C.kr135]
    p[C.kf139] = p[C.kf135]
    p[C.kr139] = p[C.kr135]
    p[C.kf140] = p[C.kf135]
    p[C.kr140] = p[C.kr135]
    p[C.kf141] = p[C.kf135]
    p[C.kr141] = p[C.kr135]
    p[C.kf142] = p[C.kf135]
    p[C.kr142] = p[C.kr135]
    p[C.kf143] = p[C.kf135]
    p[C.kr143] = p[C.kr135]
    p[C.kf144] = p[C.kf135]
    p[C.kr144] = p[C.kr135]
    p[C.kf145] = p[C.kf135]
    p[C.kr145] = p[C.kr135]
    p[C.kf146] = p[C.kf135]
    p[C.kr146] = p[C.kr135]
    p[C.kf147] = p[C.kf135]
    p[C.kr147] = p[C.kr135]
    p[C.kf149] = p[C.kf148]
    p[C.kr149] = p[C.kr148]
    p[C.kf150] = p[C.kf148]
    p[C.kr150] = p[C.kr148]
    p[C.kf151] = p[C.kf148]
    p[C.kr151] = p[C.kr148]
    p[C.kf152] = p[C.kf148]
    p[C.kr152] = p[C.kr148]
    p[C.kf153] = p[C.kf148]
    p[C.kr153] = p[C.kr148]
    p[C.kf154] = p[C.kf148]
    p[C.kr154] = p[C.kr148]
    p[C.kf155] = p[C.kf148]
    p[C.kr155] = p[C.kr148]
    p[C.kf156] = p[C.kf148]
    p[C.kr156] = p[C.kr148]
    p[C.kf157] = p[C.kf148]
    p[C.kr157] = p[C.kr148]
    p[C.kf158] = p[C.kf148]
    p[C.kr158] = p[C.kr148]
    p[C.kf159] = p[C.kf148]
    p[C.kr159] = p[C.kr148]
    p[C.kf160] = p[C.kf148]
    p[C.kr160] = p[C.kr148]
    p[C.kf162] = p[C.kf161]
    p[C.kr162] = p[C.kr161]
    p[C.kf163] = p[C.kf161]
    p[C.kr163] = p[C.kr161]
    p[C.kf164] = p[C.kf161]
    p[C.kr164] = p[C.kr161]
    p[C.kf165] = p[C.kf161]
    p[C.kr165] = p[C.kr161]
    p[C.kf166] = p[C.kf161]
    p[C.kr166] = p[C.kr161]
    p[C.kf167] = p[C.kf161]
    p[C.kr167] = p[C.kr161]
    p[C.kf168] = p[C.kf161]
    p[C.kr168] = p[C.kr161]
    p[C.kf169] = p[C.kf161]
    p[C.kr169] = p[C.kr161]
    p[C.kf170] = p[C.kf161]
    p[C.kr170] = p[C.kr161]
    p[C.kf171] = p[C.kf161]
    p[C.kr171] = p[C.kr161]
    p[C.kf172] = p[C.kf161]
    p[C.kr172] = p[C.kr161]
    p[C.kf173] = p[C.kf161]
    p[C.kr173] = p[C.kr161]
    p[C.kf175] = p[C.kf174]
    p[C.kr175] = p[C.kr174]
    p[C.kf176] = p[C.kf174]
    p[C.kr176] = p[C.kr174]
    p[C.kf177] = p[C.kf174]
    p[C.kr177] = p[C.kr174]
    p[C.kf178] = p[C.kf174]
    p[C.kr178] = p[C.kr174]
    p[C.kf179] = p[C.kf174]
    p[C.kr179] = p[C.kr174]
    p[C.kf180] = p[C.kf174]
    p[C.kr180] = p[C.kr174]
    p[C.kf181] = p[C.kf174]
    p[C.kr181] = p[C.kr174]
    p[C.kf182] = p[C.kf174]
    p[C.kr182] = p[C.kr174]
    p[C.kf183] = p[C.kf174]
    p[C.kr183] = p[C.kr174]
    p[C.kf184] = p[C.kf174]
    p[C.kr184] = p[C.kr174]
    p[C.kf185] = p[C.kf174]
    p[C.kr185] = p[C.kr174]
    p[C.kf186] = p[C.kf174]
    p[C.kr186] = p[C.kr174]
    p[C.kf187] = p[C.kf174]
    p[C.kr187] = p[C.kr174]
    p[C.kf188] = p[C.kf174]
    p[C.kr188] = p[C.kr174]
    p[C.kf189] = p[C.kf174]
    p[C.kr189] = p[C.kr174]
    p[C.kf190] = p[C.kf174]
    p[C.kr190] = p[C.kr174]
    p[C.kf191] = p[C.kf174]
    p[C.kr191] = p[C.kr174]
    p[C.kf192] = p[C.kf174]
    p[C.kr192] = p[C.kr174]
    p[C.kf194] = p[C.kf193]
    p[C.kr194] = p[C.kr193]
    p[C.kf195] = p[C.kf193]
    p[C.kr195] = p[C.kr193]
    p[C.kf196] = p[C.kf193]
    p[C.kr196] = p[C.kr193]
    p[C.kf197] = p[C.kf193]
    p[C.kr197] = p[C.kr193]
    p[C.kf198] = p[C.kf193]
    p[C.kr198] = p[C.kr193]
    p[C.kf199] = p[C.kf193]
    p[C.kr199] = p[C.kr193]
    p[C.kf200] = p[C.kf193]
    p[C.kr200] = p[C.kr193]
    p[C.kf201] = p[C.kf193]
    p[C.kr201] = p[C.kr193]
    p[C.kf202] = p[C.kf193]
    p[C.kr202] = p[C.kr193]
    p[C.kf203] = p[C.kf193]
    p[C.kr203] = p[C.kr193]
    p[C.kf204] = p[C.kf193]
    p[C.kr204] = p[C.kr193]
    p[C.kf205] = p[C.kf193]
    p[C.kr205] = p[C.kr193]
    p[C.kf206] = p[C.kf193]
    p[C.kr206] = p[C.kr193]
    p[C.kf207] = p[C.kf193]
    p[C.kr207] = p[C.kr193]
    p[C.kf208] = p[C.kf193]
    p[C.kr208] = p[C.kr193]
    p[C.kf209] = p[C.kf193]
    p[C.kr209] = p[C.kr193]
    p[C.kf210] = p[C.kf193]
    p[C.kr210] = p[C.kr193]
    p[C.kf211] = p[C.kf193]
    p[C.kr211] = p[C.kr193]
    p[C.kf213] = p[C.kf212]
    p[C.kr213] = p[C.kr212]
    p[C.kf214] = p[C.kf212]
    p[C.kr214] = p[C.kr212]
    p[C.kf215] = p[C.kf212]
    p[C.kr215] = p[C.kr212]
    p[C.kf216] = p[C.kf212]
    p[C.kr216] = p[C.kr212]
    p[C.kf217] = p[C.kf212]
    p[C.kr217] = p[C.kr212]
    p[C.kf218] = p[C.kf212]
    p[C.kr218] = p[C.kr212]
    p[C.kf219] = p[C.kf212]
    p[C.kr219] = p[C.kr212]
    p[C.kf220] = p[C.kf212]
    p[C.kr220] = p[C.kr212]
    p[C.kf221] = p[C.kf212]
    p[C.kr221] = p[C.kr212]
    p[C.kf222] = p[C.kf212]
    p[C.kr222] = p[C.kr212]
    p[C.kf223] = p[C.kf212]
    p[C.kr223] = p[C.kr212]
    p[C.kf224] = p[C.kf212]
    p[C.kr224] = p[C.kr212]
    p[C.kf226] = p[C.kf225]
    p[C.kr226] = p[C.kr225]
    p[C.kf227] = p[C.kf225]
    p[C.kr227] = p[C.kr225]
    p[C.kf228] = p[C.kf225]
    p[C.kr228] = p[C.kr225]
    p[C.kf229] = p[C.kf225]
    p[C.kr229] = p[C.kr225]
    p[C.kf230] = p[C.kf225]
    p[C.kr230] = p[C.kr225]
    p[C.kf231] = p[C.kf225]
    p[C.kr231] = p[C.kr225]
    p[C.kf232] = p[C.kf225]
    p[C.kr232] = p[C.kr225]
    p[C.kf233] = p[C.kf225]
    p[C.kr233] = p[C.kr225]
    p[C.kf234] = p[C.kf225]
    p[C.kr234] = p[C.kr225]
    p[C.kf235] = p[C.kf225]
    p[C.kr235] = p[C.kr225]
    p[C.kf236] = p[C.kf225]
    p[C.kr236] = p[C.kr225]
    p[C.kf237] = p[C.kf225]
    p[C.kr237] = p[C.kr225]
    p[C.kf238] = p[C.kf225]
    p[C.kr238] = p[C.kr225]
    p[C.kf239] = p[C.kf225]
    p[C.kr239] = p[C.kr225]
    p[C.kf240] = p[C.kf225]
    p[C.kr240] = p[C.kr225]
    p[C.kf241] = p[C.kf225]
    p[C.kr241] = p[C.kr225]
    p[C.kf243] = p[C.kf242]
    p[C.kr243] = p[C.kr242]
    p[C.kf244] = p[C.kf242]
    p[C.kr244] = p[C.kr242]
    p[C.kf245] = p[C.kf242]
    p[C.kr245] = p[C.kr242]
    p[C.kf246] = p[C.kf242]
    p[C.kr246] = p[C.kr242]
    p[C.kf247] = p[C.kf242]
    p[C.kr247] = p[C.kr242]
    p[C.kf248] = p[C.kf242]
    p[C.kr248] = p[C.kr242]
    p[C.kf249] = p[C.kf242]
    p[C.kr249] = p[C.kr242]
    p[C.kf250] = p[C.kf242]
    p[C.kr250] = p[C.kr242]
    p[C.kf251] = p[C.kf242]
    p[C.kr251] = p[C.kr242]
    p[C.kf252] = p[C.kf242]
    p[C.kr252] = p[C.kr242]
    p[C.kf253] = p[C.kf242]
    p[C.kr253] = p[C.kr242]
    p[C.kf254] = p[C.kf242]
    p[C.kr254] = p[C.kr242]
    p[C.kf255] = p[C.kf242]
    p[C.kr255] = p[C.kr242]
    p[C.kf256] = p[C.kf242]
    p[C.kr256] = p[C.kr242]
    p[C.kf257] = p[C.kf242]
    p[C.kr257] = p[C.kr242]
    p[C.kf258] = p[C.kf242]
    p[C.kr258] = p[C.kr242]
    p[C.kf269] = p[C.kf268]
    p[C.kf270] = p[C.kf268]
    p[C.kf271] = p[C.kf268]
    p[C.kf272] = p[C.kf268]
    p[C.kf273] = p[C.kf268]
    p[C.kf274] = p[C.kf268]
    p[C.kf275] = p[C.kf268]
    p[C.kf276] = p[C.kf268]
    p[C.kf277] = p[C.kf268]
    p[C.kf278] = p[C.kf268]
    p[C.kf279] = p[C.kf268]
    p[C.kf280] = p[C.kf268]
    p[C.kf281] = p[C.kf268]
    p[C.kf282] = p[C.kf268]
    p[C.kf283] = p[C.kf268]
    p[C.kf284] = p[C.kf268]
    p[C.kf285] = p[C.kf268]
    p[C.kf286] = p[C.kf268]
    p[C.kf287] = p[C.kf268]
    p[C.kf288] = p[C.kf268]
    p[C.kf289] = p[C.kf268]
    p[C.kf290] = p[C.kf268]
    p[C.kf303] = p[C.kf298]
    p[C.kr303] = p[C.kr298]
    p[C.kf304] = p[C.kf299]
    p[C.kr304] = p[C.kr299]
    p[C.kf305] = p[C.kf300]
    p[C.kr305] = p[C.kr300]
    p[C.kf306] = p[C.kf301]
    p[C.kr306] = p[C.kr301]
    p[C.kf307] = p[C.kf302]
    p[C.kr307] = p[C.kr302]
    # --------------------------------------------------------------------------

    return p, u0
end

function decode_gene2val(indiv_gene::Vector{Float64})::Vector{Float64}
    search_rgn::Matrix{Float64} = get_search_region()
    indiv::Vector{Float64} = zeros(length(indiv_gene))

    for i in eachindex(indiv_gene)
        indiv[i] = 10^(
            indiv_gene[i] * (
                search_rgn[2,i] - search_rgn[1,i]
            ) + search_rgn[1,i]
        )
    end

    return round.(indiv, sigdigits=7)
end

function encode_val2gene(indiv::Vector{Float64})
    search_rgn::Matrix{Float64} = get_search_region()
    indiv_gene::Vector{Float64} = zeros(length(indiv))

    for i in eachindex(indiv)
        indiv_gene[i] = (
            log10(indiv[i]) - search_rgn[1,i]
        ) / (
            search_rgn[2,i] - search_rgn[1,i]
        )
    end

    return indiv_gene
end

function encode_bestIndivVal2randGene(
        idx::Int64,
        best_indiv::Vector{Float64},
        p0_bounds::Vector{Float64})::Float64
    search_rgn::Matrix{Float64} = get_search_region()
    rand_gene::Float64 = (
        log10(
            best_indiv[idx] * 10^(
                rand() * log10(p0_bounds[2] / p0_bounds[1]) + log10(p0_bounds[1])
            )
        ) - search_rgn[1,idx]
    ) / (
        search_rgn[2,idx] - search_rgn[1,idx]
    )
    return rand_gene
end

function init_search_param(
        search_idx::Tuple{Array{Int64,1},Array{Int64,1}},
        p::Vector{Float64},
        u0::Vector{Float64})::Vector{Float64}
    duplicate::Vector{String} = []
    if length(search_idx[1]) != length(unique(search_idx[1]))
        for idx in findall([count(x -> x == i, search_idx[1])
                            for i in unique(search_idx[1])] .!= 1)
            push!(duplicate, C.NAMES[search_idx[1][idx]])
        end
        error(
            "Duplicate parameters (C.): $duplicate"
        )
    elseif length(search_idx[2]) != length(unique(search_idx[2]))
        for idx in findall([count(x -> x == i, search_idx[2])
                            for i in unique(search_idx[2])] .!= 1)
            push!(duplicate, V.NAMES[search_idx[2][idx]])
        end
        error(
            "Duplicate species (V.): $duplicate"
        )
    end
    search_param = zeros(
        length(search_idx[1]) + length(search_idx[2])
    )
    for (i, j) in enumerate(search_idx[1])
        @inbounds search_param[i] = p[j]
    end
    for (i, j) in enumerate(search_idx[2])
        @inbounds search_param[i + length(search_idx[1])] = u0[j]
    end

    if any(x -> x == 0.0, search_param)
        msg::String = "search_param must not contain zero."
        for idx in search_idx[1]
            if p[idx] == 0.0
                error(
                    @sprintf(
                        "`C.%s` in search_idx_params: ", C.NAMES[idx]
                    ) * msg
                )
            end
        end
        for idx in search_idx[2]
            if u0[idx] == 0.0
                error(
                    @sprintf(
                        "`V.%s` in search_idx_initials: ", V.NAMES[idx]
                    ) * msg
                )
            end
        end
    end

    return search_param
end

function conv_lin2log!(
        search_rgn::Matrix{Float64},
        search_idx::Tuple{Array{Int64,1},Array{Int64,1}})::Matrix{Float64}
    for i = 1:size(search_rgn, 2)
        if minimum(search_rgn[:,i]) < 0.0
            msg = "search_rgn[lower_bound,upper_bound] must be positive.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
            end
        elseif minimum(search_rgn[:,i]) == 0.0 && maximum(search_rgn[:,i]) != 0.0
            msg = "lower_bound must be larger than 0.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
            end
        elseif search_rgn[2,i] - search_rgn[1,i] < 0.0
            msg = "lower_bound must be smaller than upper_bound.\n"
            if i <= C.NUM
                error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
            else
                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
            end
        end
    end

    nonzero_idx::Vector{Int} = []
    for i = 1:size(search_rgn, 2)
        if search_rgn[:,i] != [0.0,0.0]
            push!(nonzero_idx, i)
        end
    end
    difference::Vector{Int} = collect(
        symdiff(
            Set(nonzero_idx),
            Set(append!(search_idx[1], C.NUM .+ search_idx[2]))
        )
    )
    if length(difference) > 0
        for idx in difference
            if idx <= C.NUM
                println(@sprintf("`C.%s`", C.NAMES[Int(idx)]))
            else
                println(@sprintf("`V.%s`", V.NAMES[Int(idx) - C.NUM]))
            end
        end
        error(
            "Set these search_params in both search_idx and search_rgn."
        )
    end

    search_rgn = search_rgn[:,nonzero_idx]

    return log10.(search_rgn)
end

ext: jl


text:
### TODO: CLEAN UP THIS FUNCTION AND DOCUMENTATION. THERE SHOULD BE NO NEED TO
### HAVE A TYPE-UNSTABLE FUNCTION HERE. (sbromberger 2017-03-26)
"""
    multiroute_flow(flow_graph, source, target[, DefaultCapacity][, flow_algorithm][, mrf_algorithm][, routes])

The generic multiroute_flow function.

The output will vary depending on the input:

- When the number of `route`s is `0`, return the set of breaking points of
the multiroute flow.
- When the number of `route`s is `1`, return a flow with a set of 1-disjoint paths
(this is the classical max-flow implementation).
- When the input is limited to a set of breaking points and a route value `k`,
return only the k-route flow.
- Otherwise, a tuple with 1) the maximum flow and 2) the flow matrix. When the
max-flow subroutine is the Boykov-Kolmogorov algorithm, the associated mincut is
returned as a third output.

When the input is a network, it requires the following arguments:

- `flow_graph`: the input graph
- `source`: the source vertex
- `target`: the target vertex
- `capacity_matrix`: matrix of edge flow capacities
- `flow_algorithm`: keyword argument for flow algorithm
- `mrf_algorithm`: keyword argument for multiroute flow algorithm
- `routes`: keyword argument for the number of routes

When the input is only the set of (breaking) points and the number of route,
it requires the following arguments:

- `breakingpoints`: vector of breaking points
- `routes`: number of routes

When the input is the set of (breaking) points, the number of routes,
and the network descriptors, it requires the following arguments:

- `breakingpoints`: vector of breaking points
- `routes`: number of routes
- `flow_graph`: the input graph
- `source`: the source vertex
- `target`: the target vertex
- `capacity_matrix`: matrix of edge flow capacities
- `flow_algorithm`: keyword argument for flow algorithm

The function defaults to the Push-relabel (classical flow) and Kishimoto
(multiroute) algorithms. Alternatively, the algorithms to be used can also
be specified through  keyword arguments. A default capacity of `1` is assumed
for each link if no capacity matrix is provided.

The `mrf_algorithm` keyword is inforced to Extended Multiroute Flow
in the following cases:

- The number of routes is non-integer
- The number of routes is 0 or non-specified

### Usage Example :
(please consult the  [`maximum_flow`](@ref) section for options about flow_algorithm
and capacity_matrix)

```jldoctest
julia> flow_graph = DiGraph(8) # Create a flow graph

julia> flow_edges = [
(1, 2, 10), (1, 3, 5),  (1, 4, 15), (2, 3, 4),  (2, 5, 9),
(2, 6, 15), (3, 4, 4),  (3, 6, 8),  (4, 7, 16), (5, 6, 15),
(5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6),  (7, 8, 10)
]

julia> capacity_matrix = zeros(Int, 8, 8) # Create a capacity matrix

julia> for e in flow_edges
    u, v, f = e
    add_edge!(flow_graph, u, v)
    capacity_matrix[u, v] = f
end

julia> f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 2) # Run default multiroute_flow with an integer number of routes = 2

julia> f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 1.5) # Run default multiroute_flow with a noninteger number of routes = 1.5

julia> points = multiroute_flow(flow_graph, 1, 8, capacity_matrix) # Run default multiroute_flow for all the breaking points values

julia> f, F = multiroute_flow(points, 1.5) # Then run multiroute flow algorithm for any positive number of routes

julia> f = multiroute_flow(points, 1.5, valueonly = true)

julia> f, F, labels = multiroute_flow(flow_graph, 1, 8, capacity_matrix, algorithm = BoykovKolmogorovAlgorithm(), routes = 2) # Run multiroute flow algorithm using Boykov-Kolmogorov algorithm as maximum_flow routine

```
"""
function multiroute_flow(
    flow_graph::AbstractGraph,                    # the input graph
    source::Integer,                        # the source vertex
    target::Integer,                        # the target vertex
    capacity_matrix::AbstractMatrix{T} =  # edge flow capacities
        DefaultCapacity(flow_graph);
    flow_algorithm::AbstractFlowAlgorithm  =    # keyword argument for algorithm
    PushRelabelAlgorithm(),
    mrf_algorithm::AbstractMultirouteFlowAlgorithm  =    # keyword argument for algorithm
    KishimotoAlgorithm(),
    routes::R = 0              # keyword argument for number of routes (0 = all values)
    ) where T where R <: Real

    # a flow with a set of 1-disjoint paths is a classical max-flow
    (routes == 1) &&
    return maximum_flow(flow_graph, source, target, capacity_matrix, flow_algorithm)

    # routes > λ (connectivity) → f = 0
    λ = maximum_flow(flow_graph, source, target, DefaultCapacity(flow_graph),
    algorithm = flow_algorithm)[1]
    (routes > λ) && return empty_flow(capacity_matrix, flow_algorithm)

    # For other cases, capacities need to be Floats
    if !(T<:AbstractFloat)
        capacity_matrix = convert(AbstractMatrix{Float64}, capacity_matrix)
    end

    # Ask for all possible values (breaking points)
    (routes == 0) &&
    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm)
    # The number of routes is a float → EMRF
    (R <: AbstractFloat) &&
    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)

    # Other calls
    return multiroute_flow(flow_graph, source, target, capacity_matrix,
    flow_algorithm, mrf_algorithm, routes)
end


text:
@defcomp climaten2ocycle begin
    # Global N2O emissions in Mt of N
    globn2o = Parameter(index=[time])

    # Atmospheric N2O concentration
    acn2o = Variable(index=[time])

    # N2O decay
    n2odecay = Variable()

    #
    lifen2o = Parameter(default = 114.0)

    # N2o pre industrial
    n2opre = Parameter(default = 285.0)

    function run_timestep(p, v, d, t)

        if is_first(t)
            v.n2odecay = 1.0 / p.lifen2o

            v.acn2o[t] = 296
        else
            # Calculate N2O concentrations
            v.acn2o[t] = v.acn2o[t - 1] + 0.2079 * p.globn2o[t] - v.n2odecay * (v.acn2o[t - 1] - p.n2opre)

            if v.acn2o[t] < 0
                error("n2o atmospheric concentration out of range")
            end
        end
    end
end


text:
# includes helper functions for the parallel solver

segment(n::Int64, u) = error("$(typeof(u)) does not implement segment")

function segment(n::Int64, u::UnitRange)
    chunks = {}

    cov = u.stop - u.start
    stride = int(floor((u.stop - u.start) / (n-1)))
    for i = 1:(n-1)
        si = u.start + (i-1)*stride
        ei = si + stride - 1
        i == (n-1) ? (ei = u.stop) : nothing
        push!(chunks, si:ei)
    end
    return chunks
end


text:
using LinearAlgebra
using SparseArrays
using Random
using Test

Random.seed!(0)

# Simply test the call to functions and their prox

test_cases_spec = [
    Dict(
        "constr" => DistL2,
        "wrong"  => [
            (IndBallL2(), -rand()),
            (IndBinary(),)
        ],
        "right" => [
            ( (IndSOC(),), randn(Float32, 10) ),
            ( (IndSOC(),), randn(Float64, 10) ),
            ( (IndNonnegative(), rand()),  randn(Float64, 10) ),
            ( (IndZero(),), randn(Float64, 10) ),
            ( (IndBox(-1, 1),), randn(Float32, 10) ),
            ( (IndBox(-1, 1),), randn(Float64, 10) ),
        ],
    ),

    Dict(
        "constr" => SqrDistL2,
        "wrong"  => [
            (IndBallL2(), -rand()),
            (IndBinary(),)
        ],
        "right" => [
            ( (IndSimplex(),), randn(Float32, 10) ),
            ( (IndSimplex(),), randn(Float64, 10) ),
            ( (IndNonnegative(), rand()), randn(Float64, 10) ),
            ( (IndZero(),), randn(Float64, 10) ),
            ( (IndBox(-1, 1),), randn(Float32, 10) ),
            ( (IndBox(-1, 1),), randn(Float64, 10) ),
        ],
    ),

    Dict(
        "constr" => ElasticNet,
        "wrong"  => [
            (-rand()),
            (-rand(), -rand()),
            (-rand(), rand()),
            (rand(), -rand())
        ],
        "right" => [
            ( (), randn(Float64, 10) ),
            ( (rand(Float32),), randn(Float32, 10) ),
            ( (rand(Float64), rand(Float64)), randn(Float64, 10) ),
            ( (rand(Float64), rand(Float64)), rand(Complex{Float64}, 20) ),
        ],
    ),

    Dict(
        "constr" => HingeLoss,
        "wrong"  => [
            (randn(10), -rand()),
        ],
        "right" => [
            ( (Int.(sign.(randn(10))), ), randn(Float32, 10) ),
            ( (Int.(sign.(randn(10))), ), randn(Float64, 10) ),
            ( (Int.(sign.(randn(20))), 0.1f0 + rand(Float32)), randn(Float32, 20) ),
            ( (Int.(sign.(randn(20))), 0.1e0 + rand(Float64)), randn(Float64, 20) ),
        ],
    ),

    Dict(
        "constr" => IndAffine,
        "right" => [
            ( (randn(10), randn()), randn(Float32, 10) ),
            ( (randn(10, 20), randn(10)), randn(20) ),
            ( (sprand(50,100, 0.1), randn(50)), randn(Float32, 100) ),
            ( (sprand(Complex{Float64}, 50, 100, 0.1), randn(50)+im*randn(50)), randn(100)+im*randn(100) ),
        ],
    ),

    Dict(
        "constr" => IndBallLinf,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (rand(Float32),), randn(Float32, 10) ),
            ( (rand(Float64),), randn(Float64, 10) ),
        ],
    ),

    Dict(
        "constr" => IndBallL0,
        "wrong"  => [
            (-4,),
        ],
        "right" => [
            ( (5,), randn(25) ),
            ( (10, ), randn(Float32, 5, 5) ),
            ( (5, ), randn(Complex{Float64}, 15) ),
        ],
    ),

    Dict(
        "constr" => IndBallL1,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (), [0.1, -0.2, 0.3, -0.39] ),
            ( (), Complex{Float64}[0.1, -0.2, 0.3, -0.39] ),
            ( (), rand(Float32, 5) ),
            ( (), rand(Float64, 5) ),
            ( (), rand(Complex{Float32}, 5) ),
            ( (Float32(3),), randn(Float32, 5) ),
            ( (Float64(5),), randn(Float64, 2, 3) ),
            ( (Float64(1),), randn(Complex{Float64}, 5) ),
        ],
    ),

    Dict(
        "constr" => IndBallL2,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (), rand(Float32, 5) ),
            ( (), rand(Float64, 5) ),
            ( (), rand(Complex{Float32}, 5) ),
            ( (Float32(3),), randn(Float32, 5) ),
            ( (Float64(5),), randn(Float64, 2, 3) ),
            ( (Float64(1),), randn(Complex{Float64}, 5) ),
        ],
    ),

    Dict(
        "constr" => IndBallRank,
        "wrong"  => [
            (-2,),
        ],
        "right" => [
            ( (1+Int(round(10*rand())),), randn(20, 50) ),
            ( (10+Int(round(5*rand())),), rand(30, 8)*rand(8,70) ),
            ( (10+Int(round(5*rand())),), randn(Float32, 5, 8) ),
            ( (10+Int(round(5*rand())),), rand(Complex{Float64}, 20, 50) ),
            ( (10+Int(round(5*rand())),), rand(Complex{Float32}, 5, 8) ),
        ],
    ),

    Dict(
        "constr" => IndBox,
        "wrong"  => [
            (+1, -1),
        ],
        "right" => [
            ( (-rand(Float32), +rand(Float32, 10)), randn(Float32, 10) ),
            ( (-rand(Float32, 10), +rand(Float32)), randn(Float32, 10) ),
            ( (-rand(Float64), +rand(Float64)), randn(Float64, 20) ),
            ( (-rand(Float64, 30), +rand(Float64, 30)), randn(Float64, 30) ),
        ],
    ),

    Dict(
        "constr" => IndExpPrimal,
        "right" => [
            ( (), randn(Float32, 3) ),
            ( (), randn(Float64, 3) ),
        ],
    ),

    Dict(
        "constr" => IndExpDual,
        "right" => [
            ( (), randn(Float32, 3) ),
            ( (), randn(Float64, 3) ),
        ],
    ),

    Dict(
        "constr" => IndFree,
        "right" => [
            ( (), randn(Float32, 5) ),
            ( (), randn(Float64, 2, 3) ),
        ],
    ),

    Dict(
        "constr" => IndGraph,
        "right" => [
            ( (sprand(50, 100, 0.2),), randn(150) ),
            ( (sprand(Complex{Float64}, 50, 100, 0.2),), randn(150)+im*randn(150) ),
            ( (rand(50, 100),), randn(150) ),
            ( (rand(Complex{Float64}, 50, 100),), randn(150)+im*randn(150) ),
            ( (rand(55, 50),), randn(105) ),
            ( (rand(Complex{Float64}, 55, 50),), randn(105)+im*randn(105) ),
        ],
    ),

    Dict(
        "constr" => IndPoint,
        "right" => [
            ( (), randn(5) ),
            ( (randn(Float32, 5),), randn(Float32, 5) ),
            ( (randn(Float64, 5),), randn(Float32, 5) ),
            ( (randn(Float64, 5),), randn(Float64, 5) ),
            ( (randn(Complex{Float32}, 5),), randn(Complex{Float32}, 5) ),
            ( (randn(Complex{Float64}, 5),), randn(Complex{Float32}, 5) ),
            ( (randn(Complex{Float64}, 5),), randn(Complex{Float64}, 5) ),
        ],
    ),

    Dict(
        "constr" => IndStiefel,
        "right" => [
            ( (), randn(Float32, 5, 3) ),
            ( (), randn(Float64, 5, 3) ),
            ( (), randn(Complex{Float32}, 5, 3) ),
            ( (), randn(Complex{Float64}, 5, 3) ),
        ],
    ),

    Dict(
        "constr" => IndZero,
        "right" => [
            ( (), randn(Float32, 5) ),
            ( (), randn(Float32, 2, 3) ),
            ( (), randn(Float64, 5) ),
            ( (), randn(Float64, 2, 3) ),
            ( (), randn(Complex{Float32}, 5) ),
            ( (), randn(Complex{Float32}, 2, 3) ),
        ],
    ),

    Dict(
        "constr" => IndHalfspace,
        "right" => [
            ( (-ones(5), -1.0), -rand(5) ),
            ( (-ones(5), 1.0), rand(5) ),
        ],
    ),

    Dict(
        "constr" => IndNonnegative,
        "right" => [
            ( (), randn(10) ),
            ( (), randn(3, 5) ),
        ],
    ),

    Dict(
        "constr" => IndNonpositive,
        "right" => [
            ( (), randn(10) ),
            ( (), randn(3, 5) ),
        ],
    ),

    Dict(
        "constr" => IndSimplex,
        "wrong"  => [
            (-rand(),)
        ],
        "right" => [
            ( (), randn(10) ),
            ( (), randn(3, 5) ),
            ( (5.0,), randn(10) ),
        ],
    ),

    Dict(
        "constr" => IndSOC,
        "right" => [
            ( (), [rand(), -rand()] ),
            ( (), [-rand(), rand()] ),
            ( (), rand(5) ),
        ],
    ),

    Dict(
        "constr" => IndRotatedSOC,
        "right" => [
            ( (), [rand(), -rand()] ),
            ( (), [-rand(), rand()] ),
            ( (), rand(5) ),
        ],
    ),

    Dict(
        "constr" => IndSphereL2,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (rand(),), randn(10) ),
            ( (sqrt(20),), randn(20) ),
            ( (1,), 1e-3*randn(10) ),
            ( (1,), randn(Complex{Float64}, 10) ),
        ],
    ),

    Dict(
        "constr" => IndPSD,
        "right" => [
            ( (), Symmetric(randn(Float32, 5, 5)) ),
            ( (), Symmetric(randn(Float64, 10, 10)) ),
            ( (), Hermitian(randn(Complex{Float32}, 5, 5)) ),
            ( (), Hermitian(randn(Complex{Float64}, 10, 10)) ),
            ( (), randn(Float32, 5, 5) ),
            ( (), randn(Float64, 10, 10) ),
            ( (), randn(Complex{Float32}, 5, 5) ),
            ( (), randn(Complex{Float64}, 10, 10) ),
            ( (), randn(15) ),
        ],
    ),

    Dict(
        "constr" => LogBarrier,
        "wrong"  => [
            (1.0, 0.0, -rand()),
        ],
        "right" => [
            ( (), rand(Float32, 5) ),
            ( (), rand(Float64, 5) ),
            ( (rand(Float32),), rand(Float32, 5) ),
            ( (rand(Float64), rand(Float64)), rand(Float64, 5) ),
            ( (rand(Float32), rand(Float32), rand(Float32)), rand(Float32, 5) ),
        ],
    ),

    Dict(
        "constr" => NormL0,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (), randn(Float32, 5) ),
            ( (), randn(Float64, 5) ),
            ( (rand(Float32)), randn(Complex{Float32}, 5) ),
        ],
    ),

    Dict(
        "constr" => NormL1,
        "wrong"  => [
            (-rand(),),
            (-rand(10),)
        ],
        "right" => [
            ( (), randn(Float32, 5) ),
            ( (), randn(Float64, 5) ),
            ( (), randn(Complex{Float32}, 5) ),
            ( (rand(Float32),), randn(Float32, 5) ),
            ( (rand(Float64, 5),), randn(Float64, 5) ),
            ( (rand(Float64, 5),), rand(Complex{Float64}, 5) ),
            ( (rand(Float32),), rand(Complex{Float32}, 5) ),
        ],
    ),

    Dict(
        "constr" => NormL2,
        "right" => [
            ( (), randn(Float32, 5), ),
            ( (), randn(Float64, 5), ),
            ( (), randn(Complex{Float32}, 5), ),
            ( (rand(Float32),), randn(Float32, 5), ),
            ( (rand(Float64),), randn(Complex{Float64}, 5), ),
        ],
    ),

    Dict(
        "constr" => NormL21,
        "right" => [
            ( (), randn(Float32, 3, 5) ),
            ( (), randn(Float64, 3, 5) ),
            ( (), randn(Complex{Float32}, 3, 5) ),
            ( (rand(Float32),), randn(Float32, 3, 5) ),
            ( (rand(Float64),), randn(Float64, 3, 5) ),
            ( (rand(Float32), 1), randn(Complex{Float32}, 3, 5) ),
            ( (rand(Float64), 2), randn(Complex{Float64}, 3, 5) ),
        ],
    ),

    Dict(
        "constr" => NormLinf,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (), randn(Float32, 5) ),
            ( (), randn(Float64, 5) ),
            ( (), randn(Complex{Float32}, 5) ),
            ( (rand(Float32),), randn(Float32, 5) ),
            ( (rand(Float32),), rand(Complex{Float32}, 5) ),
        ],
    ),

    Dict(
        "constr" => NuclearNorm,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (), rand(Float32, 5, 5) ),
            ( (), rand(Float64, 3, 5) ),
            ( (), rand(Complex{Float32}, 5, 5) ),
            ( (rand(Float32),), rand(Float32, 5, 5) ),
            ( (rand(Float64),), rand(Float64, 3, 5) ),
            ( (rand(Float32),), rand(Complex{Float32}, 5, 5) ),
            ( (sqrt(10e0),), rand(Complex{Float64}, 3, 5) ),
            ( (sqrt(10f0),), rand(Complex{Float32}, 5, 3) ),
        ],
    ),

    Dict(
        "constr" => SqrNormL2,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (), randn(Float32, 10) ),
            ( (), randn(Float64, 10) ),
            ( (rand(Float32),), randn(Float32, 10) ),
            ( (rand(Float64),), randn(Float64, 10) ),
            ( (rand(Float32, 20),), randn(Float32, 20) ),
            ( (rand(Float64, 20),), randn(Float64, 20) ),
            ( (rand(30),), rand(Complex{Float64}, 30) ),
            ( (rand(),), rand(Complex{Float64}, 50) ),
        ],
    ),

    Dict(
        "constr" => LeastSquares,
        "wrong"  => [
            (randn(3, 5), randn(4), rand()),
        ],
        "right" => [
            ( (randn(10, 25), randn(10)), randn(25) ),
            ( (randn(40, 13), randn(40), rand()), randn(13) ),
            ( (rand(Complex{Float64}, 25, 10), rand(Complex{Float64}, 25)), rand(Complex{Float64}, 10) ),
            ( (sprandn(10, 100, 0.15), randn(10), rand()), randn(100) ),
        ],
    ),

    # Dict(
    #     "constr" => SumLargest,
    #     "wrong"  => [
    #         (-1,), (0,), (1, -2.0), (2, 0.0), (2, -rand())
    #     ],
    #     "right" => [
    #         ( (), randn(10) ),
    #         ( (1,), randn(20) ),
    #         ( (5, 3.2), randn(5,10) ),
    #         ( (3, rand()), randn(8,17) ),
    #     ],
    # ),

    Dict(
        "constr" => Maximum,
        "wrong"  => [
            (-rand(),),
        ],
        "right" => [
            ( (), randn(Float32, 10) ),
            ( (), randn(Float64, 10) ),
            ( (rand(Float64),), randn(Float64, 20) ),
            ( (), randn(Float32, 5,10) ),
            ( (rand(Float32)), randn(Float32, 8,17) ),
        ],
    ),

    Dict(
        "constr" => IndBinary,
        "right" => [
            ( (), randn(10) ),
            ( (randn(), randn()), randn(10) ),
            ( (randn(), randn(10)), randn(10) ),
            ( (randn(10), randn()), randn(10) ),
            ( (randn(10), randn(10)), randn(10) ),
            ( (randn(), randn(5,5)), randn(5,5) ),
        ],
    ),

    Dict(
        "constr" => Regularize,
        "wrong"  => [
            (NormL1(), -rand()),
        ],
        "right" => [
            ( (NormL1(), rand()), randn(10) ),
            ( (NormL1(rand()), rand()), randn(5,10) ),
            ( (NormL1(), rand(), randn(20)), randn(20) ),
            ( (NormL1(rand()), rand(), randn(20)), randn(20) ),
        ],
    ),

    Dict(
        "constr" => Tilt,
        "right" => [
            ( (LeastSquares(randn(20, 10), randn(20)), randn(10)), randn(10) ),
        ],
    ),

    Dict(
        "constr" => HuberLoss,
        "wrong"  => [
            (-rand(), ),
            (rand(), -rand()),
            (-rand(), rand()),
            (-rand(), -rand())
        ],
        "right" => [
            ( (), randn(Float32, 10) ),
            ( (), randn(Float64, 10) ),
            ( (rand(Float32), ), randn(Float32, 5, 8) ),
            ( (rand(Float64), ), randn(Float64, 5, 8) ),
            ( (rand(Float64), rand(Float64)), randn(Float64, 20) ),
            ( (rand(Float64), rand(Float64)), rand(Complex{Float64}, 8, 12) ),
        ],
    ),

    Dict(
        "constr" => SumPositive,
        "right" => [
            ( (), randn(Float32, 5) ),
            ( (), randn(Float64, 5) ),
            ( (), randn(Float32, 3, 5) ),
            ( (), randn(Float64, 3, 5) ),
        ],
    ),

    Dict(
        "constr" => SeparableSum,
        "right" => [
            ( ((NormL2(2.0), NormL1(1.5), NormL2(0.5)), ), (randn(5), randn(15), randn(10)) ),
        ],
    ),

    Dict(
        "constr" => SlicedSeparableSum,
        "right" => [
            ( ((NormL2(2.0), NormL1(1.5), NormL2(0.5)), ((1:5,), (6:20,), (21:30,))), randn(30) ),
        ],
    ),
]

@testset "$(spec["constr"])" for spec in test_cases_spec
    constr = spec["constr"]

    if haskey(spec, "wrong")
        for wrong in spec["wrong"]
            @test_throws ErrorException constr(wrong...)
        end
    end

    for right in spec["right"]
        params, x = right
        f = constr(params...)

        predicates_test(f)

        ##### just call f
        fx = call_test(f, x)

        ##### compute prox with default gamma
        y, fy = prox_test(f, x)

        ##### compute prox with random gamma
        T = if typeof(x) <: Array
            eltype(x)
        elseif typeof(x) <: Tuple
            eltype(x[1])
        else
            Float64
        end
        gam = 5*rand(real(T))
        y, fy = prox_test(f, x, gam)

        ##### compute prox with multiple random gammas
        if ProximalOperators.is_separable(f)
            gam = real(T)(0.5) .+ 2 .* rand(real(T), size(x))
            y, fy = prox_test(f, x, gam)
        end

    end
end


text:
"""
    pearson_bounds(dA::UnivariateDistribution, dB::UnivariateDistribution; n::Int=7)

Compute the theoretical lower and upper Pearson correlation values for a pair of
univariate distributions.

See also: [`pearson_match`](@ref)

# Examples
```jldoctest
julia> using Distributions

julia> A = Normal(78, 10); B = LogNormal(3, 1);

julia> pearson_bounds(A, B)
(lower = -0.7628739783668452, upper = 0.762873978367046)
```
"""
function pearson_bounds(dA::UD, dB::UD; n::Int=7)
    μA = mean(dA)
    σA = std(dA)
    μB = mean(dB)
    σB = std(dB)
    _pearson_bounds(dA, dB, μA, μB, σA, σB, n)
end
function _pearson_bounds(dA::UD, dB::UD, μA, μB, σA, σB, n)
    k = 0:1:n
    a = get_coefs(dA, n)
    b = get_coefs(dB, n)

    c1 = -μA * μB
    c2 = 1 / (σA * σB)
    kab = factorial.(k) .* a .* b
    ρ_l = c1 * c2 + c2 * sum((-1) .^ k .* kab)
    ρ_u = c1 * c2 + c2 * sum(kab)

    ρ_l, ρ_u = clampcor.((ρ_l, ρ_u))
    (lower = ρ_l, upper = ρ_u)
end

"""
    pearson_bounds(margins::Vector{<:UnivariateDistribution})

Compute the theoretical lower and upper Pearson correlation values for a set of
univariate distributions.

See also: [`pearson_match`](@ref)

# Examples
```jldoctest
julia> using Distributions

julia> m = [Normal(78, 10), LogNormal(3, 1)];

julia> b = pearson_bounds(m);

julia> b.lower
2×2 Array{Float64,2}:
  1.0       -0.762874
 -0.762874   1.0

julia> b.upper
2×2 Array{Float64,2}:
 1.0       0.762874
 0.762874  1.0
```
"""
function pearson_bounds(margins::Vector{<:UD})
    d = length(margins)
    lower, upper = zeros(Float64, d, d), zeros(Float64, d, d)

    @threads for i in collect(subsets(1:d, Val{2}()))
        l, u = pearson_bounds(margins[i[1]], margins[i[2]])
        lower[i...] = l
        upper[i...] = u
    end

    lower .= Symmetric(lower)
    cor_constrain!(lower)

    upper .= Symmetric(upper)
    cor_constrain!(upper)

    (lower = lower, upper = upper)
end


text:
using CodecBase
using BenchmarkTools

# Make data.
data_bin = read(joinpath(@__DIR__, "data.bin"))
data16 = transcode(Base16Encoder(), data_bin)
data32 = transcode(Base32Encoder(), data_bin)
data64 = transcode(Base64Encoder(), data_bin)

# Define benchmarks.
suite = BenchmarkGroup()
suite["Base16"] = BenchmarkGroup()
suite["Base16"]["Encoder"] = @benchmarkable transcode(Base16Encoder(), data_bin)
suite["Base16"]["Decoder"] = @benchmarkable transcode(Base16Decoder(), data16)
suite["Base32"] = BenchmarkGroup()
suite["Base32"]["Encoder"] = @benchmarkable transcode(Base32Encoder(), data_bin)
suite["Base32"]["Decoder"] = @benchmarkable transcode(Base32Decoder(), data32)
suite["Base64"] = BenchmarkGroup()
suite["Base64"]["Encoder"] = @benchmarkable transcode(Base64Encoder(), data_bin)
suite["Base64"]["Decoder"] = @benchmarkable transcode(Base64Decoder(), data64)

# Run benchmarks and print the results.
tune!(suite)
results = run(suite, verbose=true, seconds=1)
for k1 in ["Base16", "Base32", "Base64"]
    for k2 in ["Encoder", "Decoder"]
        println("--- $(k1) $(k2) ---")
        show(STDOUT, MIME"text/plain"(), results[k1][k2])
        println()
        println()
    end
end


text:
#=
module Encodings

export Encoding
export Binary, ASCII, Latin1, UTF8      # 7/8-bit encodings
export UCS2, UCS2LE, UCS2BE, UCS2OE     # 16-bit encodings (16-bit subset of Unicode)
export UTF16, UTF16LE, UTF16BE, UTF16OE # 16-bit encodings
export UTF32, UTF32LE, UTF32BE, UTF32OE # 32-bit encodings
export BIG_ENDIAN
export native_endian, big_endian, codeunit
=#

abstract Encoding
abstract DirectIndexedEncoding <: Encoding

abstract Binary <: DirectIndexedEncoding
abstract ASCII  <: DirectIndexedEncoding
abstract Latin1 <: DirectIndexedEncoding

abstract UTF8   <: Encoding
abstract UTF16  <: Encoding
abstract UTF32  <: DirectIndexedEncoding
abstract UCS2   <: DirectIndexedEncoding

# Opposite endian encodings of 16-bit and 32-bit encodings
abstract UTF16OE <: UTF16
abstract UTF32OE <: UTF32
abstract UCS2OE  <: UCS2

# This is easier to use (and not get the ordering mixed up) than ENDIAN_BOM
const BIG_ENDIAN = reinterpret(UInt32,UInt8[1:4;])[1] == 0x01020304

if BIG_ENDIAN
    abstract UTF16BE <: UTF16
    abstract UTF32BE <: UTF32
    abstract UCS2BE  <: UCS2
    abstract UTF16LE <: UTF16OE
    abstract UTF32LE <: UTF32OE
    abstract UCS2LE  <: UCS2OE
else
    abstract UTF16LE <: UTF16
    abstract UTF32LE <: UTF32
    abstract UCS2LE  <: UCS2
    abstract UTF16BE <: UTF16OE
    abstract UTF32BE <: UTF32OE
    abstract UCS2BE  <: UCS2OE
end

native_endian{E <: Encoding}(::Type{E}) = true
native_endian{E <: UTF16OE}(::Type{E})  = false
native_endian{E <: UTF32OE}(::Type{E})  = false
native_endian{E <: UCS2OE}(::Type{E})   = false

if BIG_ENDIAN
big_endian{E <: Encoding}(::Type{E}) = native_endian(E)
else
big_endian{E <: Encoding}(::Type{E}) = !native_endian(E)
end

codeunit{E <: ASCII}(::Type{E})  = UInt8
codeunit{E <: Latin1}(::Type{E}) = UInt8
codeunit{E <: UTF8}(::Type{E})   = UInt8
codeunit{E <: UTF16}(::Type{E})  = UInt16
codeunit{E <: UCS2}(::Type{E})   = UInt16
codeunit{E <: UTF32}(::Type{E})  = UInt32

# size of code unit in bytes
Base.sizeof{E <: Encoding}(::Type{E}) = sizeof(codeunit(E))

#end


text:
mutable struct TileStates
    quads::Array{Coord, 2}
    chars::Array{Char, 2}
    rands::Array{UInt8, 2}

    TileStates() = new(Matrix{Coord}(undef, 0, 0), Matrix{Char}(undef, 0, 0), Matrix{Int}(undef, 0, 0))
end

mutable struct DrawableRoom
    map::Map
    room::Room

    fgTileStates::TileStates
    bgTileStates::TileStates

    rendering::Union{Layer, Nothing}
    layers::Array{Layer, 1}

    fillColor::Union{Nothing, colorTupleType}
end

Base.size(state::TileStates) = size(state.rands)

getTileStateSeed(name::String, package::String="", fg::Bool=false) = foldl((a, b) -> a + Int(b) * 128, collect(package * (fg ? "fg" : "bg") * name), init=0)
getTileStateSeed(room::Room, package::String="", fg::Bool=false) = getTileStateSeed(room.name, package, fg)

const drawableRoomMersenneTwister = MersenneTwister()

function updateTileStates!(room::String, package::String, states::TileStates, width::Int, height::Int, fg::Bool=false)
    seed = getTileStateSeed(room, package, fg)
    Random.seed!(drawableRoomMersenneTwister, seed)

    states.quads = fill(Coord(-1, -1), (height, width))
    states.chars = fill('0', (height, width))
    states.rands = rand(drawableRoomMersenneTwister, 0:255, (height, width))
end

updateTileStates!(room::Room, package::String, states::TileStates, width::Int, height::Int, fg::Bool=false) = updateTileStates!(room.name, package, states, width, height, fg)

function getDrawingLayers()
    return Layer[
        Layer("bgParallax"),
        Layer("bgTiles", clearOnReset=false),
        Layer("bgDecals"),
        Layer("entities"),
        Layer("fgTiles", clearOnReset=false),
        Layer("fgParallax"),
        Layer("fgDecals"),
        Layer("triggers"),

        Layer("tools", selectable=false),
        Layer("all", dummy=true),
        Layer("objtiles", dummy=true)
    ]
end

function roomVisible(camera::Camera, width::Int, height::Int, room::Room)
    actuallX = camera.x / camera.scale
    actuallY = camera.y / camera.scale

    actuallWidth = width / camera.scale
    actuallHeight = height / camera.scale

    cameraRect = Rectangle(actuallX, actuallY, actuallWidth, actuallHeight)
    roomRect = Rectangle(room.position[1], room.position[2], room.size[1], room.size[2])

    return checkCollision(cameraRect, roomRect)
end

redrawRenderingLayer(renderingLayer::Layer, layers::Array{Layer, 1}) = renderingLayer.redraw || any(map(layer -> layer.redraw && !layer.dummy, layers))

function getRoomBackgroundColor(room::Room)::colorTupleType
    if 0 <= room.color < length(colors.background_room_color_coded_fill)
        return colors.background_room_color_coded_fill[room.color + 1]

    else
        return colors.background_room_fill
    end
end

function DrawableRoom(map::Map, room::Room)
    return DrawableRoom(
        map,
        room,

        TileStates(),
        TileStates(),

        Layer("rendering"),
        getDrawingLayers(),

        nothing
    )
end

function destroy(dr::DrawableRoom)
    deleteSurface(dr.rendering.surface)

    for layer in dr.layers
        deleteSurface(layer.surface)
    end
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
module HistoryLogger

using Printf
using RecipesBase
using Statistics
using POLO
using POLO: AbstractLogger
using POLO.Utility: LogLevel, Value, Decision, Gradient, Full, SimpleLogger
import POLO.log

struct History{L <: POLO.AbstractLogger} <: AbstractLogger
    fhist::Vector{Float64}
    inner::L

    function (::Type{History})()
        return new{SimpleLogger.Logger{None}()}(Vector{Float64}(), SimpleLogger.None())
    end

    function (::Type{History})(::Type{L}) where L <: LogLevel
        return new{SimpleLogger.Logger{L}}(Vector{Float64}(), SimpleLogger.Logger(L))
    end

    function (::Type{History})(logger::AbstractLogger)
        return new{typeof(logger)}(Vector{Float64}(), logger)
    end
end

function log(logger::History,k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector)
    push!(logger.fhist,fval)
    log(logger.inner,k,fval,x,g)
end

@recipe function f(logger::History)
    length(logger.fhist) > 0 || error("No solution data. Has solver been run?")
    fmin = minimum(logger.fhist)
    fmax = maximum(logger.fhist)
    increment = std(logger.fhist)

    linewidth --> 4
    linecolor --> :black
    tickfontsize := 14
    tickfontfamily := "sans-serif"
    guidefontsize := 16
    guidefontfamily := "sans-serif"
    titlefontsize := 22
    titlefontfamily := "sans-serif"
    xlabel := "k"
    ylabel := "f"
    ylims --> (fmin-increment,fmax+increment)
    xlims --> (1,length(logger.fhist)+1)
    xticks --> 1:5:length(logger.fhist)
    yticks --> fmin:increment:fmax
    xformatter := (d) -> @sprintf("%.1f",d)
    yformatter := (d) -> begin
        if abs(d) <= sqrt(eps())
            "0.0"
        elseif (log10(abs(d)) < -2.0 || log10(abs(d)) > 3.0)
            @sprintf("%.4e",d)
        elseif log10(abs(d)) > 2.0
            @sprintf("%.1f",d)
        else
            @sprintf("%.2f",d)
        end
    end

    @series begin
        seriescolor --> :black
        1:1:length(logger.fhist),logger.fhist
    end
end

end


text:
"""
    AbstractSystem

Abstract supertype for all system types.
"""
abstract type AbstractSystem end

"""
    statedim(s::AbstractSystem)

Returns the dimension of the state space of system `s`.
"""
function statedim end

"""
    stateset(s::AbstractSystem)

Returns the set of allowed states of system `s`.
"""
function stateset end

"""
    inputdim(s::AbstractSystem)

Returns the dimension of the input space of system `s`.
"""
function inputdim end

"""
    inputset(s::AbstractSystem)

Returns the set of allowed inputs of system `s`.
"""
function inputset end

"""
    AbstractDiscreteSystem

Abstract supertype for all discrete system types.
"""
abstract type AbstractDiscreteSystem <: AbstractSystem end

"""
    AbstractContinuousSystem

Abstract supertype for all continuous system types.
"""
abstract type AbstractContinuousSystem <: AbstractSystem end

"""
    islinear(s::AbstractSystem)

Specifies if the dynamics of system `s` is specified by linear equations.

### Notes

We adopt the notion from [Section 2.7, 1]. For example, the system with inputs
``x' = f(t, x, u) = A x + B u`` is linear, since the function ``f(t, ⋅, ⋅)`` is
linear in ``(x, u)`` for each ``t ∈ \\mathbb{R}``. On the other hand,
``x' = f(t, x, u) = A x + B u + c`` is affine but not linear, since it is not
linear in ``(x, u)``.

This function uses the information of the type, not the value. So, if a system
type allows an instance that is not linear, it returns `false` by default.
For example, polynomial systems can be nonlinear; hence `islinear`
returns `false`.

[1] Sontag, Eduardo D. *Mathematical control theory: deterministic finite dimensional
systems.* Vol. 6. Springer Science & Business Media, 2013.
"""
function islinear(::AbstractSystem) end

"""
    isaffine(s::AbstractSystem)

Specifies if the dynamics of system `s` is specified by affine equations.

### Notes

An affine system is the composition of a linear system and a translation.
See [`islinear(::AbstractSystem)`](@ref) for the notion of linear system adopted
in this library.

This function uses the information of the type, not the value. So, if a system
type allows an instance that is not affine, it returns `false` by default.
For example, polynomial systems can be nonlinear; hence `isaffine` is `false`.
"""
function isaffine(::AbstractSystem) end

"""
    AbstractMap

Abstract supertype for all map types.
"""
abstract type AbstractMap end

"""
    outputdim(m::AbstractMap)

Returns the dimension of the output space of the map `m`.
"""
function outputdim end

"""
    outputmap(s::SystemWithOutput)

Returns the output map of a system with output.
"""
function outputmap end

"""
    islinear(m::AbstractMap)

Specifies if the map `m` is linear or not.

### Notes

A map is *linear* if it preserves the operations of scalar multiplication and
vector addition.
"""
function islinear(::AbstractMap) end

"""
    isaffine(m::AbstractMap)

Specifies if the map `m` is affine or not.

### Notes

An affine map is the composition of a linear map and a translation.
See also [`islinear(::AbstractMap)`](@ref).
"""
function isaffine(::AbstractMap) end

"""
    apply(m::AbstractMap, args...)

Apply the rule specified by the map to the given arguments.
"""
function apply end


text:
module TestBasics

import Schedulers
using ArgCheck: @check
using Schedulers.Internal: @record
using Test

using ..Utils: fetch_finished

function check_simple_spawn(scheduler)
    @record(:check_simple_spawn_begin, scheduler)
    local t1, t2
    Schedulers.open(scheduler) do _
        t1 = Schedulers.spawn(() -> 1)
        t2 = Schedulers.spawn(() -> fetch(t1) + 2)
        wait(t2)
    end
    @record(:check_simple_spawn_end, scheduler)
    @test fetch_finished(t1) == 1
    @test fetch_finished(t2) == 3
end

function test_simple_spawn()
    @testset "$scheduler" for scheduler in [Schedulers.workstealing, Schedulers.prioritized]
        check_simple_spawn(scheduler())
    end
end

end  # module


text:
sorttable(T; ordering=<, column=1, reverse=false) =
  sort(T, by = t -> t[column], lt = reverse ? (a,b) -> ordering(b,a) : ordering)
sorttable(T, ordering=<, column=1, reverse=false) =
  sorttable(T, ordering=ordering, column=column, reverse=reverse)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Unzip

import Base:
    axes,
    _collect,
    empty,
    getindex,
    push!,
    push_widen,
    setindex!,
    setindex_widen_up_to,
    similar,
    size
using Base:
    EltypeUnknown,
    HasEltype,
    HasLength,
    isabstracttype,
    isvatuple,
    IteratorEltype,
    IteratorSize,
    @propagate_inbounds,
    @pure,
    SizeUnknown,
    tail

@inline function all_unrolled(call)
    true
end
@inline function all_unrolled(call, item, rest...)
    if call(item)
        all_unrolled(call, rest...)
    else
        false
    end
end

@inline same_axes() = true
function same_axes(call, first_column, rest...)
    all_unrolled(let first_axes = axes(first_column)
        column -> axes(column) == first_axes
    end, rest...)
end

struct Rows{Row, Dimensions, Columns} <: AbstractArray{Row, Dimensions}
    columns::Columns
end

@propagate_inbounds function Rows{Row, Dimension}(
    columns::Columns,
) where {Row, Dimension, Columns}
    @boundscheck if !same_axes(columns...)
        throw(DimensionMismatch("All arguments to `Rows` must have the same axes"))
    end
    Rows{Row, Dimension, Columns}(columns)
end

model_column(columns) = first(columns)
model_column(::Tuple{}) = 1:0
model_column(rows::Rows) = model_column(rows.columns)

@propagate_inbounds function Rows(columns)
    Rows{Tuple{map(eltype, columns)...}, ndims(model_column(columns))}(columns)
end

axes(rows::Rows, dimensions...) = axes(model_column(rows), dimensions...)
size(rows::Rows, dimensions...) = size(model_column(rows), dimensions...)

@propagate_inbounds function getindex(rows::Rows, an_index...)
    map(let an_index = an_index
        @propagate_inbounds function (column)
            column[an_index...]
        end
    end, rows.columns)
end

@propagate_inbounds function setindex!(rows::Rows, row, an_index...)
    map(let an_index = an_index
        @propagate_inbounds function (column, value)
            column[an_index...] = value
        end
    end, rows.columns, row)
    nothing
end

function push!(rows::Rows, row)
    map(push!, rows.columns, row)
    nothing
end

# can we do any better?
@pure function can_guess_column_types(_)
    false
end
@pure function can_guess_column_types(row_type::DataType)
    !(isabstracttype(row_type) || (row_type.name == Tuple.name && isvatuple(row_type)))
end
@pure function val_fieldtypes(row_type)
    if can_guess_column_types(row_type)
        map(Val, (row_type.types...,))
    else
        ()
    end
end

function similar(rows::Rows, ::Type{ARow}, dimensions::Dims) where {ARow}
    @inbounds Rows(map(
        let model = model_column(rows), dimensions = dimensions
            function (::Val{Value},) where {Value}
                similar(model, Value, dimensions)
            end
        end,
        val_fieldtypes(ARow),
    ))
end

function empty(column::Rows{OldRow}, ::Type{NewRow} = OldRow) where {OldRow, NewRow}
    similar(column, NewRow)
end

function widen_column(
    ::HasLength,
    new_length,
    an_index,
    column::Array{Element},
    item::Item,
) where {Element, Item <: Element}
    @inbounds column[an_index] = item
    column
end
function widen_column(::HasLength, new_length, an_index, column::Array, item)
    setindex_widen_up_to(column, item, an_index)
end

function widen_column(
    ::SizeUnknown,
    new_length,
    an_index,
    column::Array{Element},
    item::Item,
) where {Element, Item <: Element}
    push!(column, item)
    column
end
function widen_column(::SizeUnknown, new_length, an_index, column::Array, item)
    push_widen(column, item)
end

function widen_column(
    iterator_size,
    new_length,
    an_index,
    ::Missing,
    item::Item,
) where {Item}
    new_column = Array{Union{Missing, Item}}(missing, new_length)
    @inbounds new_column[an_index] = item
    new_column
end
function widen_column(iterator_size, new_length, an_index, ::Missing, ::Missing)
    Array{Missing}(missing, new_length)
end

get_new_length(::SizeUnknown, rows, an_index) = an_index
get_new_length(::HasLength, rows, an_index) = length(rows)

zip_missing(::Tuple{}, ::Tuple{}) = ()
zip_missing(::Tuple{}, longer) = map(second_one -> (missing, second_one), longer)
zip_missing(longer, ::Tuple{}) = map(first_one -> (first_one, missing), longer)
function zip_missing(tuple1, tuple2)
    (first(tuple1), first(tuple2)), zip_missing(tail(tuple1), tail(tuple2))...
end

function widen_columns(iterator_size, rows, row, an_index = length(rows) + 1)
    columns = rows.columns
    @inbounds Rows(map(
        let iterator_size = iterator_size,
            new_length = get_new_length(iterator_size, rows, an_index),
            an_index = an_index

            function (column_item,)
                widen_column(iterator_size, new_length, an_index, column_item...)
            end
        end,
        zip_missing(rows.columns, row),
    ))
end

push_widen(rows::Rows, row) = widen_columns(SizeUnknown(), rows, row)

function setindex_widen_up_to(rows::Rows, row, an_index)
    widen_columns(HasLength(), rows, row, an_index)
end

"""
    unzip(rows)

Collect into columns.

```jldoctest
julia> using Unzip

julia> using Base: Generator

julia> using Test: @inferred

julia> stable(x) = (x, x + 0.0, x, x + 0.0, x, x + 0.0);

julia> @inferred unzip(Generator(stable, 1:4))
([1, 2, 3, 4], [1.0, 2.0, 3.0, 4.0], [1, 2, 3, 4], [1.0, 2.0, 3.0, 4.0], [1, 2, 3, 4], [1.0, 2.0, 3.0, 4.0])

julia> unstable(x) =
           if x == 2
               (x, x + 0.0, x, x + 0.0)
           else
               (x, x + 0.0)
           end;

julia> unzip(Generator(unstable, 1:3))
([1, 2, 3], [1.0, 2.0, 3.0], Union{Missing, Int64}[missing, 2, missing], Union{Missing, Float64}[missing, 2.0, missing])

julia> unzip(Iterators.filter(row -> true, Generator(unstable, 1:3)))
([1, 2, 3], [1.0, 2.0, 3.0], Union{Missing, Int64}[missing, 2, missing], Union{Missing, Float64}[missing, 2.0, missing])

julia> unzip(Iterators.filter(row -> false, Generator(unstable, 1:4)))
()

julia> unzip(x for x in [(1, 2), ("a", "b")])
(Any[1, "a"], Any[2, "b"])

julia> unzip([(a=1, b=2), (a="a", b="b")])
(Any[1, "a"], Any[2, "b"])
```
"""
function unzip(rows)
    iterator_eltype = IteratorEltype(rows)
    _collect(
        (@inbounds Rows(())),
        rows,
        # if we can't guess the column types, ignore the eltype
        if iterator_eltype isa HasEltype && !can_guess_column_types(eltype(rows))
            EltypeUnknown()
        else
            iterator_eltype
        end,
        IteratorSize(rows)
    ).columns
end

export unzip

end


text:
using Clustering, DataFrames

n = 1_000

srand(1)

data = randn(n, 2)

centers = [0.0 0.0;
           5.0 5.0;
           10.0 0.0;
           15.0 -5.0;]

assignments = zeros(Int64, n)

for i = 1:n
    assignments[i] = rand(1:4)
    data[i, :] += centers[assignments[i], :]
end

results = dp_means(data, 50.0)

df = DataFrame()
df["Label"] = assignments
df["Cluster"] = results.assignments
by(df, ["Label", "Cluster"], nrow)

# Should really look at active clusters, not created clusters


text:
@testset "IfThenElse2Reification: $(fct_type), $(T)" for fct_type in ["vector of variables", "vector affine function"], T in [Int, Float64]
    base_model = if T == Int
        IntReificationEqualToModel{T}()
    elseif T == Float64
        FloatReificationEqualToModel{T}()
    else
        @assert false
    end
    mock = MOIU.MockOptimizer(base_model)
    model = COIB.IfThenElse2Reification{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Reification{MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.IfThenElse{MOI.LessThan{T}, MOI.LessThan{T}, MOI.LessThan{T}},
    )

    if T == Int
        x_1, _ = MOI.add_constrained_variable(model, MOI.Integer())
        x_2, _ = MOI.add_constrained_variable(model, MOI.Integer())
        x_3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    elseif T == Float64
        x_1 = MOI.add_variable(model)
        x_2 = MOI.add_variable(model)
        x_3 = MOI.add_variable(model)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x_1, x_2, x_3])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x_1, x_2, x_3])
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.IfThenElse(MOI.LessThan(one(T)), MOI.LessThan(zero(T)), MOI.LessThan(zero(T))))

    @test MOI.is_valid(model, x_1)
    @test MOI.is_valid(model, x_2)
    @test MOI.is_valid(model, x_3)
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.IfThenElse}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.IfThenElse{MOI.LessThan{T}, MOI.LessThan{T}, MOI.LessThan{T}}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Reification),
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 3
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == 3
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Reification}()) == 3
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == 2

        @test MOI.get(bridge, MOI.ListOfVariableIndices()) == [bridge.var_condition, bridge.var_true, bridge.var_false]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}()) == [bridge.var_condition_bin, bridge.var_true_bin, bridge.var_false_bin]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Reification}()) == [bridge.con_reif_condition, bridge.con_reif_true, bridge.con_reif_false]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == [bridge.con_if, bridge.con_else]
    end

    @testset "Set of variables" begin
        @test MOI.is_valid(model, bridge.var_condition)
        @test MOI.is_valid(model, bridge.var_condition_bin)
        @test MOI.get(model, MOI.ConstraintFunction(), bridge.var_condition_bin) == bridge.var_condition
        @test MOI.get(model, MOI.ConstraintSet(), bridge.var_condition_bin) == MOI.ZeroOne()

        @test MOI.is_valid(model, bridge.var_true)
        @test MOI.is_valid(model, bridge.var_true_bin)
        @test MOI.get(model, MOI.ConstraintFunction(), bridge.var_true_bin) == bridge.var_true
        @test MOI.get(model, MOI.ConstraintSet(), bridge.var_true_bin) == MOI.ZeroOne()

        @test MOI.is_valid(model, bridge.var_false)
        @test MOI.is_valid(model, bridge.var_false_bin)
        @test MOI.get(model, MOI.ConstraintFunction(), bridge.var_false_bin) == bridge.var_false
        @test MOI.get(model, MOI.ConstraintSet(), bridge.var_false_bin) == MOI.ZeroOne()
    end

    @testset "Reification of the condition" begin
        @test MOI.is_valid(model, bridge.con_reif_condition)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_reif_condition)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_reif_condition) == CP.Reification(MOI.LessThan(one(T)))

        t1 = f.terms[1]
        @test t1.output_index == 1
        @test t1.scalar_term.coefficient === one(T)
        @test t1.scalar_term.variable == bridge.var_condition

        t2 = f.terms[2]
        @test t2.output_index == 2
        @test t2.scalar_term.coefficient === one(T)
        @test t2.scalar_term.variable == x_1
    end

    @testset "Reification of the true constraint" begin
        @test MOI.is_valid(model, bridge.con_reif_true)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_reif_true)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_reif_true) == CP.Reification(MOI.LessThan(zero(T)))

        t1 = f.terms[1]
        @test t1.output_index == 1
        @test t1.scalar_term.coefficient === one(T)
        @test t1.scalar_term.variable == bridge.var_true

        t2 = f.terms[2]
        @test t2.output_index == 2
        @test t2.scalar_term.coefficient === one(T)
        @test t2.scalar_term.variable == x_2
    end

    @testset "Reification of the false constraint" begin
        @test MOI.is_valid(model, bridge.con_reif_false)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_reif_false)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_reif_false) == CP.Reification(MOI.LessThan(zero(T)))

        t1 = f.terms[1]
        @test t1.output_index == 1
        @test t1.scalar_term.coefficient === one(T)
        @test t1.scalar_term.variable == bridge.var_false

        t2 = f.terms[2]
        @test t2.output_index == 2
        @test t2.scalar_term.coefficient === one(T)
        @test t2.scalar_term.variable == x_3
    end

    @testset "If" begin
        @test MOI.is_valid(model, bridge.con_if)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_if)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_if) == MOI.GreaterThan(zero(T))

        t1 = f.terms[1]
        @test t1.coefficient === -one(T)
        @test t1.variable == bridge.var_condition

        t2 = f.terms[2]
        @test t2.coefficient === one(T)
        @test t2.variable == bridge.var_true
    end

    @testset "Else" begin
        @test MOI.is_valid(model, bridge.con_else)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_else)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_else) == MOI.GreaterThan(one(T))

        t1 = f.terms[1]
        @test t1.coefficient === one(T)
        @test t1.variable == bridge.var_condition

        t2 = f.terms[2]
        @test t2.coefficient === one(T)
        @test t2.variable == bridge.var_false
    end
end


text:
using JuMP, Gurobi
include("../basics/main.jl")
include("../basics/validations.jl")
include("../JuMP/JuMP_models.jl")
include("../heuristics/construction_algorithms.jl")
include("../heuristics/local_search.jl")

function GRASP_final_local_search(CTS::Constants)
    tasks_by_w = OperationalShippingProblem(beta, CTS)
    final_update = "false"
    best_makespan = CTS.H
    best_LS = init_ls(CTS)
    for iter = 1:1000
        makespan = CTS.H
        LS = init_ls(CTS)
        QC = init_qc(CTS)
        TIME = init_timer(CTS)

        while LS.tasks_left > 0 && TIME.period < TIME.horizon_plan
            cranes_status = "idle"
            while cranes_status != "All cranes are busy" && cranes_status != "Next time period" && cranes_status != "LS is completed"
                cranes_status = GRASP("dist", tasks_by_w, bj, LS, TIME, QC, CTS)
            end
            next_time_period(TIME, QC);
        end

        if check_solution(LS, QC, CTS) == true
            makespan = total_makespan(LS, CTS)
            if makespan < best_makespan
                best_makespan = makespan
                best_LS = LS
            end

            #perform local_search
            final_LS = init_ls(CTS)
            final_makespan = CTS.H
            for iter_ls = 1:5
                min_q, moves, TIME, QC = get_final_state(best_LS, CTS)
                if moves == 2
                    break
                end
                final_LS = remove_final_tasks(moves, best_LS, CTS)
                while final_LS.tasks_left > 0 && TIME.period < TIME.horizon_plan
                    cranes_status = "idle"
                    while cranes_status != "All cranes are busy" && cranes_status != "Next time period" && cranes_status != "LS is completed"
                        cranes_status = ls_final(min_q, tasks_by_w, bj, final_LS, TIME, QC, CTS)
                    end
                    next_time_period(TIME, QC);
                end

                if check_solution(final_LS, QC, CTS) == true
                    final_makespan = total_makespan(final_LS, CTS)
                    if final_makespan <= best_makespan
                        best_makespan = final_makespan
                        best_LS = deepcopy(final_LS)
                        final_update = "true"
                    end
                end
            end
        end
    end

    return(final_update, best_makespan, best_LS)
end

final_update, best_makespan, best_LS = GRASP_final_local_search(CTS)
plot_solution(best_LS, best_makespan, CTS)

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
function powerstoragevariables(m::JuMP.AbstractModel, devices_netinjection:: A, devices::Array{T,1}, time_periods::Int64) where {A <: JumpExpressionMatrix, T <: PowerSystems.Storage}

    on_set = [d.name for d in devices if d.available]
    t = 1:time_periods

    pstin = @variable(m, pstin[on_set,t])
    pstout = @variable(m, pstout[on_set,t])

    devices_netinjection = varnetinjectiterate!(devices_netinjection,  pstin, pstout, t, devices)

    return pstin, pstout, devices_netinjection
end

function energystoragevariables(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64) where T <: PowerSystems.Storage

    on_set = [d.name for d in devices if d.available]
    t = 1:time_periods

    ebt = @variable(m, ebt[on_set,t] >= 0.0)

    return ebt
end