text:
using JuLIP, NBodyIPs, IPFitting, StaticArrays, Plots

include(homedir() * "/Dropbox/PIBmat/W_Data/W.jl")
E0 = W.get_E0()

db = LsqDB(homedir() * "/scratch/nbodyips/W_4BBL_2s_med")

db

dataweights = Dict("E" => 50.0, "F" => 1.0, "V" => 1.0)
configweights = Dict("dia" => (1.0),
                     "surface_110" => 1.0,
                     "surface_001" => 1.0,
                     "surface_111" => 1.0,
                     )

r0 = rnn(:Si)
rcuts = unique(cutoff.(db.basis))
rcut4, rcut3, rcut2 = (sort(rcuts)...,)

reg = [ BAReg(2, 0.5*r0, 0.85*r0, creg = 1.0),
        BAReg(3, 0.7*r0, 0.85*r0, creg = 0.1),
        BAReg(4, 0.7*r0, 0.85*r0, creg = 0.1),
        ]

IPreg, inforeg = lsqfit( db; E0 = E0, dataweights = dataweights,
                           configweights = configweights,
                           regularisers = reg )

V20 = IPreg.components[2]
V21 = IPreg.components[3]
r = range(0.8*r0, stop=rcut2, length=200)
P = plot(r, V20.Vr.(r) );
plot!(P, r, V21.Vr.(r)  );
display(P)


text:
using jInv.Mesh
using EikonalInv
using jInv.ForwardShare
using jInv.Utils

using Base.Test

function testSensMatVec(pFor,m0)

	println("\t\tchecking derivatives")

	dobs, = getData(vec(m0),pFor)

	function testFun(m,v=[])
		dobs, = getData(vec(m),pFor)
		if isempty(v)
			return vec(dobs)
		else
			JTv = getSensMatVec(vec(v),vec(m),pFor)
			return vec(dobs), JTv
		end
	end

	chkDer, = jInv.Utils.checkDerivative(testFun,vec(m0))
	@test chkDer

	(nd,nm) = getSensMatSize(pFor)
	v0 = randn(nm)
	w0 = randn(nd)

	J  = getSensMat(vec(m0),pFor)
	J1 = J*v0;
	J2 = J'*w0;

	J1mf = getSensMatVec(v0,vec(m0),pFor)
	J2mf = getSensTMatVec(w0,vec(m0),pFor)

	println("\t\ttest matrix free sensitivities")
	@test norm(J1mf-J1)/norm(J1) < 1e-10
	@test norm(J2mf-J2)/norm(J2) < 1e-10

	println("\t\tadjoint test")
	t1 = dot(w0,J1mf)
	t2 = dot(v0,J2mf)

	p = abs(t1-t2)/abs(t1);
	@test p < 1e-10
	if p > 1e-10
		println("Inconsistency in sensitivity: ",p);
	end

end
############################################################################################
# get 2D mesh with uneven number of cells
domain = [1. 4.0 .4 4.2]
n      = [17;16;]
M      = getRegularMesh(domain,n)
# sources / receivers on top edge
idx    = reshape(collect(1:prod(n+1)),tuple(n+1...))
ib     = idx[1:2:end,1];
I      = speye(prod(n+1))
Q      = I[:,vec(ib)]
R      = copy(Q)

println("2D - test low order scheme")
pFor = getEikonalInvParam(M,Q,R,false)
m0   = rand(prod(n+1))+10.0

try
	testSensMatVec(pFor,vec(m0))
catch
	warn("Sensitivity for low-order scheme and mesh size n=$n FAILED!")
end

############################################################################################
println("2D - test high order scheme")
pFor = getEikonalInvParam(M,Q,R,true)
m0   = rand(prod(n+1))+10.0
try
	testSensMatVec(pFor,vec(m0))
catch
	warn("Sensitivity for higher-order scheme and mesh size n=$n FAILED!")
end

n      = [16;32;]
M      = getRegularMesh(domain,n)
# sources / receivers on top edge
idx    = reshape(collect(1:prod(n+1)),tuple(n+1...))
ib     = idx[1:2:end,1];
I      = speye(prod(n+1))
Q      = I[:,vec(ib)]
R      = copy(Q)
println("2D - test low order scheme")
pFor = getEikonalInvParam(M,Q,R,false)
m0   = rand(prod(n+1))+10.0
try
	testSensMatVec(pFor,vec(m0))
catch
	warn("Sensitivity for lower-order scheme and mesh size n=$n FAILED!")
end

println("2D - test high order scheme")
pFor = getEikonalInvParam(M,Q,R,true)
m0   = rand(prod(n+1))+10.0
try
	testSensMatVec(pFor,vec(m0))
catch
	warn("Sensitivity for higher-order scheme and mesh size n=$n FAILED!")
end
############################################################################################

# get 3D mesh
domain = [1. 4.0 .4 4.2 0 1.2]
n      = [16;16;8]
M      = getRegularMesh(domain,n)
# sources / receivers on top edge
idx    = reshape(collect(1:prod(n+1)),tuple(n+1...))
ib     = idx[1:6:end,1,1:4:end];
I      = speye(prod(n+1))
Q      = I[:,vec(ib)]
R      = copy(Q)
println("3D - test low order scheme")
pFor = getEikonalInvParam(M,Q,R,false)
m0   = rand(prod(n+1))+10.0
try
	testSensMatVec(pFor,vec(m0))
catch
	warn("Sensitivity for low-order scheme and mesh size n=$n FAILED!")
end

println("3D - test high order scheme")
pFor = getEikonalInvParam(M,Q,R,true)
m0   = rand(prod(n+1))+10.0
try
	testSensMatVec(pFor,vec(m0))
catch
	warn("Sensitivity for higher-order scheme and mesh size n=$n FAILED!")
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_ApplicationGatewayAuthenticationCertificatePropertiesFormat.jl")
include("model_ApplicationGatewayAvailableWafRuleSetsResult.jl")
include("model_ApplicationGatewayBackendAddress.jl")
include("model_ApplicationGatewayBackendAddressPoolPropertiesFormat.jl")
include("model_ApplicationGatewayBackendHealth.jl")
include("model_ApplicationGatewayBackendHealthHttpSettings.jl")
include("model_ApplicationGatewayBackendHealthPool.jl")
include("model_ApplicationGatewayBackendHealthServer.jl")
include("model_ApplicationGatewayBackendHttpSettingsPropertiesFormat.jl")
include("model_ApplicationGatewayConnectionDraining.jl")
include("model_ApplicationGatewayFirewallDisabledRuleGroup.jl")
include("model_ApplicationGatewayFirewallRule.jl")
include("model_ApplicationGatewayFirewallRuleGroup.jl")
include("model_ApplicationGatewayFirewallRuleSetPropertiesFormat.jl")
include("model_ApplicationGatewayFrontendIPConfigurationPropertiesFormat.jl")
include("model_ApplicationGatewayFrontendPortPropertiesFormat.jl")
include("model_ApplicationGatewayHttpListenerPropertiesFormat.jl")
include("model_ApplicationGatewayIPConfigurationPropertiesFormat.jl")
include("model_ApplicationGatewayListResult.jl")
include("model_ApplicationGatewayPathRulePropertiesFormat.jl")
include("model_ApplicationGatewayProbePropertiesFormat.jl")
include("model_ApplicationGatewayPropertiesFormat.jl")
include("model_ApplicationGatewayRequestRoutingRulePropertiesFormat.jl")
include("model_ApplicationGatewaySku.jl")
include("model_ApplicationGatewaySslCertificatePropertiesFormat.jl")
include("model_ApplicationGatewaySslPolicy.jl")
include("model_ApplicationGatewayUrlPathMapPropertiesFormat.jl")
include("model_ApplicationGatewayWebApplicationFirewallConfiguration.jl")
include("model_Resource.jl")
include("model_SubResource.jl")
include("model_ApplicationGateway.jl")
include("model_ApplicationGatewayAuthenticationCertificate.jl")
include("model_ApplicationGatewayBackendAddressPool.jl")
include("model_ApplicationGatewayBackendHttpSettings.jl")
include("model_ApplicationGatewayFirewallRuleSet.jl")
include("model_ApplicationGatewayFrontendIPConfiguration.jl")
include("model_ApplicationGatewayFrontendPort.jl")
include("model_ApplicationGatewayHttpListener.jl")
include("model_ApplicationGatewayIPConfiguration.jl")
include("model_ApplicationGatewayPathRule.jl")
include("model_ApplicationGatewayProbe.jl")
include("model_ApplicationGatewayRequestRoutingRule.jl")
include("model_ApplicationGatewaySslCertificate.jl")
include("model_ApplicationGatewayUrlPathMap.jl")
include("model_NetworkInterfaceIPConfiguration.jl")


text:
module Kalman

export State,
AbstractState,
UnscentedState,
Observation,
KalmanFilter,
LinearKalmanFilter,
BasicKalmanFilter,
Model,
ObservationModel,
LinearModel,
LinearObservationModel,
predict,
update,
predictupdate,
predict!,
update!,
predictupdate!,
BasicExtendedKalmanFilter,
NonlinearModel,
NonlinearObservationModel,
ExtendedKalmanFilter,
AdditiveUnscentedKalmanFilter,
UnscentedKalmanFilter,
AdditiveUnscentedModel,
AdditiveUnscentedObservationModel,
sigma,
AugmentedUnscentedModel,
AugmentedUnscentedObservationModel,
AugmentedUnscentedKalmanFilter,
augment

include("bkf.jl")
include("filter.jl")
include("ekf.jl")
include("unscented.jl")
include("augmented.jl")

# Nullable Observations work but only for Julia v0.4
if VERSION >= v"0.4.0-"
    include("nullable.jl")
end

end # module


text:
module ego_model
export EgoModel

using Knet

push!(LOAD_PATH, pwd());
using models: Dense

struct EgoModel
    dt::Float32
    front_wb::Param
    rear_wb::Param
    steer_gain::Param
    brake_accel::Param
    throt_accel::Dense

    function EgoModel(dt=1/4)
        _atype = Knet.CUDA.functional() ? KnetArray{Float32} : Array{Float32}
        front_wb = Param(_atype(ones(1)))
        rear_wb = Param(_atype(ones(1)))

        steer_gain = Param(_atype(ones(1)))
        brake_accel = Param(_atype(zeros(1)))

        throt_accel = Dense(1,1)
        new(dt,front_wb,rear_wb,steer_gain,brake_accel, throt_accel)
    end
end

(m::EgoModel)(locs, yaws, spds, acts) =
begin
    steer = selectdim(acts, length(size(acts)), 1:1)
    throt = selectdim(acts, length(size(acts)), 2:2)
    brake = convert(KnetArray{UInt8}, selectdim(acts, length(size(acts)), 3:3))

    # accel = #TODO: do acceleration calculation
    wheel = m.steer_gain .* steer

    beta = atan.(m.rear_wb/(m.front_wb + m.rear_wb) .* wheel)
    next_locs = locs + spds .* cat(cos.(yaws+beta), sin.(yaws+beta);dims=length(size(beta)))
    next_yaws = yaws + spds ./ m.rear_wb .* sin.(beta) * m.dt
    next_spds = spds + dt .* accel
    return next_locs, next_yaws, relu(next_spds)
end
end


text:
using CorrelationView
using Base.Test

@test true


text:
mutable struct Shader <: Resource
	program::GLuint
	uniforms::Dict{Symbol, UniformVar}
	samplers::Vector{Symbol}
	blocks::Vector{UniformBlock}
	setupMaterial::Function
	id::Symbol
	attribType::DataType
	renderer::Renderer

	Shader() = new(0, Dict{Symbol, UniformVar}(), Symbol[], UniformBlock[])
end

isvalid(shader::Shader) = shader.program != 0

function init(shader::Shader, renderer::Renderer, path::AbstractString, setupMaterial::Function = identity; id::Symbol = Symbol(path))
	local vsSource, psSource
	open(path * ".vert") do f
		vsSource = vcat(renderer.shaderPreamble, read(f))
	end

	open(path * ".frag") do f
		psSource = vcat(renderer.shaderPreamble, read(f))
	end
	init(shader, renderer, vsSource, psSource, setupMaterial, id = id)
end

function init(shader::Shader, renderer::Renderer, vs::AbstractArray{UInt8}, ps::AbstractArray{UInt8}, setupMaterial::Function = identity; id::Symbol = :shader)
	@assert !isvalid(shader)

	vertexShader = compileshader(GL_VERTEX_SHADER, vs, id)
	if vertexShader != 0
		fragmentShader = compileshader(GL_FRAGMENT_SHADER, ps, id)
		if fragmentShader != 0
			shader.program = glCreateProgram()

			glAttachShader(shader.program, vertexShader)
			glAttachShader(shader.program, fragmentShader)

			glLinkProgram(shader.program)

			linkSuccess = GLint[0]
			glGetProgramiv(shader.program, GL_LINK_STATUS, linkSuccess)
			if linkSuccess[1] != GL_TRUE
				msg = get_program_info_log(shader.program)
				@info("Error linking shader program id $id\n$msg")
				glDeleteProgram(shader.program)
				shader.program = 0
			end

			glDeleteShader(fragmentShader)
		end
		glDeleteShader(vertexShader)
	end

	if isvalid(shader)
		shader.setupMaterial = setupMaterial
		init_resource(shader, renderer, id)
		initblocks(shader)
		inituniforms(shader)

		initattributes(shader)
	end
	shader
end

function done(shader::Shader)
	if isvalid(shader)
		doneblocks(shader)
		glDeleteProgram(shader.program)
		shader.program = 0
		empty!(shader.uniforms)
		empty!(shader.samplers)
		remove_renderer_resource(shader)
	end
end

function apply(shader::Shader)
	@assert isvalid(shader)

	# possibly use glValidateProgram first?

	glUseProgram(shader.program)

	for block in shader.blocks
		apply(block)
	end
end

function get_info_log(glObj::GLuint, get::Function, getInfoLog::Function)
	logLength = GLint[0]
	get(glObj, GL_INFO_LOG_LENGTH, logLength)

	message = Array{UInt8}(undef, logLength[1])
	getInfoLog(glObj, logLength[1], C_NULL, message)
	String(message)
end

get_shader_info_log(shaderObj::GLuint) = get_info_log(shaderObj, glGetShaderiv, glGetShaderInfoLog)
get_program_info_log(programObj::GLuint) = get_info_log(programObj, glGetProgramiv, glGetProgramInfoLog)

function compileshader(shaderType::UInt32, source::AbstractArray{UInt8}, id::Symbol)
	shaderObj = glCreateShader(shaderType)

	srcArray = Ptr{UInt8}[pointer(source)]
	lenArray = GLint[length(source)]
	glShaderSource(shaderObj, 1, srcArray, lenArray)

	glCompileShader(shaderObj)

	compileSuccess = GLint[0]
	glGetShaderiv(shaderObj, GL_COMPILE_STATUS, compileSuccess)

	if compileSuccess[1] != GL_TRUE
		msg = get_shader_info_log(shaderObj)
		shaderTypeName = shaderType == GL_VERTEX_SHADER ? "VERTEX" : (shaderType == GL_FRAGMENT_SHADER ? "FRAGMENT" : "$shaderType")
		@info("Error compiling shader type $shaderTypeName with id $id\n$msg")

		glDeleteShader(shaderObj)
		shaderObj = 0
	end

	return shaderObj
end

function initblocks(shader::Shader)
	@assert isvalid(shader)
	@assert isempty(shader.blocks)

	blockCount = GLint[0]
	glGetProgramiv(shader.program, GL_ACTIVE_UNIFORM_BLOCKS, blockCount)
	@assert glGetError() == GL_NO_ERROR

	for i::GLuint = 0:blockCount[1]-1
		block = UniformBlock()
		init(block, shader.program, i)
		@assert isvalid(block)
		push!(shader.blocks, block)
	end
end

function get_block_id(shader::Shader, index::GLint)
	if index >= 0
		found = findfirst(b -> b.index == index, shader.blocks)
		found == nothing && return -1
		return found
	else
		return -1
	end
end

function inituniforms(shader::Shader)
	@assert isvalid(shader)
	@assert isempty(shader.uniforms)
	@assert isempty(shader.samplers)

	val = GLint[0]
	glGetProgramiv(shader.program, GL_ACTIVE_UNIFORMS, val)
	@assert glGetError() == GL_NO_ERROR

	uniformIndices = GLuint[i for i in 0:val[1]-1]
	count = length(uniformIndices)
	@assert count == val[1]

	uniformTypes = Array{GLint}(undef, count)
	glGetActiveUniformsiv(shader.program, count, uniformIndices, GL_UNIFORM_TYPE, uniformTypes)
	@assert glGetError() == GL_NO_ERROR

	uniformOffsets = Array{GLint}(undef, count)
	glGetActiveUniformsiv(shader.program, count, uniformIndices, GL_UNIFORM_OFFSET, uniformOffsets)
	@assert glGetError() == GL_NO_ERROR

	uniformArraySizes = Array{GLint}(undef, count)
	glGetActiveUniformsiv(shader.program, count, uniformIndices, GL_UNIFORM_SIZE, uniformArraySizes)
	@assert glGetError() == GL_NO_ERROR

	uniformArrayStrides = Array{GLint}(undef, count)
	glGetActiveUniformsiv(shader.program, count, uniformIndices, GL_UNIFORM_ARRAY_STRIDE, uniformArrayStrides)
	@assert glGetError() == GL_NO_ERROR

	uniformMatrixStrides = Array{GLint}(undef, count)
	glGetActiveUniformsiv(shader.program, count, uniformIndices, GL_UNIFORM_MATRIX_STRIDE, uniformMatrixStrides)
	@assert glGetError() == GL_NO_ERROR

	uniformBlockIndices = Array{GLint}(undef, count)
	glGetActiveUniformsiv(shader.program, count, uniformIndices, GL_UNIFORM_BLOCK_INDEX, uniformBlockIndices)
	@assert glGetError() == GL_NO_ERROR

	glGetProgramiv(shader.program, GL_ACTIVE_UNIFORM_MAX_LENGTH, val)
	@assert glGetError() == GL_NO_ERROR
	nameBuf = Array{UInt8}(undef, val[1])

	# set program so we can set sampler uniforms
	glUseProgram(shader.program)
	@assert glGetError() == GL_NO_ERROR
	nameLen = GLsizei[0]
	varSize = GLsizei[0]
	varType = GLenum[0]
	for i in 1:count
		glGetActiveUniform(shader.program, uniformIndices[i], length(nameBuf), nameLen, varSize, varType, nameBuf)
		@assert glGetError() == GL_NO_ERROR

		blockId = get_block_id(shader, uniformBlockIndices[i])
		loc = -1
		if blockId <= 0
			loc = glGetUniformLocation(shader.program, nameBuf)
			@assert glGetError() == GL_NO_ERROR
			@assert loc >= 0
		end

		var = UniformVar(Symbol(unsafe_string(pointer(nameBuf))),
						 gl2jltype(uniformTypes[i]),
						 loc,
						 uniformOffsets[i],
						 uniformArraySizes[i],
						 uniformArrayStrides[i],
						 uniformMatrixStrides[i],
						 blockId)
		shader.uniforms[var.name] = var

		if var.varType <: SamplerType
			push!(shader.samplers, var.name)
			# set the sampler index once
			# sampler index can only be set via glUniform1i
			glUniform1i(var.location, length(shader.samplers))
			@assert glGetError() == GL_NO_ERROR
		end
	end

	glUseProgram(0)
end

function doneblocks(shader::Shader)
	for block in shader.blocks
		done(block)
	end
	empty!(shader.blocks)
end

function initattributes(shader::Shader)
	@assert isvalid(shader)

	val = GLint[0]

	glGetProgramiv(shader.program, GL_ACTIVE_ATTRIBUTES, val)
	@assert glGetError() == GL_NO_ERROR
	attribCount = val[1]

	glGetProgramiv(shader.program, GL_ACTIVE_ATTRIBUTE_MAX_LENGTH, val)
	@assert glGetError() == GL_NO_ERROR
	nameBuf = Array{UInt8}(undef, val[1])

	typeBuf = GLenum[0]
	sizeBuf = GLint[0]

	typeFields = Array{Expr}(undef, attribCount)
	for i = 1:attribCount
		glGetActiveAttrib(shader.program, i-1, length(nameBuf), C_NULL, sizeBuf, typeBuf, nameBuf)
		@assert glGetError() == GL_NO_ERROR

		location = glGetAttribLocation(shader.program, nameBuf) + 1
		@assert glGetError() == GL_NO_ERROR
		@assert 1 <= location <= attribCount

		# todo: add support for attribute arrays (via fixed size array types)
		if sizeBuf[1] != 1
			error("Vertex attribute arrays aren't supported in shader program $(shader.id), vertex attribute $nameBuf")
		end

		nameSym = Symbol(unsafe_string(pointer(nameBuf)))
		jlType = gl2jltype(typeBuf[1])

		typeFields[location] = :($nameSym::$jlType)
	end

	typeSym = Symbol("Vert#" * join(map(string, typeFields), "#"))
	typeExpr = quote
		struct $typeSym
			$(typeFields...)
		end
		$typeSym
	end

	shader.attribType = eval(typeExpr)
end

function get_sampler_index(shader::Shader, name::Symbol)
	found = findfirst(isequal(name), shader.samplers)
	found == nothing && return -1
	return found
end

function setuniform(shader::Shader, uniform::Symbol, value)
	@assert isvalid(shader)
	@assert gl_get_current_program() == shader.program

	if haskey(shader.uniforms, uniform)
		glUniform(shader.uniforms[uniform].location, value)
		return true
	end

	@info("Attempt to set inexistent uniform $uniform in shader $(shader.name)")
	return false
end

function setuniform(shader::Shader, uniform::Symbol, tex::AbstractTexture)
	@assert isvalid(shader)
	# @assert gl_get_current_program() == shader.program
	@assert isvalid(tex)

	texIndex = get_sampler_index(shader, uniform)
	if texIndex != 0
		@assert haskey(shader.uniforms, uniform)
		# apply the texture state to the texture unit we set on initialization
		apply(tex, texIndex)

		return true
	end

	return false
end


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

export _wsave

_wsave(s, fig::Figure; dpi::Int=200) = fig.savefig(s, bbox_inches="tight", dpi=dpi)


text:
# WAH Elements
# ============
#
# The individual elements of a WAHVector.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/WAHVectors.jl/blob/master/LICENSE.md

# The WAHElement type and constants
# ---------------------------------

primitive type WAHElement 32 end
Base.convert(::Type{WAHElement}, x::UInt32) = reinterpret(WAHElement, x)
Base.convert(::Type{UInt32}, x::WAHElement) = reinterpret(UInt32, x)

WAHElement(x::UInt32) = convert(WAHElement, x)
UInt32(x::WAHElement) = convert(UInt32, x)
function WAHElement(value::UInt32, nruns::UInt32)
    @assert nruns ≤ WAH_MAX_NRUNS "nruns is too big"
    return WAHElement(((0x00000002 + value) << 30) + nruns)
end

const WAH_LITERAL_ZEROS = convert(WAHElement, 0x00000000)
const WAH_LITERAL_ONES = convert(WAHElement, 0x7FFFFFFF)

const WAH_FULL_ZEROS = convert(WAHElement, 0xBFFFFFFF)
const WAH_FULL_ONES = convert(WAHElement, 0xFFFFFFFF)
const WAH_MAX_NRUNS = 0x3FFFFFFF

# Basic general operations on WAHElements
# ---------------------------------------

Base.:(==)(x::WAHElement, y::WAHElement) = convert(UInt32, x) == convert(UInt32, y)

# Basic operations
#-----------------

"""
    isruns(x::WAHElement)

Return `true` if the element of the WAH compressed bit array `x`
represents a compressed run of words. Otherwise return `false`.
"""
isruns(x::WAHElement) = convert(UInt32, x) ≥ 0x80000000

"""
    is_zeros_runs(x::WAHElement)

Return `true` if the WAH Element represents a compressed run of all zero words.
"""
is_zeros_runs(x::WAHElement) = 0xC0000000 > convert(UInt32, x) ≥ 0x80000000

"""
    is_ones_runs(x::WAHElement)

Return `true` if the WAH Element represents a compressed run of all one words.
"""
is_ones_runs(x::WAHElement) = convert(UInt32, x) ≥ 0xC0000000

"""
    isliteral(x::WAHElement)

Return `true` if the element of the WAH compressed bit array `x`
represents a literal word rather than a compressed run of words, in
which case this function will return `false`.
"""
isliteral(x::WAHElement) = convert(UInt32, x) < 0x80000000

"""
    nwords(x::WAHElement)

Get the number of words represented by this single WAHElement.
"""
@inline function nwords(x::WAHElement)
    if isruns(x)
        return nruns(x)
    else
        return 0x00000001
    end
end

"""
    runval(x::WAHElement)

Get the runvalue of a compressed WAHElement.

!!! warning
    Output will not make sense if you use it on a WAH element that
    represents a literal word (i.e. `isliteral(x)` returns `true`) rather
    than a number of compressed words.
"""
runval(x::WAHElement) = (convert(UInt32, x) & 0x40000000) >> 30

"""
    nruns(x::WAHElement)

Get the number of words represented by this single WAHElement, assuming it
is a compressed element.

!!! warning
    Output will not make sense if you use it on a WAH element that
    represents a literal word (i.e. `isliteral(x)` returns `true`) rather
    than a number of compressed words.
"""
nruns(x::WAHElement) = convert(UInt32, x) & WAH_MAX_NRUNS

"""
    nfree(x::WAHElement)

Return how many more words this WAHElement can compress into it.

!!! warning
    Output will not make sense if you use it on a WAH element that
    represents a literal word (i.e. `isliteral(x)` returns `true`) rather
    than a number of compressed words.
"""
nfree(x::WAHElement) = WAH_MAX_NRUNS - nruns(x)

"""
    isfull(x::WAHElement)

Check whether this WAHElement contains as many words as it can.

!!! warning
    Output will not make sense if you use it on a WAH element that
    represents a literal word (i.e. `isliteral(x)` returns `true`) rather
    than a number of compressed words.
"""
isfull(x::WAHElement) = (x == WAH_FULL_ZEROS) || (x == WAH_FULL_ONES)

matchingfills(x::WAHElement, y::WAHElement) = (convert(UInt32, x) >> 30) == (convert(UInt32, y) >> 30)

hasroom(x::WAHElement) = nruns(x) < WAH_MAX_NRUNS
hasroom(x::WAHElement, required::UInt32) = (nruns(x) + required) < WAH_MAX_NRUNS

increment_nruns_unsafe(x::WAHElement) = WAHElement(convert(UInt32, x) + 0x00000001)
increment_nruns_unsafe(x::WAHElement, y::UInt32) = WAHElement(convert(UInt32, x) + y)

@inline function runfill(x::WAHElement)
    return ifelse(isruns(x), ifelse(is_ones_runs(x), 0x7FFFFFFF, 0x00000000), convert(UInt32, x))
end

nbits(x::WAHElement) = UInt64(31) * UInt64(nwords(x))

function Base.show(io::IO, element::WAHElement)
    rf = string(runfill(element), base = 16, pad = 2)
    if isruns(element)
        println(io, "$(nruns(element)) x $rf")
    else
        println(io, "1 x $rf")
    end
end


text:
using ITensors,
      Test

@testset "not" begin
  i = Index(2,"i")
  j = Index(2,"j")
  k = Index(2,"k")

  A = randomITensor(i,j,k)

  Ap = prime(A,not("j"))

  @test hassameinds(Ap,(i',j,k'))

  At = addtags(A,"x",not("k"))

  @test hassameinds(At,(addtags(i,"x"),addtags(j,"x"),k))

  Ap2 = prime(A,2,not(i))

  @test hassameinds(Ap2,(i,j'',k''))

  Ap3 = prime(A,3,not(i,k))

  @test hassameinds(Ap3,(i,j''',k))

  At2 = settags(A,"y",not(IndexSet(j,k)))

  @test hassameinds(At2,(settags(i,"y"),j,k))
end

nothing


text:
abstract type ExponentialFamily end
function bregmanDivergence()
end
function update_parameters!()
end
function dimension()
end
function sourceParameters()
end
function sourceVector()
end
function source_to_expectation()
end
include("WatsonDistribution.jl")
include("MultivariateGaussian.jl")


text:
export AbstractDVIDArray, DVIDBase

abstract AbstractDVIDArray <: AbstractArray

immutable DVIDBase
    address     ::IPAddr
    port        ::Integer
    node       ::String
end


text:
using SparseProbVec
using MTD

Omeg = reshape([0.5, 0.5, 0.4, 0.6, 0.9, 0.1, 0.2, 0.8], (2,2,2))
Omeg = reshape([0.5, 0.5, 0.4, 0.6, 0.9, 0.1, 0.2, 0.8, 0.55, 0.45, 0.4, 0.6, 0.85, 0.15, 0.2, 0.8], fill(2,3+1)...)
Omeg = [0.2, 0.5, 0.3]

Q0 = reduce_intercept(Omeg)
A = deepcopy(Q0[3])

# function reduce_ord2(A::Array{T}) where T <: Real
#     order = ndims(A) - 1
#     K = size(A)[1]
#     order == 2 || error("reduce_ord2 requires three dimensional array.")
#
#     lam_cand = Vector{T}(undef, K) # in generalization, K will be replaced with number of lag combos
#     take_cand = Array{T}(undef, (K, K, order)) # in generalization, first current_order dimensions will be the "take candidates", the remaining will index lag combos
#
#     for ℓ = 1:order # in generalization, we'll somehow iterate over all combinations
#         # Mins = reshape( minimum(A, dims=1+setdiff(1:2, ℓ)[1]), fill(K, order)... ) # in generalization, I think this will need to be an array with one less dimension than A
#         Mins = minimum(A, dims=1+setdiff(1:2, ℓ)[1]) # in generalization, I think this will need to be an array with one less dimension than A
#         sums_Mins = sum(Mins, dims=1)
#         lam_cand[ℓ] = minimum(sums_Mins)
#         take_cand[:,:,ℓ] = reshape( Mins ./ sums_Mins , fill(K, order)...)
#     end
#
#     lam_sel, lag_sel_indx = findmax(lam_cand)
#
#     if abs(lam_sel) < 1.0e-6
#         lam_sel = 0.0
#         leftover = deepcopy(A)
#         take = zeros(T, (K,K))
#         lag_sel_indx = 0
#     else
#         leftover = A .- lam_sel .* take_cand[:,:,lag_sel_indx]
#         take = take_cand[:,:,lag_sel_indx]
#     end
#
#     leftover_sums = vec(sum(leftover, dims=1))
#     maximum(abs.( leftover_sums .- leftover_sums[1] )) < 1.0e-6 || error("All leftover marginal sums should be equal.")
#
#     return lam_sel, lag_sel_indx, take, leftover, leftover_sums[1]
# end

# Q1 = reduce_ord2(A)
# A1 = deepcopy(Q1[4])
# reduce_ord2(A1)

## tests of reduce_gen
Omeg = TransTens(reshape([0.5, 0.5, 0.4, 0.6, 0.9, 0.1, 0.2, 0.8, 0.55, 0.45, 0.4, 0.6, 0.85, 0.15, 0.2, 0.8], fill(2,3+1)...))
Q0 = reduce_intercept(Omeg)
Juno.@enter reduce_gen(Q0[3], 1)
Q1a = reduce_gen(Q0[3], 1)
reduce_gen(Q1a[4], 1)
Juno.@enter reduce_gen(Q1a[4], 2)
Q2a = reduce_gen(Q1a[4], 2)
reduce_gen(Q2a[4], 2)
reduce_gen(Q2a[4], 1)

## Other tests
reduce_transTens(Omeg)
reduce_transTens(Omeg, toler=0.5)

Omeg = reshape([0.8, 0.2, 0.4, 0.6, 0.8, 0.2, 0.4, 0.6], (2,2,2)) # only lag 1 matters
Omeg = reshape([0.7, 0.3, 0.7, 0.3, 0.4, 0.6, 0.4, 0.6], (2,2,2)) # only lag 2 matters
Omeg = reshape([0.7, 0.3, 0.65, 0.35, 0.4, 0.6, 0.4, 0.6], (2,2,2)) # lag 1 matters when lag 2 is first state, doesn't matter when lag 2 is second state

r = reduce_transTens(TransTens(Omeg), leftoversum_tol=1.0e-6)
print(r)

# Juno.@enter reduce_transTens(TransTens(Omeg))

## rDirichlet
K = 2
R = 5
Omeg = reshape(hcat([ rDirichlet(fill(1.0, K)) for i = 1:K^R ]...), fill(K, R+1)...)
out = reduce_transTens(TransTens(Omeg, sum1tol=1.0e-6), leftoversum_tol=1.0e-6)
print(out)
sum(out.levels[:,2])

reduce_transTens(TransTens(out[3][5])) # what we were looking for

## from data
# using BSON: @load, parse
# BSON.parse("/Users/mjheiner/Documents/research/MMTD/MMTDpaper_codeExamples/data/sim1_sets.bson")
# @load "/Users/mjheiner/Documents/research/MMTD/MMTDpaper_codeExamples/data/sim1_sets.bson" y
#
# K = length(unique(y[1]))
# L = 4
# N = zeros(Int, fill(K, L+1)...)
# for i = 1:length(y)
#     N .+= count_trans_L(y[i], K, L)
# end
#
# Omeg_hat = TransTens(transTens_MLE(N))
# out = reduce_transTens(Omeg_hat)
# print(out)

## simulation 1 scenario (not the same as used in paper simulations)
using Random
Random.seed!(11)
L = 4
M = 3
K = 3
λ_indx = build_λ_indx(L, M)
Λ = [0.0, 0.0, 0.0, 1.0]
λ = [ zeros(Float64, λ_indx.lens[m]) for m in 1:M ]
λ[3][3] = 1.0
println(λ_indx.indxs[3][3])
Q0 = zeros(Float64, K)
Q = [ reshape( hcat([ rDirichlet(ones(K)) for k in 1:(K)^m ]...), fill(K, m+1)... ) for m in 1:M ]

Omeg = TransTens(transTensor_mmtd(L, M, K, λ_indx, Λ, λ, Q0, Q))

out = reduce_transTens(Omeg)
print(out)

## simulation 2 scenario (not the same as used in simulations)
Random.seed!(12)
L = 5
M = 5
K = 2
λ_indx = build_λ_indx(L, M)
Λ = [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
λ = [ zeros(Float64, λ_indx.lens[m]) for m in 1:M ]
λ[5][1] = 1.0
println(λ_indx.indxs[5][1])
Q0 = zeros(Float64, K)
Q = [ reshape( hcat([ rDirichlet(ones(K)) for k in 1:(K)^m ]...), fill(K, m+1)... ) for m in 1:M ]
Omeg = TransTens(transTensor_mmtd(L, M, K, λ_indx, Λ, λ, Q0, Q))

out = reduce_transTens(Omeg)
sum(out[1][:,2])

## reduce sum to unity standard in reduce functions
Q2 = reshape([0.32727557, 0.20496972, 0.21370286, 0.25405185, 0.31204655, 0.35311113, 0.17880579, 0.15603652, 0.27311296, 0.27393458, 0.17974408, 0.27320838, 0.29620478, 0.14803316, 0.13574267, 0.42001939, 0.55194678, 0.24776699, 0.11573846, 0.08454777, 0.24961685, 0.32725009, 0.22171401, 0.20141904, 0.29946564, 0.36172226, 0.16753204, 0.17128006, 0.26105888, 0.21806586, 0.22125031, 0.29962494, 0.42724184, 0.27502159, 0.13734542, 0.16039115, 0.26375264, 0.26357236, 0.21962539, 0.25304961, 0.10465893, 0.10479424, 0.10476206, 0.68578476, 0.11517821, 0.11295583, 0.11774202, 0.65412393, 0.28004061, 0.21793946, 0.17852386, 0.32349607, 0.43282710, 0.16505096, 0.26887101, 0.13325093, 0.26559577, 0.35463422, 0.17626712, 0.20350288, 0.12663684, 0.33712198, 0.25674244, 0.27949874],
    fill(4,2+1)...) # pink posterior mean Q2 when p(Lam2) > .5
Q2 = reshape([0.33311281, 0.19747092, 0.21188314, 0.25753312, 0.29564221, 0.37250940, 0.17511323, 0.15673517, 0.27234339, 0.27394530, 0.18228286, 0.27142845, 0.28392956, 0.13990997, 0.12187194, 0.45428853, 0.59164182, 0.22542352, 0.10477350, 0.07816116, 0.24486780, 0.33018902, 0.21656417, 0.20837900, 0.28265837, 0.38362433, 0.16448888, 0.16922841, 0.26011627, 0.21857911, 0.21123421, 0.31007042, 0.46091674, 0.26167019, 0.11828426, 0.15912881, 0.25596878, 0.26842541, 0.21516935, 0.26043645, 0.09823094, 0.08611348, 0.08652709, 0.72912849, 0.10172308, 0.10267278, 0.10688165, 0.68872249, 0.29055474, 0.20131754, 0.18229348, 0.32583424, 0.46557378, 0.16005667, 0.25361717, 0.12075238, 0.26228341, 0.35979494, 0.17587809, 0.20204356, 0.10855796, 0.35744539, 0.25716761, 0.27682904],
    fill(4,2+1)...) # pink posterior mean Q2 when p(Lam2) > .8

sum(Q2, dims=1)
levels, configs, Q, involv = reduce_transTens(TransTens(Q2))
levels
configs
involv
Q[4]

L = 3
R = 2
K = 2
lam_indx = build_λ_indx(L, R)
Lam = [0.3, 0.5, 0.2]
lLam = log.(Lam)
lam = [ [0.6, 0.3, 0.1], [0.0, 0.2, 0.8]]
llam = [log.(lam[i]) for i = 1:length(lam)]
Q0 = [0.6, 0.4]
lQ0 = log.(Q0)
Q = [ [0.5 0.2; 0.5 0.8], reshape([0.8, 0.2, 0.2, 0.8, 0.3, 0.7, 0.1, 0.9], 2,2,2) ]
lQ = [log.(Q[i]) for i = 1:length(Q)]

Omeg = transTensor_mmtd(L, R, K, lam_indx, Lam, lam, Q0, Q, login=false)
Omeg = transTensor_mmtd(L, R, K, lam_indx, lLam, llam, lQ0, lQ, login=true)

dcp = reduce_transTens(TransTens(Omeg))
print(dcp)

max_stars_repo_name: mheiner/MTD

max_issues_repo_name: mheiner/MTD

max_forks_repo_name: mheiner/MTD


text:
import SimpleDiffEq

"""
    ContinuousSystemSetup(resampling = NoResampling, dt::Int = 0.01, Ttr::Int = 0,
        sample_step::Int = 1, n_pts::Int = 100,
        diffeq = (alg = SimpleDiffEq.SimpleATsit5(), abstol = 1.0e-6, reltol = 1.0e-6),
        source, target)

Setup for causality analysis of continuous dynamical systems.

## Mandatory keyword arguments

- **`source`**:  The variable(s) to use as `source` when calling `causality`. Usually integer(s).

- **`target`**:  The variable(s) to use as `target` when calling `causality`. Usually integer(s).

## Optional keywords

- **`resampling`**: An instance of a resampling scheme. Defaults to `NoResampling()`.

- **`dt::Int`**: The time step when iterating the system. Defaults to 0.01.

- **`Ttr::Int`**: The number of transient iterations before starting sampling. Defaults to 0.

- **`sample_step::Int`**: The sampling step in the final orbit. If `sample_step` > 1, then the
    orbit is generated until time `T*sample_step`, after which each `sample_step`th point is
    drawn to generate the final orbit.

- **`diffeq`**: Arguments propagated to init DifferentialEquations.jl. Defaults to
    `diffeq = (alg = SimpleDiffEq.SimpleATsit5(), abstol = 1.0e-6, reltol = 1.0e-6)`.

- **`n_pts::Int`**: The number of points in the orbit. Defaults to 100.

"""
Base.@kwdef struct ContinuousSystemSetup{R} <: DynamicalSystemSetup where R
    """ An instance of a resampling scheme (defaults to `NoResampling()`)"""
    resampling::R = NoResampling()

    """ The time step """
    dt = 0.01

    """ The number of transient iterations (in units of `dt`) before starting sampling. """
    Ttr::Int = 0

    """
    The sampling step. If `sample_step` > 1, then the orbit is generated until time
    T*sample_step, after which each `sample_step`th point is drawn to generate the
    final orbit.
    """
    sample_step::Int = 1

    """ The noise level, expressed as a fraction of the standard deviation of each variable. """
    noise_level = 0

    """ The number of points in the orbit (in units of `dt`). Defaults to 100. """
    n_pts::Int = 100

    """ The variable(s) to use as `source` """
    source

    """ The variable(s) to use as `target` """
    target
end

export ContinuousSystemSetup


text:
export Interval

using LinearAlgebra
using RecipesBase

"""
    Interval <: AbstractShape{1,2}

Fields: `start`, `stop`, `origin`, `a` (length of interval)
"""
struct Interval <: AbstractShape{1,2}
    start::Point{1,Cartesian}
    stop::Point{1,Cartesian}
    origin::Point{1,Cartesian}
    a::Float64
end

"""
    (::Interval)(::Point{1}) -> ::Bool

evaluates to true if the point is between `start` and `stop`.

Alternately, `(::Interval)(x,y,z...)`
"""
(i::Interval)(p::Point) = i.start ≤ p.x ≤ i.stop
(i::Interval)(x::Real,y...) = i(Point(x,y...))

"""
    Interval(start, stop; [reference=:center]) -> interval

1-Dimensional Interval.

`ref` determines where `origin` of the interval is relative to `start` & `stop`
"""
function Interval(start,stop;reference::Symbol=:center)
    p1, p2 = Point(start), Point(stop)
    return Interval(p1,p2,generate_origin(p1,p2,reference),norm(stop-start))
end

function generate_origin(start::Point{1},stop::Point{1},ref::Symbol)
    left,right = minmax(start,stop)
    if Base.sym_in(ref,(:left,:l,:Left,:L))
        origin = Point(left)
    elseif Base.sym_in(ref,(:right,:r,:Right,:R))
        origin = Point(right)
    else
        origin = Point((left+right)/2)
    end
    return origin
end

import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_DARK

function Base.show(io::IO,i::Interval)
    printstyled(io, "Interval",color=PRINTED_COLOR_DARK)
    print(io," (")
    printstyled(io,i.start.x,color=PRINTED_COLOR_NUMBER)
    print(io,", ")
    printstyled(io,i.stop.x,color=PRINTED_COLOR_NUMBER)
    print(io,")")
end

################################################################################
# Plotting

@recipe function f(i::Interval)
    ylims --> (-1,1)
    legend --> false
    [i.start.x,i.stop.x],[0,0]
end


text:
@testset "350.intersection-of-two-arrays-ii.jl" begin
    @test intersection_of_arrays([1, 2, 2, 1], [2, 2]) == [2, 2]
    @test intersection_of_arrays([4, 9, 5], [9, 4, 9, 8, 4]) == [4, 9]
    @test intersection_of_arrays([2], [1]) == Int[]
    @test intersection_of_arrays([2], [2]) == [2]
end


text:
RESET = "\033[0m"
FAINT_COLOR = "\033[90m"
PASS_COLOR = "\033[32m"
PASS_LIGHT_COLOR = "\033[92m"
FAILED_COLOR = "\033[31m"
FAILED_LIGHT_COLOR = "\033[91m"

DESCRIPTION_ERROR_MESSAGE = "PANIC! Got an error while called describe"

HOUR = 3600_000
MINUTES = 60_000
SECOND = 1000

toMilis(n) = int(n * 10e2)

function report(desc::Description)
end

function summaryReport(passes::Array{Test,1},errors::Array{Union(Error,DescriptionError),1},elapsed)
  println('\n',PASS_LIGHT_COLOR, length(passes), " passing ",FAINT_COLOR, "(", showTime(elapsed), ")")
  length(errors) == 0 || println(FAILED_LIGHT_COLOR, length(errors), " failing")
  println(RESET)
  length(errors) == 0  || fullFailedReport(errors)
end

function showTime(n)
  if n > HOUR
    string(n/HOUR)[1:3] * "h"
  elseif n > MINUTES
    string(n/MINUTES)[1:3] * "m"
  elseif n > SECOND
    string(n/SECOND)[1:3] * "s"
  else
    string(n) * "ms"
  end
end

function colorTime(elapsed::Int)
  if elapsed != 0
    color = elapsed > 500 ? 91 : (elapsed > 100 ? 93 : 90)
    " \e[$(color)m($(showTime(elapsed)))$RESET"
  else
    ""
  end
end

function passReport(test::Test, elapsed::Int)
  print(PASS_COLOR,'.', RESET)
end

function failedReport(err::DescriptionError)
  println(FAILED_COLOR,
    "\nPANIC! Got an error while called describe: ",err.desc.name, RESET
  )
end

function failedReport(err::Error)
  print(FAILED_COLOR,'.',RESET)
end

function fullFailedReport(errors::Array{Union(Error,DescriptionError),1})
  for i in 1:length(errors)
    err = errors[i]
    print(PAD, i, ") ")
    if isa(err, DescriptionError)
      print(DESCRIPTION_ERROR_MESSAGE, " - ", err.desc.name)
    else
      print(err.test.desc.name, " - ", err.test.name)
    end
    println(":", FAILED_COLOR)
    dump(err.err)
    println(RESET)
  end
end


text:
using PyCall
using JLD2

load = pyimport("numpy").load

datadir = "/Users/aza/Projects/TTS/data/LJSpeech-1.1"
savepath = joinpath(datadir, "melspectrograms.jld2")
npy_mels_pectrograms_dir = joinpath(datadir, "melspectrograms")

melspectrograms = Dict{String,Matrix{Float32}}()
for (root, _, files) ∈ walkdir(npy_mels_pectrograms_dir)
   for file ∈ files
      key = first(splitext(file))
      npyfilepath = joinpath(root, file)
      value = load(npyfilepath)
      melspectrograms[key] = permutedims(value)
   end
end

@save savepath melspectrograms


text:
using SymEGRSSMatrices
using LinearAlgebra
using Test

import SymEGRSSMatrices: spline_kernel

@testset "SymEGRSSMatrix" begin
    include("test_SymEGRSSMatrix.jl")
end

@testset "SymEGRSSCholesky" begin
    include("test_SymEGRSSCholesky.jl")
end

@testset "SymEGRQSMatrix" begin
    include("test_SymEGRQSMatrix.jl")
end

@testset "SymEGRQSCholesky" begin
    include("test_SymEGRQSCholesky.jl")
end

max_stars_repo_name: mipals/SymEGRSSMatrices


text:
"""
# Usage

    solve(prob::FODESystem, h, GL())

Use Grunwald Letnikov difference method to solve system of system of FODE.

### Reference

```tex
@INPROCEEDINGS{8742063,
author={Clemente-López, D. and Muñoz-Pacheco, J. M. and Félix-Beltrán, O. G. and Volos, C.},
booktitle={2019 8th International Conference on Modern Circuits and Systems Technologies (MOCAST)},
title={Efficient Computation of the Grünwald-Letnikov Method for ARM-Based Implementations of Fractional-Order Chaotic Systems},
year={2019},
doi={10.1109/MOCAST.2019.8742063}}
```

Python version by https://github.com/DClementeL/Grunwald_Letnikov
"""
struct GL <: FractionalDiffEqAlgorithm end
# Grunwald Letnikov discretization method dispatch for FODESystem
# struct GLWithMemory <: FractionalDiffEqAlgorithm end

function solve(prob::FODESystem, h, ::GL)
    @unpack f, α, u0, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    hα = h^α[1]
    n::Int64 = floor(Int64, (T-t0)/h)+1
    l = length(u0)

    # Initialize solution
    result = zeros(Float64, length(u0), n)
    result[:, 1] = u0

    # generating generalized binomial Cα
    Cα = zeros(Float64, n)
    Cα[1] = 1
    @fastmath @inbounds @simd for j in range(2, n, step=1)
        Cα[j] = (1-(1+α[1])/(j-1))*Cα[j-1]
    end

    du = zeros(Float64, l)

    @fastmath @inbounds @simd for k in range(2, n, step=1)
        summation = zeros(Float64, length(u0))

        @fastmath @inbounds @simd for j in range(1, k-1, step=1)
            for i in eachindex(summation)
                summation[i] += Cα[j+1]*result[i, k-j]
            end
        end

        f(du, result[:, k-1], nothing, t0+(k-1)*h)
        result[:, k] .= hα.*du .- summation
    end
    return result'
end


text:
#Import Turing, Distributions and DataFrames
using Turing, Distributions, DataFrames, Distributed

# Import MCMCChain, Plots, and StatsPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots

# Set a seed for reproducibility.
using Random
Random.seed!(12);

# Turn off progress monitor.
Turing.setprogress!(false)

theta_noalcohol_meds = 1    # no alcohol, took medicine
theta_alcohol_meds = 3      # alcohol, took medicine
theta_noalcohol_nomeds = 6  # no alcohol, no medicine
theta_alcohol_nomeds = 36   # alcohol, no medicine

# no of samples for each of the above cases
q = 100

#Generate data from different Poisson distributions
noalcohol_meds = Poisson(theta_noalcohol_meds)
alcohol_meds = Poisson(theta_alcohol_meds)
noalcohol_nomeds = Poisson(theta_noalcohol_nomeds)
alcohol_nomeds = Poisson(theta_alcohol_nomeds)

nsneeze_data = vcat(
    rand(noalcohol_meds, q),
    rand(alcohol_meds, q),
    rand(noalcohol_nomeds, q),
    rand(alcohol_nomeds, q),
)
alcohol_data = vcat(zeros(q), ones(q), zeros(q), ones(q))
meds_data = vcat(zeros(q), zeros(q), ones(q), ones(q))

df = DataFrame(;
    nsneeze=nsneeze_data,
    alcohol_taken=alcohol_data,
    nomeds_taken=meds_data,
    product_alcohol_meds=meds_data .* alcohol_data,
)
df[sample(1:nrow(df), 5; replace=false), :]

#Data Plotting

p1 = Plots.histogram(
    df[(df[:, :alcohol_taken] .== 0) .& (df[:, :nomeds_taken] .== 0), 1];
    title="no_alcohol+meds",
)
p2 = Plots.histogram(
    (df[(df[:, :alcohol_taken] .== 1) .& (df[:, :nomeds_taken] .== 0), 1]);
    title="alcohol+meds",
)
p3 = Plots.histogram(
    (df[(df[:, :alcohol_taken] .== 0) .& (df[:, :nomeds_taken] .== 1), 1]);
    title="no_alcohol+no_meds",
)
p4 = Plots.histogram(
    (df[(df[:, :alcohol_taken] .== 1) .& (df[:, :nomeds_taken] .== 1), 1]);
    title="alcohol+no_meds",
)
plot(p1, p2, p3, p4; layout=(2, 2), legend=false)

# Convert the DataFrame object to matrices.
data = Matrix(df[:, [:alcohol_taken, :nomeds_taken, :product_alcohol_meds]])
data_labels = df[:, :nsneeze]
data

# # Rescale our matrices.
data = (data .- mean(data; dims=1)) ./ std(data; dims=1)

# Bayesian poisson regression (LR)
@model function poisson_regression(x, y, n, σ²)
    b0 ~ Normal(0, σ²)
    b1 ~ Normal(0, σ²)
    b2 ~ Normal(0, σ²)
    b3 ~ Normal(0, σ²)
    for i in 1:n
        theta = b0 + b1 * x[i, 1] + b2 * x[i, 2] + b3 * x[i, 3]
        y[i] ~ Poisson(exp(theta))
    end
end;

# Retrieve the number of observations.
n, _ = size(data)

# Sample using NUTS.

num_chains = 4
m = poisson_regression(data, data_labels, n, 10)
chain = sample(m, NUTS(200, 0.65), MCMCThreads(), 2_500, num_chains; discard_adapt=false)

gelmandiag(chain)

# Taking the first chain
c1 = chain[:, :, 1]

# Calculating the exponentiated means
b0_exp = exp(mean(c1[:b0]))
b1_exp = exp(mean(c1[:b1]))
b2_exp = exp(mean(c1[:b2]))
b3_exp = exp(mean(c1[:b3]))

print("The exponent of the meaned values of the weights (or coefficients are): \n")
println("b0: ", b0_exp)
println("b1: ", b1_exp)
println("b2: ", b2_exp)
println("b3: ", b3_exp)
print("The posterior distributions obtained after sampling can be visualised as :\n")

plot(chain)

# Note the standard deviation before removing the warmup samples
describe(chain)

# Removing the first 200 values of the chains.
chains_new = chain[201:2500, :, :]
describe(chains_new)

plot(chains_new)

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
function molecular_mass(chemical_formula)
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
module Dimacs10

import CodecBzip2
using ..Io: parse_uints
using ...GraphIdx: BiAdjacentIndex

"""
    open_dimacs10(f, fname; as_process = false)

Open file `fname` to perform reading with function `f`.
If `fname` ends with `".bz2"`, start a corresponding decompressor.
If `fname == "-"`, read from `stdin`.
"""
function open_dimacs10(f::Function, fname::String;
                        as_process::Bool=false)
    if endswith(fname, ".bz2")
        if as_process
            open(`bzip2 -d -c $fname`) do io
                return f(io)
            end
        else
            open(fname) do bz2
                io = CodecBzip2.Bzip2DecompressorStream(bz2)
                return f(io)
            end
        end
    else
        if fname == "-"
            f(stdin)
        else
            @assert ispath(fname)
            open(f, fname, "r")
        end
    end
end

"""
    parse_dimacs10_header(io)

Read lines until lines does not start with comment character "%" and
then return the first two integers (number of nodes, number of edges).
"""
function parse_dimacs10_header(io::IO)::Tuple{Int,Int}
    line = ""
    while startswith(line, "%") || line == ""
        line = readline(io)
    end
    n, m = map(x-> parse(Int, x), split(line))
    n, m
end

"""
    parse_dimacs10_idx(io)

Parse a Dimacs10 file into a bidirectional index.
"""
function parse_dimacs10_idx(io::IO)::BiAdjacentIndex
    n, m, head, tail = parse_dimacs10_edges(io, Int)
    return BiAdjacentIndex(n, head, tail)
end

"""
    parse_dimacs10_edges(io)

Return the graph as `n, m, head, tail`.
"""
function parse_dimacs10_edges(io::IO, ::Type{T} = Int32) where {T<:Integer}
    n, m = parse_dimacs10_header(io)
    head = Vector{T}()
    tail = Vector{T}()
    sizehint!(head, m)
    sizehint!(tail, m)
    local i::Ref{T} = Ref{T}(1)
    parse_uints(io) do u::UInt, last::Bool
        if i[] < u
            push!(head, T(i[]))
            push!(tail, T(u))
        end
        if last
            i[] += 1
        end
    end
    return n, m, head, tail
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
@inline space(A::AbstractArray) = Tuple(TrivialVectorSpace(n) for n in size(A))
@inline space(A::AbstractArray, l::Int) = TrivialVectorSpace(size(A, l))

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
# #############################################################################
# #############################################################################
#Julia wrappers for C solvers, and ray solver

#This file contains the Julia wrappers and functions for running the internal solvers written in C (using the GSL library). The following Julia function is provided for the ray solver:

#   raypath,bounceindices,bouncepoints = rayevolve_gsl(bnd::Boundary, idx::RefractiveIndex,init::Array{Float64,1};tmax::Float64=200.0,bouncemax::Int64=500,reltol::Float64=1.0e-12,abstol::Float64=1.0e-12)
#   Computes the trajectory and bounce points in a cavity with boundary <bnd> and index distribution <idx>, of a ray with initial conditions <init>, which is an array of initial (<r>,<theta>,<phi>) or (radial position, angular position, angle of propagation) values or the ray. Runs the ODE solver with relative and absolute tolerances <reltol> and <abstol> (default 1.0e-12 each), and halts computation when pathlength reaches <tmax> (default 200.0) or number of bounces reaches <bouncemax> (default 500).
#   Returns the 2D array <raypath> of (<r>,<theta>) positions traversed in the ray's path as its rows, the 1D array <bouncesindices> of row indices in <raypath> that bounces occur, and the 2D array <bouncepoints> with the (<theta>,<chi>) or (angular position,angle of incidence) values of each bounce as its rows.
#   This means that <raypath>[<bouncesindices>,2] == <bouncepoints>[:,1].

# #############################################################################
# #############################################################################
#Initiate

#Dependencies
#require("basic.jl")
#require("boundary.jl")
#require("refractive_index.jl")

# #############################################################################
# #############################################################################
#Interface for Julia callbacks in C

#Generate C-callable functions wrapping the relevant julia functions for boundary and refractive index types
function rfunc_c(bndptr::Ptr{Void},theta::Cdouble)
    bnd::Boundary = unsafe_pointer_to_objref(bndptr)
    return rfunc(bnd,theta)::Cdouble
end
function rsys!_c(bndptr::Ptr{Void},theta::Cdouble,resultptr::Ptr{Cdouble})
    bnd::Boundary = unsafe_pointer_to_objref(bndptr)
    result::Array{Float64,1} = pointer_to_array(resultptr,2)
    rsys!(bnd,theta,result)
    return nothing
end
function nfunc_c(idxptr::Ptr{Void},r::Cdouble,theta::Cdouble)
    idx::RefractiveIndex = unsafe_pointer_to_objref(idxptr)
    return nfunc(idx,r,theta)::Cdouble
end
function nderiv!_c(idxptr::Ptr{Void},r::Cdouble,theta::Cdouble,resultptr::Ptr{Cdouble})
    idx::RefractiveIndex = unsafe_pointer_to_objref(idxptr)
    result::Array{Float64,1} = pointer_to_array(resultptr,3)
    nderiv!(idx,r,theta,result)
    return nothing
end

#Get C function pointers
const rfunc_p = cfunction(rfunc_c,Cdouble,(Ptr{Void},Cdouble))
const rsys!_p = cfunction(rsys!_c,Void,(Ptr{Void},Cdouble,Ptr{Cdouble}))
const nfunc_p = cfunction(nfunc_c,Cdouble,(Ptr{Void},Cdouble,Cdouble))
const nderiv!_p = cfunction(nderiv!_c,Void,(Ptr{Void},Cdouble,Cdouble,Ptr{Cdouble}))

# #############################################################################
# #############################################################################
#Wrapper for GSL Ray Tracing

#Ray traversal (GSL ODE solver)
#Computes the trajectory and bounce points of a ray
function rayevolve_gsl(bnd::Boundary,idx::RefractiveIndex,init::Array{Float64,1};tmax::Float64=200.0,bouncemax::Int64=500,reltol::Float64=1.0e-12,abstol::Float64=1.0e-12)

    #Convert initial (r,theta,phi) conditions to the vector (r,theta,pr,ptheta)
    # Notation:
    # (r,theta) are polar position coordinates
    # phi is the angle the direction of travel makes with the horizontal
    # (pr,ptheta) are "momentum" coordinates corresponding to (r,theta)
    # chi is the angle of incidence at each bounce
    const n::Float64 = nfunc(idx,init[1],init[2])
    const pr0::Float64 = n*cos(init[3]-init[2])
    const ptheta0::Float64 = n*init[1]*sin(init[3]-init[2])

    #Preallocate results arrays -----------------------------------------------
    #Estimate length required
    const prealloc::Int64 = 250*ceil(Int64,tmax) #Length to preallocate raypath array
    #Initialize arrays for the 1st and 2nd columns of raypath matrix storing
    #(r,theta) positions along the path of the ray as rows
    raypath_r::Array{Float64,1} = Array(Float64,prealloc) #raypath r value
    raypath_theta::Array{Float64,1} = Array(Float64,prealloc) #raypath theta value
    #Initialize arrays for the 1st and 2nd columns of the bounce information matrix
    #storing the corresponding row of the raypath array and angle of incidence (chi)
    #values in each row
    bounceindices::Array{Int64,1} = Array(Int64,bouncemax) #bounce indices
    bouncepts_chi::Array{Float64,1} = Array(Float64,bouncemax) #bounce chi value
    #Record actual number of rows of the raypath and bounces arrays here
    lengths::Array{Int64,1} = Array(Int64,2)

    #Make C call --------------------------------------------------------------
    ccall((:rayevolve,"libcavchaos.so"), Void,
        #Tuple of argument types
        (Ptr{Cdouble},Ptr{Cdouble},Ptr{Clong},Ptr{Cdouble},Ptr{Clong},
         Cdouble,Cdouble,Cdouble,Cdouble,
         Cdouble,Cint,Cdouble,Cdouble,
         Any,Ptr{Void},Ptr{Void},Any,Ptr{Void},Ptr{Void}),
        #Argument values
        raypath_r,raypath_theta,bounceindices,bouncepts_chi,lengths, #Storage arrays
        init[1],init[2],pr0,ptheta0, #Initial conditions
        tmax,bouncemax,reltol,abstol, #Simulation parameters
        bnd,rfunc_p,rsys!_p,idx,nfunc_p,nderiv!_p) #Cavity parameters

    #Truncate or extract data from allocated arrays for output
    #raypath position array
    raypath::Array{Float64,2} = Array(Float64,lengths[1],2)
    for i=1:lengths[1]
        raypath[i,1] = raypath_r[i]
        raypath[i,2] = raypath_theta[i]
    end
    #bounce points
    bouncepoints::Array{Float64,2} = Array(Float64,lengths[2],2)
    for i=1:lengths[2]
        bouncepoints[i,1] = raypath_theta[bounceindices[i]]
        bouncepoints[i,2] = mod(bouncepts_chi[i]+pi,2*pi)-pi
    end
    #indices in raypath corresponding to bounces
    resize!(bounceindices,lengths[2])

    #Return (raypath, bounceindices, bouncepoints), where the rows of raypath are the (r,theta) positions along the ray's path, bounceindices contains the indices of raypath that bounces occur, and the rows of bouncepoints are the (theta,chi) values of each bounce
    return (raypath, bounceindices, bouncepoints)

end


text:
using TimeSeriesClustering
using Test
using JLD2
using Random
# make sure to put Random.seed!() before every evaluation of run_clust, because most of them use random number generators in clustering.jl
# Random.seed!() should give the same result even on different machines:https://discourse.julialang.org/t/is-rand-guaranteed-to-give-the-same-sequence-of-random-numbers-in-every-machine-with-the-same-seed/12344

# put DAM_GER, DAM_CA, CEP GER1, CEP GER18 all years into data.jl file that just tests kmeans with all data inputs

@load normpath(joinpath(dirname(@__FILE__),"reference_generation","run_clust.jld2")) reference_results

Random.seed!(1111)
@testset "run_clust $data" for data in [:CEP_GER1,:CEP_GER18] begin
    ts_input_data = load_timeseries_data(data)
    #mr: method, representation, n_init
    mr = [["kmeans","centroid",1000],
    ["kmeans","medoid",1000],
    ["kmedoids","centroid",1000],
    ["kmedoids","medoid",1000],
    ["hierarchical","centroid",1],
    ["hierarchical","medoid",1]]
    @testset "method=$method + representation=$repr" for (method,repr,n_init) in mr begin
        # somehow the following passes if I use julia runtest.jl, but does not pass if
        # I use ] test TimeSeriesClustering.
        #@testset "default" begin
        #    Random.seed!(1111)
        #        ref = reference_results["$data-$method-$repr-default"]
        #        t = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=n_init)
        #        test_ClustResult(t,ref)
        #    catch
        #        @test reference_results["$data-$method-$repr-default"] == "not defined"
        #    end
        #end
        @testset "n_clust=1" begin
            Random.seed!(1111)
            try
                ref = reference_results["$data-$method-$repr-n_clust1"]
                t = run_clust(ts_input_data;method=method,representation=repr,n_clust=1,n_init=n_init)
                test_ClustResult(t,ref)
            catch
                @test reference_results["$data-$method-$repr-n_clust1"] == "not defined"
            end
        end
        @testset "n_clust=N" begin
            Random.seed!(1111)
            try
                ref = reference_results["$data-$method-$repr-n_clustK"]
                t= run_clust(ts_input_data;method=method,representation=repr,n_clust=ts_input_data.K,n_init=n_init)
                test_ClustResult(t,ref)
            catch
                @test reference_results["$data-$method-$repr-n_clustK"] == "not defined"
            end
        end
      end
    end
  end
end

# Use the same data for all subsequent tests
data = :CEP_GER1
ts_input_data = load_timeseries_data(data)

using Cbc
optimizer = Cbc.Optimizer
method = "kmedoids_exact"
repr = "medoid"
# kmedoids exact: only run for small system because cbc does not solve for large system
# no seed needed because kmedoids exact solves globally optimal
@testset "$method-$repr-$data" begin
    @testset "default" begin
        ref = reference_results["$data-$method-$repr-default"]
        t = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=1,kmexact_optimizer=optimizer)
        test_ClustResult(t,ref)
    end
    @testset "n_clust=1" begin
        ref = reference_results["$data-$method-$repr-n_clust1"]
        t = run_clust(ts_input_data;method=method,representation=repr,n_clust=1,n_init=1,kmexact_optimizer=optimizer)
        test_ClustResult(t,ref)
    end
    @testset "n_clust=1" begin
        ref = reference_results["$data-$method-$repr-n_clustK"]
        t = run_clust(ts_input_data;method=method,representation=repr,n_clust=ts_input_data.K,n_init=1,kmexact_optimizer=optimizer)
        test_ClustResult(t,ref)
    end
end

@testset "MultiClustAtOnce" begin
    method = "hierarchical"
    repr = "centroid"
    Random.seed!(1111)
    ref_array = reference_results["$data-$method-$repr-MultiClust"]
    t_array = run_clust(ts_input_data,[1,5,ts_input_data.K];method=method,representation=repr,n_init=1)
    for i = 1:length(t_array)
        test_ClustResult(t_array[i],ref_array[i])
    end
end

@testset "ClustResultAll" begin
    method = "hierarchical"
    repr = "medoid"
    Random.seed!(1111)
    ref_all = reference_results["$data-$method-$repr-ClustResultAll"]
    t_all = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=10,get_all_clust_results=true)
    test_ClustResult(t_all,ref_all)
    for i = 1:length(t_all.centers_all)
        @test all(t_all.centers_all[i] .≈ ref_all.centers_all[i])
        @test all(t_all.weights_all[i] .≈ ref_all.weights_all[i])
        @test all(t_all.clustids_all[i] .≈ ref_all.clustids_all[i])
        @test all(t_all.cost_all[i] .≈ ref_all.cost_all[i])
        @test all(t_all.iter_all[i] .≈ ref_all.iter_all[i])
    end
end

@testset "AttributeWeighting" begin
    method = "hierarchical"
    repr = "centroid"
    Random.seed!(1111)
    attribute_weights=Dict("solar"=>1.0, "wind"=>2.0, "el_demand"=>3.0)
    ref = reference_results["$data-$method-$repr-AttributeWeighting"]
    t = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=1,attribute_weights=attribute_weights)
    test_ClustResult(t,ref)
end


text:
# already in AAM
## function viewshape(img::Image, lms::Matrix{Float64})
##     imgc, img2 = ImageView.view(img)
##     for i=1:size(lms, 1)
##         annotate!(imgc, img2, AnnotationPoint(lms[i, 2], lms[i, 1], shape='.',
##                                               size=4))
##     end
##     imgc, img2
## end
## viewshape(mat::Matrix{Float64}, lms::Shape) = viewshape(convert(Image, mat), lms)

function nviewall(imgs, padding=10)
    @assert(length(imgs) > 0, "Need at least one image")
    h, w = size(imgs[1])
    n = length(imgs)
    rows = int(floor(sqrt(n)))
    cols = int(ceil(n / rows))
    halfpad = div(padding, 2)
    dat = zeros(rows * (h + padding), cols * (w + padding))
    for i=1:n
        ## wt = W[i, :]
        ## wt = reshape(wt, length(wt))
        ## wim = map_nonzeros(IMSIZE, wt, nzs)
        img = imgs[i]
        img = img ./ (maximum(img) - minimum(img))
        r = div(i - 1, cols) + 1
        c = rem(i - 1, cols) + 1
        dat[(r-1)*(h+padding)+halfpad+1 : r*(h+padding)-halfpad,
            (c-1)*(w+padding)+halfpad+1 : c*(w+padding)-halfpad] = imgs[i]
    end
    nview(dat)
    return dat
end

function normalize(img)
    mn, mx = minimum(img), maximum(img)
    nimg = (img .- mn) ./ (mx - mn)
    nimg
end

nview(img) = ImageView.view(normalize(img))


text:
using stoosc
rc("font", family="arial")

# Dictionary to store all the numerical solutions
solutions = Dict()

# Limit cycle oscillator without noise
model = load_model("model02", "ode")
model = set_parameter_value(model, ["B", "d", "τ", "I"], [10, 2, 7.63, 0])
model = set_initial_conditions(model, [0, 1])
model = set_timespan(model, (0.0, 6.0))
solution = simulate_population(model)
solutions["limit_cycle_ode"] = Dict(
    "t" => solution.time,
    "x" => solution.mean[:, 1]
)

# Noise-induced oscillator without noise
model = set_parameter_value(model, ["B", "d", "τ", "I"], [1, -0.1, 6.2, 0])
solution = simulate_population(model)
solutions["noise_induced_ode"] = Dict(
    "t" => solution.time,
    "x" => solution.mean[:, 1]
)

# Limit cycle oscillator with noise
model = load_model("model02", "sde")
model = set_parameter_value(model, ["B", "d", "τ", "I", "σ"], [10, 2, 7.63, 0, 0.2])
model = set_initial_conditions(model, [0, 1])
model = set_timespan(model, (0.0, 6.0))
solution = simulate_population(model)
solutions["limit_cycle_sde"] = Dict(
    "t" => solution.time,
    "x" => solution.mean[:, 1]
)

# Noise-induced with noise
model = set_parameter_value(model, ["B", "d", "τ", "I",  "σ"], [1, -0.1, 6.2, 0, 0.2])
solution = simulate_population(model)
t = solution.time
x = solution.mean[:, 1]
solutions["noise_induced_sde"] = Dict(
    "t" => solution.time,
    "x" => solution.mean[:, 1]
)

## ==
fig, ax_arr = subplots(1, 2, figsize=(7, 1.5))

ax_arr[1].plot(solutions["limit_cycle_ode"]["t"], solutions["limit_cycle_ode"]["x"], color="black")
ax_arr[1].plot(solutions["limit_cycle_sde"]["t"], solutions["limit_cycle_sde"]["x"], color="blue")
ax_arr[1].set_title("Limit cycle oscillator", loc="left", pad=0)

ax_arr[2].plot(solutions["noise_induced_ode"]["t"], solutions["noise_induced_ode"]["x"], color="black")
ax_arr[2].plot(solutions["noise_induced_sde"]["t"], solutions["noise_induced_sde"]["x"], color="blue")
ax_arr[2].set_title("Noise-induced oscillator", loc="left", pad=0)

for i = 1:2
    ax_arr[i].set_xlabel("Time", labelpad=0)
    ax_arr[i].set_ylabel("X",  labelpad=0)
    ax_arr[i].set_xlim(0, 6)
    ax_arr[i].legend(["σ = 0", "σ = 0.2"], edgecolor="black", framealpha=1.0, loc="best", ncol=2)
end

fig.tight_layout(pad=0.3)
# Save the figure
save_figure(fig, "van_der_pol_traces.svg")

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
using Pkg
Pkg.activate(pwd())

# # Functions

function plus(x,y)
    x + y
end

#+

plus(2, 3)
plus(2, -3)

#+

function plus(x,y)
    return x + y
end

#+

function plus(x, y)
    return x + y
    println("I am a useless line of code!!")
end

#+

plus(4, 5)
plus(3, -5)

#+

function powers(x)
    return x, x^2, x^3, x^4
end

#+

ps = powers(2)
typeof(ps)

#+

x1, x2, x3, x4 = powers(2)
x3

# ### Exercise:
# Write function `power(x::Real, p::Integer)` that for a number $x$ and a (possibly
# negative) integer $p$ computes $x^p$ without using the `^` operator. Use only basic
# arithmetic operators `+`, `-`, `*`, `/` and the `if` condition. The annotation
# `p::Integer` ensures that the input `p` is always an integer.
#
# **Hint:** use recursion.
#
# ---
# ### Solution:

# ---
#
# ## One-line functions

plus(x, y) = x + y

#+

plus(4, 5)
plus(3, -5)

#+

f(φ) = -4sin(φ - π/12)

#+

g(x) = (x -= 1; x *= 2; x)

#+

g(3)

# ### Exercise:
# Write a one-line function that returns `true` if the input argument is an even number and
# `false` otherwise.
#
# **Hint:** use modulo function and [ternary operator](@ref Ternary-operator) `?`.
#
# ---
# ### Solution:

# ---
#
# ## Optional arguments

hello(x = "world") = println("Hello $(x).")

#+

hello()
hello("people")

#+

powers(x, y = x*x, z = y*x, v = z*x) = x, y, z, v

#+

powers(2)
powers(2, 3)

#+

f(x = 1, y = x) = (x, y)
g(x = y, y = 1) = (x, y)

#+

f()
g()

# ### Exercise:
# Write a function which computes the value of the following quadratic form
# $$
# q_{a,b,c}(x,y) = ax^2 + bxy + cy^2,
# $$
# where $a, b, c, x \in \mathbb{R}$. Use optional arguments to set default values for
# parameters
# $$
# a = 1, \quad b = 2a, \quad c = 3(a + b).
# $$
# What is the function value at point $(4, 2)$ for default parameters? What is the
# function value at the same point if we use $c = 3$?
#
# ---
# ### Solution:

# ---
#
# ## Keyword arguments

linear(x; a = 1, b = 0) = a*x + b

#+

linear(2)
linear(2; a = 2)
linear(2; b = 4)
linear(2; a = 2, b = 4)

#+

linear(b = 4, 2, a = 2) # If you use this, you will burn in hell 🔥🔥🔥🔥

#+

a, b = 2, 4

linear(2; a = a, b = b)
linear(2; a, b)

# ### Exercise:
# Write a probability density function for the
# [Gaussian distribution](https://en.wikipedia.org/wiki/Normal_distribution)
# $$
# f_{\mu, \sigma}(x) = \frac{1}{\sigma \sqrt{ 2\pi }} \exp\left\{ -\frac{1}{2} \left(
# \frac{x - \mu}{\sigma} \right) ^2 \right\},
# $$
# where $\mu \in \mathbb{R}$ and $\sigma^2 > 0$. Use keyword arguments to obtain the
# standardized normal distribution ($\mu = 0$ and $\sigma = 1$). Check that the inputs
# are correct.
#
# Bonus: verify that this function is a probability density function, i.e., its integral
# equals 1.
#
# ---
# ### Solution:

# ---
#
# ## Variable number of arguments

nargs(x...) = println("Number of arguments: ", length(x))

#+

nargs()
nargs(1, 2, "a", :b, [1,2,3])

#+

args = (1, 2, 3)

nargs(args[1], args[2], args[3])
nargs(args...)
nargs(args)

#+

nargs(1:100)
nargs(1:100...)

#+

nargs([1,2,3,4,5])
nargs([1,2,3,4,5]...)

#+

roundmod(x, y; kwargs...) = round(mod(x, y); kwargs...)

#+

roundmod(12.529, 5)
roundmod(12.529, 5; digits = 2)
roundmod(12.529, 5; sigdigits = 2)

# ### Exercise:
# Write a function `wrapper`, that accepts a number and applies one of `round`, `ceil` or
# `floor` functions based on the keyword argument `type`. Use the function to solve the
# following tasks:
# - Round `1252.1518` to the nearest larger integer and convert the resulting value
# to `Int64`.
# - Round `1252.1518` to the nearest smaller integer and convert the resulting value
# to `Int16`.
# - Round `1252.1518` to `2` digits after the decimal point.
# - Round `1252.1518` to `3` significant digits.
#
# ---
# ### Solution:

# ---
#
# ## Anonymous functions

h1 = function (x)
    x^2 + 2x - 1
end

#+

h2 = x ->  x^2 + 2x - 1

#+

using Plots

f(x,a) = (x + a)^2
plot(-1:0.01:1, x -> f(x,0.5))

#+

map(x -> x^2 + 2x - 1, [1,3,-1])

#+

map([1,3,-1], [2,4,-2]) do x, y
    println("x = $(x), y = $(y)")
    return x + y
end

#+

function f(x, y)
    println("x = $(x), y = $(y)")
    return x + y
end

map(f, [1,3,-1], [2,4,-2])

# ## Dot syntax for vectorizing functions

x = [0, π/2, 3π/4];
A = zeros(length(x));

for (i, xi) in enumerate(x)
    A[i] = sin(xi)
end

A

#+

A = [sin(xi) for xi in x]

#+

A = sin.(x)

#+

plus(x::Real, y::Real) = x + y

#+

plus(1,3)
plus(1.4,2.7)

#+

x = [1,2,3,4]; # column vector

plus(x, x)
plus.(x, x)

#+

broadcast(plus, x, x)

#+

y = [1 2 3 4]; # row vector

plus.(x, y)
plus.(x, 1)

# ## Function composition and piping

(sqrt ∘ +)(3, 6) # equivalent to sqrt(3 + 6)
(sqrt ∘ abs ∘ sum)([-3, -6, -7])  # equivalent to sqrt(abs(sum([-3, -6, -7])))

#+

[-3, -6, -7] |> sum |> abs |> sqrt
[-4, 9, -16] .|> abs .|> sqrt

#+s

["a", "list", "of", "strings"] .|> [uppercase, reverse, titlecase, length]


text:
using Test, DrWatson

cd()

path = "test project"
name = "lala"

initialize_project(path, force = true)

@test projectname() == path
@test typeof(findproject(@__DIR__)) == String
for p in DrWatson.DEFAULT_PATHS
    @test ispath(joinpath(path, p))
end

@test ispath(projectdir("data"))

@test isfile(joinpath(path, ".gitignore"))
@test isfile(joinpath(path, "README.md"))
@test isfile(joinpath(path, "Project.toml"))

@test_throws ErrorException initialize_project(path, name)

initialize_project(path, name; force = true, authors = ["George", "Nick"])

@test projectname() == name
for p in DrWatson.DEFAULT_PATHS
    @test ispath(joinpath(path, p))
end
@test isfile(joinpath(path, ".gitignore"))
@test isfile(joinpath(path, "README.md"))
@test isfile(joinpath(path, "Project.toml"))
z = read((path*"/Project.toml"), String)
@test occursin("[\"George\", \"Nick\"]", z)

initialize_project(path, name; force = true, authors = "Sophia", git = false)
@test !isdir(joinpath(path, ".git"))
z = read((path*"/Project.toml"), String)
@test occursin("[\"Sophia\"]", z)

cd(path)
@test findproject(pwd()) == pwd()
cd()

rm(path, recursive = true, force = true)
@test !isdir(path)


text:
# Test that Mocking works at various scopes

# Global scope
global_scope() = "foo"

# The @mock macro is essentially a no-op
@test (@mock global_scope()) == global_scope()

# Create a patched version of func() and return the alternative
# version at call sites using the @mock macro
patch = (@patch global_scope() = "bar")
apply(patch) do
    @test (@mock global_scope()) != global_scope()
end

# The @mock macro should return to the original behaviour
@test (@mock global_scope()) == global_scope()

# Local scope within a function
function scope_test()
    function_scope() = "foo"
    @test (@mock function_scope()) == function_scope()

    patch = @patch function_scope() = "bar"
    apply(patch) do
        @test (@mock function_scope()) != function_scope()
    end

    @test (@mock function_scope()) == function_scope()
end
scope_test()

# Local scope within a let-block
let let_scope
    let_scope() = "foo"
    @test (@mock let_scope()) == let_scope()

    patch = @patch let_scope() = "bar"
    apply(patch) do
        @test (@mock let_scope()) != let_scope()
    end

    @test (@mock let_scope()) == let_scope()
end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2021 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

function calculate_transcription_kinetic_limit_array(t::Float64, x::Array{Float64,1},
    parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    kinetic_limit_array = Array{Float64,1}()
    model_parameter_array = parameter_dictionary["model_parameter_array"]
    model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]
    system_array = parameter_dictionary["system_concentration_array"]
    eX = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["transcription_elongation_rate"].parameter_value)       # default units: nt/s
    LX = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["characteristic_transcript_length"].parameter_value)    # default units: nt
    k_cat_characteristic = (eX/LX)

    # helper function -
    r(kcat, L_char, L, polymerase, 𝛕, K, species) = kcat*(L_char/L)*polymerase*(species/(𝛕*K+(1+𝛕)*species))

    # alias the model species -
    gene_venus = x[1]
	mRNA_venus = x[2]
	P_venus = x[3]

    # alias the system species -
    RNAP = system_array[1]
	RIBOSOME = system_array[2]
	P_σ70 = system_array[3]

    # alias the system parameters -
    W_gene_venus = model_parameter_array[model_parameter_index_map[:W_gene_venus]]
	W_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:W_gene_venus_P_σ70]]
	K_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:K_gene_venus_P_σ70]]
	n_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:n_gene_venus_P_σ70]]
	K_gene_venus = model_parameter_array[model_parameter_index_map[:K_gene_venus]]
	𝛕_gene_venus = model_parameter_array[model_parameter_index_map[:𝛕_gene_venus]]
	K_P_venus = model_parameter_array[model_parameter_index_map[:K_P_venus]]
	𝛕_P_venus = model_parameter_array[model_parameter_index_map[:𝛕_P_venus]]
	𝛳_mRNA_venus = model_parameter_array[model_parameter_index_map[:𝛳_mRNA_venus]]
	𝛳_P_venus = model_parameter_array[model_parameter_index_map[:𝛳_P_venus]]

    # compute the transcription kinetic limit array -
    push!(kinetic_limit_array, r(k_cat_characteristic,LX,972,RNAP,𝛕_gene_venus,K_gene_venus,gene_venus))

    # return -
    return kinetic_limit_array
end

function calculate_translation_kinetic_limit_array(t::Float64, x::Array{Float64,1},
    parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    kinetic_limit_array = Array{Float64,1}()
    system_array = parameter_dictionary["system_concentration_array"]
    model_parameter_array = parameter_dictionary["model_parameter_array"]
    model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]
    eL = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["translation_elongation_rate"].parameter_value)         # default units: aa/s
    LL = parse(Float64, parameter_dictionary["biophysical_parameters_dictionary"]["characteristic_protein_length"].parameter_value)       # default units: aa
    k_cat_characteristic = (eL/LL)

    # helper function -
    r(kcat, L_char, L, ribosome, 𝛕, K, species) = kcat*(L_char/L)*ribosome*(species/(𝛕*K+(1+𝛕)*species))

    # alias the model species -
    gene_venus = x[1]
	mRNA_venus = x[2]
	P_venus = x[3]

    # alias the system species -
    RNAP = system_array[1]
	RIBOSOME = system_array[2]
	P_σ70 = system_array[3]

    # alias the system parameters -
    W_gene_venus = model_parameter_array[model_parameter_index_map[:W_gene_venus]]
	W_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:W_gene_venus_P_σ70]]
	K_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:K_gene_venus_P_σ70]]
	n_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:n_gene_venus_P_σ70]]
	K_gene_venus = model_parameter_array[model_parameter_index_map[:K_gene_venus]]
	𝛕_gene_venus = model_parameter_array[model_parameter_index_map[:𝛕_gene_venus]]
	K_P_venus = model_parameter_array[model_parameter_index_map[:K_P_venus]]
	𝛕_P_venus = model_parameter_array[model_parameter_index_map[:𝛕_P_venus]]
	𝛳_mRNA_venus = model_parameter_array[model_parameter_index_map[:𝛳_mRNA_venus]]
	𝛳_P_venus = model_parameter_array[model_parameter_index_map[:𝛳_P_venus]]

    # compute the translation kinetic limit array -
    push!(kinetic_limit_array, r(k_cat_characteristic,LL,320,RIBOSOME,𝛕_P_venus,K_P_venus,mRNA_venus))

    # return -
    return kinetic_limit_array
end

function calculate_dilution_degradation_array(t::Float64, x::Array{Float64,1},
    parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    μ = parameter_dictionary["specific_growth_rate"]
    degradation_dilution_array = Array{Float64,1}()
    model_parameter_array = parameter_dictionary["model_parameter_array"]
    model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]

    # alias the model species -
    gene_venus = x[1]
	mRNA_venus = x[2]
	P_venus = x[3]

    # alias the system parameters -
    W_gene_venus = model_parameter_array[model_parameter_index_map[:W_gene_venus]]
	W_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:W_gene_venus_P_σ70]]
	K_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:K_gene_venus_P_σ70]]
	n_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:n_gene_venus_P_σ70]]
	K_gene_venus = model_parameter_array[model_parameter_index_map[:K_gene_venus]]
	𝛕_gene_venus = model_parameter_array[model_parameter_index_map[:𝛕_gene_venus]]
	K_P_venus = model_parameter_array[model_parameter_index_map[:K_P_venus]]
	𝛕_P_venus = model_parameter_array[model_parameter_index_map[:𝛕_P_venus]]
	𝛳_mRNA_venus = model_parameter_array[model_parameter_index_map[:𝛳_mRNA_venus]]
	𝛳_P_venus = model_parameter_array[model_parameter_index_map[:𝛳_P_venus]]

    # formulate the degradation and dilution terms -
    	push!(degradation_dilution_array, (μ + 𝛳_mRNA_venus)*mRNA_venus)
	push!(degradation_dilution_array, (μ + 𝛳_P_venus)*P_venus)

    # return -
    return degradation_dilution_array
end