text:
using DelimitedFiles

D21=Dict{Int64,Char}(1 =>'-',2 =>'A',3 =>'C',4 =>'D',5 =>'E',6 =>'F',7 =>'G',8 =>'H',9 =>'I',10 =>'K',11 =>'L',12 =>'M',13 =>'N',14=>'P',15=>'Q',16 =>'R',17=>'S',18 =>'T',19 =>'V',20 =>'W',21 =>'Y')
Di21=Dict{Char,Int64}('-'=>1,'X'=>1,'B'=>1,'J'=>1,'U'=>1,'Z'=>1,'O'=>1,'A'=>2,'C'=>3,'D'=>4,'E'=>5,'F'=>6,'G'=>7,'H'=>8,'I'=>9,'K'=>10,'L'=>11,'M'=>12,'N'=>13,'P'=>14,'Q'=>15,'R'=>16,'S'=>17,'T'=>18,'V'=>19,'W'=>20,'Y'=>21)

function readmsanum(infile::String ; format=1, header=false)
	Y = Array{Float64,2}(undef,0,0)
	try
		if header
			Y = readdlm(infile, Int64, skipstart=1)
		else
			Y = readdlm(infile, Int64)
		end
	catch err
		println("inputoutput.jl - readmsanum: readdlm failed on file $infile. The alignment may not be of the correct format.")
		error(err)
	end

	if format==0
		Y .+= 1
	elseif format==1
		if findmin(Y)[1] == 0
			error("inputoutput.jl - readmsanum: file $infile contains a 0 value. `format` should be set to 0.")
		end
	end

	return Y
end

function transform_seq(seq::Any,D::Dict{Char,Int64})
    q=length(seq)
    let_seq=Array{Int64}(undef,q)
    for i in 1:q
       let_seq[i]=D[seq[i]]
    end
  return let_seq
end

#read the MSA as fasta file and transform to a MSA as a matrix of numbers in {1,...,21}.
function translate_fasta_to_num_matrix(msa_fasta::String)
    pfam_sa = readlines(msa_fasta);
    M=Int(length(pfam_sa)/2)
    L=length(pfam_sa[2])
	let_msa=Array{Int64,2}(undef,M,L);
    index=1
	for i in 1:length(pfam_sa)
        if i%2==0
         let_msa[index,:].=transform_seq(pfam_sa[i],Di21)
            index+=1
		end

	end
return let_msa
end

function export_fasta_file(outfile::String,msa::Array{Char,2})
(N,L)=size(msa)
fname=outfile
open(fname,"w") do file
      for i in 1:N
       write(file,">sequence_$i")
       write(file,"\n")
       write(file,String(msa[i,:]))
       write(file,"\n")
       end

end

end

#transform MSA from num to letters
function transform_MSA_fasta(msa::Array{Int64,2})
	(N,L)=size(msa)
	let_msa=Array{Char}(undef,N,L);
	for i in 1:N
		for j in 1:L
		 let_msa[i,j]=D21[msa[i,j]]
        end

	end
return let_msa
end


text:
# Interpolates the x, y points to estimate the value at the point "t".
# val (out) ...... the estimated function value at the point "t".
# t (in).......... the "x" value at which to estimate the value
# X (array, in) .. the x-values of the points
# Y (array, in) .. the y-values of the points
# n (in) ......... the length of the arrays X and Y
# REAL(8), INTENT(in) :: t
# INTEGER, INTENT(in) :: n
# REAL(8), INTENT(in) :: X(n), Y(n)
# REAL(8), INTENT(out) :: val
# REAL(8) :: f, denum
# INTEGER :: i, j
function lagrange_interp(X, Y, xx)
    N = size(X,1)
    retval = 0.0
    for j in 1:N
        f = 1.0
        den = 1.0
        for i in 1:N
            if i == j
                continue
            end
            f = f * ( xx - X[i] )
            den = den * ( X[j] - X[i] )
        end
        retval = retval + Y[j] * f/den
    end
    return retval
end

ext: jl


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

name = "SCOTCH"
version = v"6.1.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://gitlab.inria.fr/scotch/scotch/-/archive/v$(version)/scotch-v$(version).tar.gz","4e54f056199e6c23d46581d448fcfe2285987e5554a0aa527f7931684ef2809e"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/scotch*
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/native_build.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Makefile.patch"
if [[ "${target}" == *apple* || "${target}" == *freebsd* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/OSX_FreeBSD.patch"
fi
if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/Windows.patch"
fi
cd src
make scotch
make esmumps
make prefix=$prefix install
install_license ../LICENSE_en.txt
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=Sys.iswindows)

# The products that we will ensure are always built
products = [
    LibraryProduct("libscotcherr", :libscotcherr),
    LibraryProduct("libscotcherrexit", :libscotcherrexit),
    LibraryProduct("libscotchmetis", :libscotchmetis),
    LibraryProduct("libscotch", :libscotch),
    LibraryProduct("libesmumps", :libesmumps)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9.1.0", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("_add_friends.jl")
include("_add_enemies.jl")
include("_add_exp_score.jl")
include("start_ata.jl")
include("add_bank.jl")
include("_add_irt.jl")
include("add_constraints.jl")
include("add_overlap.jl")
include("psychometrics/load_parameters.jl")
include("compute_estimated_iif.jl")
include("add_obj_fun_maximin_minimax.jl")
include("add_obj_fun_cc_maximin.jl")
include("add_obj_fun_soyster_maximin.jl")
include("add_obj_fun_de_jong_maximin.jl")
include("add_obj_fun_robust_maximin.jl")
include("_group_by_friends.jl")
include("load_design.jl")


text:
using Jets, JetPackTransforms, Test, Wavelets
# warning! This only tests the default Haar wavelet
@testset "dwt" begin
    # 1D transform Haar
    println("1D Haar")
    A = JopDwt(Float64, 128)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 2D transform Haar
    println("2D Haar")
    A = JopDwt(JetSpace(Float64, 128, 128))
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 3D transform HAar
    println("3D Haar")
    A = JopDwt(Float64, 128, 128, 128)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # db1
    wt = wavelet(WT.db1)
    # 1D transform
    println("1D $(wt.name)")
    A = JopDwt(Float64, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 2D transform
    println("2D $(wt.name)")
    A = JopDwt(Float64, 128, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 3D transform
    println("3D $(wt.name)")
    A = JopDwt(Float64, 128, 128, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # db6
    wt = wavelet(WT.db6)
    # 1D transform
    println("1D $(wt.name)")
    A = JopDwt(Float64, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 2D transform
    println("2D $(wt.name)")
    A = JopDwt(Float64, 128, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 3D transform
    println("3D $(wt.name)")
    A = JopDwt(Float64, 128, 128, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # db10
    wt = wavelet(WT.db10)
    # 1D transform
    println("1D $(wt.name)")
    A = JopDwt(Float64, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 2D transform
    println("2D $(wt.name)")
    A = JopDwt(Float64, 128, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)

    # 3D transform
    println("3D $(wt.name)")
    A = JopDwt(Float64, 128, 128, 128; wt=wt)
    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox(lhs, rhs, rtol=1e-7)
end


text:
using Insolation
using Plots
using Dates

using CLIMAParameters
using CLIMAParameters.Planet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# Difference in NH and SH zenith angles at time x in given year
function zdiff(x, year)
    date = xtomarchdate(x,year)
    theta_s, dist = daily_zenith_angle(date, -45., param_set, milankovitch=true)
    theta_n, dist = daily_zenith_angle(date, 45., param_set, milankovitch=true)
    return theta_n - theta_s
end

# x is date relative to March 1, with 1.00 representing March 1 00:00
function xtomarchdate(x, year)
    basedate = Dates.DateTime(year, 3, 1)
    deltat = Dates.Second(round((x-1)*Planet.day(param_set)))
    return basedate + deltat
end

# Earth-Sun distance
function edist(x, year)
    date = xtojandate(x,year)
    _, dist = daily_zenith_angle(date, 0., param_set, milankovitch=true)
    return dist/orbit_semimaj(param_set)
end

# x is date relative to Jan 1, with 1.00 representing Jan 1 00:00
function xtojandate(x, year)
    basedate = Dates.DateTime(year, 1, 1)
    deltat = Dates.Second(round((x-1)*Planet.day(param_set)))
    date = basedate + deltat
    return date
end


text:
"""
    scurve(x, a, b, c)

Sigmoid "S"-shaped curve.

# Arguments
- `x::Real`: input
- `a::Real`: determines the centre level
- `b::Real`: determines the amplitude of the curve
- `c::Real`: determines the steepness or "stepwiseness" of the curve.
             The higher c the sharper the function. A negative c reverses the function.
"""
function scurve(x, a, b, c)
    s = one(x) / (b + exp(-c * (x - a)))
    return s
end

"""
    rainfall_interception_gash(cmax, e_r, canopygapfraction, precipitation, canopystorage, maxevap)

Interception according to the Gash model (for daily timesteps). `cmax` is the maximum canopy storage and
`e_r` is the ratio of the average evaporation from the wet canopy and the average precipitation intensity on
a saturated canopy.
"""
function rainfall_interception_gash(
    cmax,
    e_r,
    canopygapfraction,
    precipitation,
    canopystorage,
    maxevap,
)
    # TODO: add other rainfall interception method (lui)
    # TODO: include subdaily Gash model
    # Hack for stemflow (pt)
    pt = 0.1 * canopygapfraction
    pfrac = max((1.0 - canopygapfraction - pt), 0.0)
    p_sat = (-cmax / e_r) * log(1.0 - min((e_r / pfrac), 1.0))
    p_sat = isinf(p_sat) ? 0.0 : p_sat

    # large storms P > P_sat
    largestorms = precipitation > p_sat

    iwet = largestorms ? (pfrac * p_sat) - cmax : precipitation * pfrac
    isat = largestorms ? (e_r) * (precipitation - p_sat) : 0.0
    idry = largestorms ? cmax : 0.0
    itrunc = 0.0

    stemflow = pt * precipitation

    throughfall = precipitation - iwet - idry - isat - itrunc - stemflow
    interception = iwet + idry + isat + itrunc

    # Now corect for area without any Interception (say open water Cmax -- zero)
    cmaxzero = cmax <= 0.0
    throughfall = cmaxzero ? precipitation : throughfall
    interception = cmaxzero ? 0.0 : interception
    stemflow = cmaxzero ? 0.0 : stemflow

    # Now corect for maximum potential evap
    overestimate = interception > maxevap ? interception - maxevap : 0.0
    interception = min(interception, maxevap)

    # Add surpluss to the thoughdfall
    throughfall = throughfall + overestimate

    return throughfall, interception, stemflow, canopystorage

end

"""
    rainfall_interception_modrut(precipitation, potential_evaporation, canopystorage, canopygapfraction, cmax)

Interception according to a modified Rutter model. The model is solved explicitly and there is no
drainage below `cmax`.
"""
function rainfall_interception_modrut(
    precipitation,
    potential_evaporation,
    canopystorage,
    canopygapfraction,
    cmax,
)

    p = canopygapfraction
    pt = 0.1 * p

    # Amount of p that falls on the canopy
    pfrac = max((1.0 - p - pt), 0.0) * precipitation

    # S cannot be larger than Cmax, no gravity drainage below that
    dd = canopystorage > cmax ? canopystorage - cmax : 0.0
    canopystorage = canopystorage - dd

    # Add the precipitation that falls on the canopy to the store
    canopystorage = canopystorage + pfrac

    # Now do the Evap, make sure the store does not get negative
    dc = -1.0 * min(canopystorage, potential_evaporation)
    canopystorage = canopystorage + dc

    leftover = potential_evaporation + dc
    # Amount of evap not used

    # Now drain the canopy storage again if needed...
    d = canopystorage > cmax ? canopystorage - cmax : 0.0
    canopystorage = canopystorage - d

    # Calculate throughfall
    throughfall = dd + d + p * precipitation
    stemflow = precipitation * pt

    # Calculate interception, this is NET Interception
    netinterception = precipitation - throughfall - stemflow
    interception = -dc

    return netinterception, throughfall, stemflow, leftover, interception, canopystorage

end

"""
    acttransp_unsat_sbm(rootingdepth, ustorelayerdepth, sumlayer, restpotevap, sum_actevapustore, c, usl, θₛ, θᵣ, hb, ust::Bool = false)

Compute actual transpiration for unsaturated zone.
If `ust` is `true`, the whole unsaturated store is available for transpiration.

# Arguments
- `rootingdepth`
- `ustorelayerdepth`
- `sumlayer` (depth (z) of upper boundary unsaturated layer)
- `restpotevap` (remaining evaporation)
- `sum_actevapustore` (cumulative actual transpiration (more than one unsaturated layers))
- `c` (Brooks-Corey coefficient)
- `usl` (thickness of unsaturated zone)
- `θₛ`
- `θᵣ`
- `hb` (air entry pressure)
- `ust`

# Output
- `ustorelayerdepth`
- `sum_actevapustore`
- `restpotevap`
"""
function acttransp_unsat_sbm(
    rootingdepth,
    ustorelayerdepth,
    sumlayer,
    restpotevap,
    sum_actevapustore,
    c,
    usl,
    θₛ,
    θᵣ,
    hb,
    ust::Bool = false,
)

    # AvailCap is fraction of unsat zone containing roots
    if ust
        availcap = ustorelayerdepth * 0.99
    else
        if usl > 0.0
            availcap = min(1.0, max(0.0, (rootingdepth - sumlayer) / usl))
        else
            availcap = 0.0
        end
    end

    maxextr = availcap * ustorelayerdepth

    # Next step is to make use of the Feddes curve in order to decrease actevapustore when soil moisture values
    # occur above or below ideal plant growing conditions (see also Feddes et al., 1978). h1-h4 values are
    # actually negative, but all values are made positive for simplicity.
    h1 = hb  # cm (air entry pressure)
    h2 = 100.0  # cm (pF 2 for field capacity)
    h3 = 400.0  # cm (pF 3, critical pF value)
    h4 = 15849.0  # cm (pF 4.2, wilting point)

    # According to Brooks-Corey
    par_lambda = 2.0 / (c - 3.0)
    if usl > 0.0
        vwc = ustorelayerdepth / usl
    else
        vwc = 0.0
    end
    vwc = max(vwc, 0.0000001)
    head = hb / (pow(((vwc) / (θₛ - θᵣ)), (1.0 / par_lambda)))  # Note that in the original formula, thetaR is extracted from vwc, but thetaR is not part of the numerical vwc calculation
    head = max(head, hb)

    # Transform h to a reduction coefficient value according to Feddes et al. (1978).
    # For now: no reduction for head < h2 until following improvement (todo):
    #       - reduction only applied to crops
    if head <= h1
        alpha = 1.0
    elseif head >= h4
        alpha = 0.0
    elseif (head < h2) && (head > h1)
        alpha = 1.0
    elseif (head > h3) && (head < h4)
        alpha = 1.0 - (head - h3) / (h4 - h3)
    else
        alpha = 1.0
    end

    actevapustore = (min(maxextr, restpotevap, ustorelayerdepth)) * alpha
    ustorelayerdepth = ustorelayerdepth - actevapustore
    restpotevap = restpotevap - actevapustore
    sum_actevapustore = actevapustore + sum_actevapustore

    return ustorelayerdepth, sum_actevapustore, restpotevap
end

"""
    infiltration(avail_forinfilt, pathfrac, cf_soil, tsoil, infiltcapsoil, infiltcappath, ustorecapacity, modelsnow::Bool, soilinfreduction::Bool)

Soil infiltration based on infiltration capacity soil `infiltcapsoil`, infiltration capacity compacted area
`infiltcappath` and capacity unsatured zone `ustorecapacity`. The soil infiltration capacity can be adjusted
in case the soil is frozen (`modelsnow` and `soilinfreduction` is `true`).

"""
function infiltration(
    avail_forinfilt,
    pathfrac,
    cf_soil,
    tsoil,
    infiltcapsoil,
    infiltcappath,
    ustorecapacity,
    modelsnow::Bool,
    soilinfreduction::Bool,
)
    # First determine if the soil infiltration capacity can deal with the amount of water
    # split between infiltration in undisturbed soil and compacted areas (paths)
    soilinf = avail_forinfilt * (1.0 - pathfrac)
    pathinf = avail_forinfilt * pathfrac
    if modelsnow && soilinfreduction
        bb = 1.0 / (1.0 - cf_soil)
        soilinfredu = scurve(tsoil, Float(0.0), bb, Float(8.0)) + cf_soil
    else
        soilinfredu = 1.0
    end
    max_infiltsoil = min(infiltcapsoil * soilinfredu, soilinf)
    max_infiltpath = min(infiltcappath * soilinfredu, pathinf)
    infiltsoilpath = min(max_infiltpath + max_infiltsoil, max(0.0, ustorecapacity))

    if max_infiltpath + max_infiltsoil > 0.0
        infiltsoil =
            max_infiltsoil *
            min(1.0, max(0.0, ustorecapacity) / (max_infiltpath + max_infiltsoil))
        infiltpath =
            max_infiltpath *
            min(1.0, max(0.0, ustorecapacity) / (max_infiltpath + max_infiltsoil))
    else
        infiltsoil = 0.0
        infiltpath = 0.0
    end

    infiltexcess = (soilinf - max_infiltsoil) + (pathinf - max_infiltpath)

    return infiltsoilpath, infiltsoil, infiltpath, soilinf, pathinf, infiltexcess
end

"""
    unsatzone_flow_layer(usd, kv_z, l_sat, c)

Assuming a unit head gradient, the transfer of water from an unsaturated store layer `usd` is controlled by the
vertical saturated hydraulic conductivity `kv_z` (bottom layer or water table), the effective saturation
degree of the layer (ratio `usd` and `l_sat`), and a Brooks-Corey power coefficient `c`.
"""
function unsatzone_flow_layer(usd, kv_z, l_sat, c)
    if usd <= 0.0
        return 0.0, 0.0
    end
    sum_ast = 0.0
    # first transfer soil water > maximum soil water capacity layer (iteration is not required because of steady theta (usd))
    st = kv_z * min(pow(usd / l_sat, c), 1.0)
    st_sat = max(0.0, usd - l_sat)
    usd -= min(st, st_sat)
    sum_ast = sum_ast + min(st, st_sat)
    ast = max(min(st - min(st, st_sat), usd), 0.0)
    # number of iterations (to reduce "overshooting") based on fixed maximum change in soil water per iteration step (0.2 mm / model timestep)
    its = Int(cld(ast, 0.2))
    for _ = 1:its
        st = (kv_z / its) * min(pow(usd / l_sat, c), 1.0)
        ast = min(st, usd)
        usd -= ast
        sum_ast += ast
    end

    return usd, sum_ast
end

"""
    unsatzone_flow_sbm(ustorelayerdepth, soilwatercapacity, satwaterdepth, kv_z, usl, θₛ, θᵣ)

The transfer of water from the unsaturated store `ustorelayerdepth` to the saturated store `satwaterdepth`
is controlled by the vertical saturated hydraulic conductivity `kv_z` at the water table and the ratio between
`ustorelayerdepth` and the saturation deficit (`soilwatercapacity` minus `satwaterdepth`). This is the
original Topog_SBM vertical transfer formulation.

"""
function unsatzone_flow_sbm(
    ustorelayerdepth,
    soilwatercapacity,
    satwaterdepth,
    kv_z,
    usl,
    θₛ,
    θᵣ,
)

    sd = soilwatercapacity - satwaterdepth
    if sd <= 0.00001
        ast = 0.0
    else
        st = kv_z * min(ustorelayerdepth, usl * (θₛ - θᵣ)) / sd
        ast = min(st, ustorelayerdepth)
        ustorelayerdepth = ustorelayerdepth - ast
    end

    return ustorelayerdepth, ast

end

"""
    snowpack_hbv(snow, snowwater, precipitation, temperature, tti, tt, ttm, cfmax, whc)

HBV type snowpack modeling using a temperature degree factor.
All correction factors (RFCF and SFCF) are set to 1.
The refreezing efficiency factor is set to 0.05.

# Arguments
- `snow` (snow storage)
- `snowwater` (liquid water content in the snow pack)
- `precipitation` (throughfall + stemflow)
- `temperature`
- `tti` (snowfall threshold interval length)
- `tt` (threshold temperature for snowfall)
- `ttm` (melting threshold)
- `cfmax` (degree day factor, rate of snowmelt)
- `whc` (Water holding capacity of snow)
- `rfcf` correction factor for rainfall
- `sfcf` correction factor for snowfall
- `cfr` refreeing efficiency constant in refreezing of freewater in snow

# Output
- `snow`
- `snowwater`
- `snowmelt`
- `rainfall` (precipitation that occurs as rainfall)
- `snowfall` (precipitation that occurs as snowfall)
"""
function snowpack_hbv(
    snow,
    snowwater,
    precipitation,
    temperature,
    tti,
    tt,
    ttm,
    cfmax,
    whc;
    rfcf = 1.0,
    sfcf = 1.0,
    cfr = 0.05,
)

    # fraction of precipitation which falls as rain
    rainfrac = if iszero(tti)
        Float(temperature > tt)
    else
        frac = (temperature - (tt - tti / 2.0)) / tti
        min(frac, 1.0)
    end
    rainfrac = max(rainfrac, 0.0)

    # fraction of precipitation which falls as snow
    snowfrac = 1.0 - rainfrac
    # different correction for rainfall and snowfall
    snowfall = snowfrac * sfcf * precipitation  # snowfall depth
    rainfall = rainfrac * rfcf * precipitation  # rainfall depth
    # potential snow melt, based on temperature
    potsnowmelt = temperature > ttm ? cfmax * (temperature - ttm) : 0.0
    # potential refreezing, based on temperature
    potrefreezing = temperature < ttm ? cfmax * cfr * (ttm - temperature) : 0.0
    # actual refreezing
    refreezing = temperature < ttm ? min(potrefreezing, snowwater) : 0.0

    # no landuse correction here
    snowmelt = min(potsnowmelt, snow)  # actual snow melt
    snow = snow + snowfall + refreezing - snowmelt  # dry snow content
    snowwater = snowwater - refreezing  # free water content in snow
    maxsnowwater = snow * whc  # max water in the snow
    snowwater = snowwater + snowmelt + rainfall  # add all water and potentially supersaturate the snowpack
    rainfall = max(snowwater - maxsnowwater, 0.0)  # rain + surpluss snowwater
    snowwater = snowwater - rainfall

    return snow, snowwater, snowmelt, rainfall, snowfall
end

"""
    glacier_hbv(glacierfrac, glacierstore, snow, temperature, tt, cfmax, g_sifrac, Δt)

HBV-light type of glacier modelling.
First, a fraction of the snowpack is converted into ice using the HBV-light
model (fraction between 0.001-0.005 per day).
Glacier melting is modelled using a temperature degree factor and only
occurs if the snow cover < 10 mm.

# Arguments
- `glacierFrac` fraction covered by glaciers [-]
- `glacierstore` volume of the glacier [mm] w.e.
- `snow` snow pack on top of glacier [mm]
- `temperature` air temperature [°C]
- `tt` temperature threshold for ice melting [°C]
- `cfmax` ice degree-day factor in [mm/(°C/day)]
- `g_sifrac` fraction of the snow turned into ice [-]
- `Δt` model timestep [s]

# Output
- `snow`
- `snow2glacier`
- `glacierstore`
- `glaciermelt`

"""
function glacier_hbv(glacierfrac, glacierstore, snow, temperature, tt, cfmax, g_sifrac, Δt)

    # Fraction of the snow transformed into ice (HBV-light model)
    snow2glacier = g_sifrac * snow
    snow2glacier = glacierfrac > 0.0 ? snow2glacier : 0.0

    # Max conversion to 8mm/day
    snow2glacier = min(snow2glacier, 8.0 * (Δt / basetimestep))

    snow = snow - (snow2glacier * glacierfrac)
    glacierstore = glacierstore + snow2glacier

    # Potential snow melt, based on temperature
    potmelt = temperature > tt ? cfmax * (temperature - tt) : 0.0

    # actual Glacier melt
    glaciermelt = snow < 10.0 ? min(potmelt, glacierstore) : 0.0
    glacierstore = glacierstore - glaciermelt

    return snow, snow2glacier, glacierstore, glaciermelt

end


text:
using TOML
export PDSpan, PDRelation
export readpdfanno

mutable struct PDSpan
    page::Int
    range::Range
    label::String
    text::String
end

mutable struct PDRelation
    head::String
    tail::String
    label::String
end

mutable struct PDFAnno
    spans::Vector{PDSpan}
end

function readpdfanno(filename::String)
    dict = TOML.parsefile(filename)
    if haskey(dict, "spans")
        spans = map(dict["spans"]) do d
            r = d["textrange"]
            PDSpan(d["page"], r[1]:r[2], d["label"], d["text"])
        end
    else
        spans = PDSpan[]
    end
    PDFAnno(spans)
end


text:
module Interp
using PyCall
export meshgrid, spline_interp2, nearest_interp2, linear_interp1

@pyimport scipy.interpolate as scii

function meshgrid(side_x,side_y)
	x = repmat(reshape([side_x],(1,length(side_x))) ,length(side_y),1)
	y = repmat(reshape([side_y],(length(side_y),1)) ,1,length(side_x))
	x,y
end

function  perodic_padd_xy(x,y,pad_proportion = 0.05)
	maxx = x[1,end];
	minx = x[1,1];
	maxy = y[end,1];
	miny = y[1,1];
	Deltx = x[2,2]-x[1,1];
	Delty = y[2,1]-y[1,1];

	npad = round(pad_proportion*size(x,1));
	middle_x_section =    [x[:,(end-npad+1):end] .- maxx .+ minx .- Deltx x x[:,1:npad] .- minx .+ maxx .+ Deltx];
	full_x = [middle_x_section[1:npad,:]; middle_x_section; middle_x_section[1:npad,:] ];

	middle_y_section = [ y[(end-npad+1):end,:] .- maxy .+ miny .- Delty ; y ; y[1:npad,:] .- miny .+ maxy .+ Delty ];
	full_y =  [middle_y_section[:,1:npad] middle_y_section middle_y_section[:,1:npad] ];

	full_x, full_y
end

function  perodic_padd_z(z,pad_proportion = 0.05)
	npad = int(pad_proportion*size(z,1));
	z_right = z[:, (end-npad+1):end]
	z_left = z[:, 1:npad];
	z_top = z[1:npad,:];
	z_bottom = z[(end-npad+1):end,:];
	z_top_left = z[1:npad, 1:npad];
	z_bottom_left = z[(end-npad+1):end,1:npad];
	z_top_right = z[1:npad, (end-npad+1):end];
	z_bottom_right = z[(end-npad+1):end, (end-npad+1):end];
	full_z = [z_bottom_right z_bottom z_bottom_left;z_right z z_left;z_top_right z_top z_top_left];
	full_z
end
function  perodic_padd(x,y,z,pad_proportion = 0.05)
	full_x,full_y = perodic_padd_xy(x,y,pad_proportion)
	full_z = perodic_padd_z(z,pad_proportion)
	full_x,full_y,full_z
end

function  spline_interp2(xtmp,ytmp,ztmp,xi,yi, pad_proportion = 0.02)
	x,y,z = perodic_padd(xtmp,ytmp,ztmp,pad_proportion)
	iz = scii.RectBivariateSpline(x[1,:],y[:,1],z,s=0)[:ev];
	#iz = scii.RectBivariateSpline(x[1,:],y[:,1],z)[:ev];
	zi = iz(yi[:],xi[:]);
	return reshape(zi,size(xi))
end

function  nearest_interp2(xtmp,ytmp,ztmp,xi,yi; pad_proportion = 0.05)
	x,y,z = perodic_padd(xtmp,ytmp,ztmp,pad_proportion)
	xmin = x[1,1];
	ymin = x[1,1];
	xdelt = x[2,2]-x[1,1];
	ydelt = y[2,1]-y[1,1];

	zi = Array(Float64,size(xi));
	ind_x =  int(round((xi .- xmin)./xdelt .+ 1));
	ind_y =  int(round((yi .- ymin)./ydelt .+ 1));
	for j=1:size(xi,2)
		for i=1:size(xi,1)
			zi[i,j] = z[ind_y[i,j], ind_x[i,j]];
		end
	end
	return zi
end

 function linear_interp1(x,z,xi)
	# assumes that x is a monotonic grid
	xitmp = xi[:]
	deltx = x[2] - x[1]
	xi_upindex = int( ceil(xitmp/deltx) .- x[1]/deltx .+ 1 )
	xi_downindex = xi_upindex .- 1
	idown = find(xi_downindex .< 1)
	xi_upindex[idown] = 2
	xi_downindex[idown] = 1
	iup = find(xi_upindex .> endof(x))
	xi_upindex[iup]=endof(x)
	xi_downindex[iup]=endof(x)-1

	slope=(z[xi_upindex] - z[xi_downindex]) ./ (x[xi_upindex] - x[xi_downindex])
	answer = slope .* (xitmp .- x[xi_downindex]) .+ z[xi_downindex]
	reshape(answer,size(xi))
end

end

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
using PGFPlots

using CSV
using DataFrames
using Missings
using Query
using DocOpt

doc = """

Usage:

"""

args = docopt(doc)

function only(a)
    @assert length(a) == 1
    return first(a)
end

# filename = Pkg.dir("Multilane", "data/all_gaps_Friday_13_Apr_10_17.csv")
filename = Pkg.dir("Multilane", "data/all_gaps_Monday_16_Apr_15_05.csv")

data = CSV.read(filename)

data[:time_to_lane] = data[:steps_to_lane].*data[:dt]
data[:brakes_per_step] = data[:nb_brakes]./data[:n_steps]

combined = by(data, [:cor, :solver, :lambda]) do df
    n = size(df, 1)
    return DataFrame(reached_lane=sum(df[:terminal].=="lane")/n,
                     had_brakes=sum(df[:nb_brakes].>=1)/n,
                    )
end

for cordf in groupby(combined, [:cor])
    @show sort(cordf, cols=[:solver])
    cor = only(unique(cordf[:cor]))
    plts = Plots.Plot[]
    for df in groupby(cordf, [:solver])
        df = sort(df, cols=[:lambda])
        name = only(unique(df[:solver]))
        plt = Plots.Linear(1.0.-df[:had_brakes], df[:reached_lane], legendentry=name)
        push!(plts, plt)
    end
    a = Axis(plts,
             title="Correlation $cor",
             ylabel="Fraction Successful",
             xlabel="Fraction Safe",
             legendPos="south west"
            )

    file = tempname()*".pdf"
    save(file, a)
    run(`xdg-open $file`)
end


text:
using Test, MLBStatsAPI

@testset "game" begin
    gamepk = 566279
    params = Dict(
        "gamePk" => gamepk
    )
    result = game(params)
    @test result["gamePk"] == gamepk
end

@testset "schedulemlb" begin
    params = Dict(
        "startDate" => "2018-07-01",
        "endDate" => "2018-07-02",
        "sportId" => 1
    )
    result = schedulemlb(params)
    @test result["totalGames"] == 25
end

@testset "attendance" begin
    params = Dict(
        "teamId" => 100
    )
    result = attendance(params)
    records = result["records"]
    for record in records
        @test record["team"]["name"] == "Georgia Tech Yellow Jackets"
    end
end

@testset "awards" begin
    params = Dict("awardId" => "MLBHOF", "season" => 2019)
    result = awards(params)
    recps = result["awards"]
    for recp in recps
        @test recp["name"] == "Hall Of Fame"
    end
end

@testset "draft" begin
    params = Dict("year" => 2019)
    result = draft(params)
    queriedsize = result["queriedSize"]
    @test queriedsize == 5134
end

@testset "people" begin
    params = Dict("personIds" => 661255)
    result = people(params)
    ps = result["people"]
    for p in ps
        @test p["nameFirstLast"] == "Yoan Lopez"
    end
end

@testset "standings" begin
    params = Dict("leagueId" => 103, "season" => 2019)
    result = standings(params)
    records = result["records"]
    for record in records
        leagueid = record["league"]["id"]
        @test leagueid == 103
    end
end

@testset "teams" begin
    params = Dict("season" => 2019)
    result = teams(params)
    tms = result["teams"]
    @test length(tms) == 620
end

@testset "meta" begin
    params = Dict("type" => "awards")
    result = meta(params)
    @test haskey(result, "awards")
end


text:
"""
Function to simulate a galaxy depending on type of parameters specified given the sample size.
"""
function simulate_galaxy(p::NFWParameters{G}, samplesize::T; x0::G = 10.0^-8, αt::G = 1.0, rate::Bool = false) where {G <: AbstractFloat, T <: Integer}
    function nfwopt(x::Array{G, 1})
        return -(x[1] * x[1] * x[3] * profile_density(x[1], x[2], x[3], p, αt = αt))
    end
    #function nfwopt(x::Array{G, 1})
    #    return -(2.0 * log(x[1]) +  log(x[3]) + log_profile_density(x[1], x[2], x[3], p, αt = αt))
    #end

    minProb = Optim.optimize(nfwopt, [0.1, 1.0, 1.0], method = Optim.NelderMead())
    fmax = max(0.0,  -1.1 * minProb.minimum)

    vmax0 = escape_velocity(x0, p)
    sampledr = Array{Float64}(samplesize)
    sampledvr = Array{Float64}(samplesize)
    sampledvt = Array{Float64}(samplesize)
    if rate
        tested = 0
    end
    accepted = 0
    while accepted < samplesize
        x = (p.rlim  / p.rs) * rand()
        vr = vmax0 * rand()
        vt = vmax0 * rand()
        u = rand()
        #u = log(rand())
        sampleDensity = x * x * vt * profile_density(x, vr, vt, p, αt = αt)
        #sampleDensity = 2.0 * log(x) + log(vt) + log_profile_density(x, vr, vt, p, αt = αt)
        if rate
            tested += 1
        end
        if sampleDensity >= (u * fmax)
            #if sampleDensity >= (u + fmax)
            accepted += 1
            sampledr[accepted] = x * p.rs
            sampledvr[accepted] =  vr
            sampledvt[accepted] = vt
        end
    end
    if rate
        return SphericalGalaxy(sampledr, sampledvr, sampledvt, samplesize), accepted / tested
    end
    return SphericalGalaxy(sampledr, sampledvr, sampledvt, samplesize)
end

function simulate_galaxy(p::SFWParameters{G}, samplesize::T; x0::G = 10.0^-8, αt::G = 1.0, rate::Bool = false) where {G <: AbstractFloat, T <: Integer}
    function sfwopt(x::Array{G, 1})
        return -(x[1] * x[1] * x[3] * profile_density(x[1], x[2], x[3], p, αt = αt))
    end

    minProb = Optim.optimize(sfwopt, [0.1, 1.0, 1.0], method = Optim.NelderMead())
    fmax = max(0.0,  -1.1 * minProb.minimum)

    vmax0 = escape_velocity(x0, p)
    sampledr = Array{Float64}(samplesize)
    sampledvr = Array{Float64}(samplesize)
    sampledvt = Array{Float64}(samplesize)
    if rate
        tested = 0
    end
    accepted = 0
    while accepted < samplesize
        x = (p.rlim  / p.rs) * rand()
        vr = vmax0 * rand()
        vt = vmax0 * rand()
        u = rand()
        sampleDensity = x * x * vt * profile_density(x, vr, vt, p, αt = αt)
        if rate
            tested += 1
        end
        if sampleDensity >= (u * fmax)
            accepted += 1
            sampledr[accepted] = x * p.rs
            sampledvr[accepted] =  vr
            sampledvt[accepted] = vt
        end
    end
    if rate
        return SphericalGalaxy(sampledr, sampledvr, sampledvt, samplesize), accepted / tested
    end
    return SphericalGalaxy(sampledr, sampledvr, sampledvt, samplesize)
end

function simulate_galaxy(p::NFWParameters{G}, mμ::G, mσ::G, samplesize::T; x0::G = 10.0^-8, αt::G = 1.0, rate::Bool = false) where {G <: AbstractFloat, T <: Integer}
    metal = randn(samplesize) .* mσ .+ mμ
    if rate
        sg, acptrt = simulate_galaxy(p, samplesize, x0 = x0, αt = αt, rate = rate)
        return MetallicSphericalGalaxy(sg, metal), acptrt
    else
        return MetallicSphericalGalaxy(simulate_galaxy(p, samplesize, x0 = x0, αt = αt, rate = rate), metal)
    end
end

function simulate_galaxy(p::SFWParameters{G}, mμ::G, mσ::G, samplesize::T; x0::G = 10.0^-8, αt::G = 1.0, rate::Bool = false) where {G <: AbstractFloat, T <: Integer}
    metal = randn(samplesize) .* mσ .+ mμ
    if rate
        sg, acptrt = simulate_galaxy(p, samplesize, x0 = x0, αt = αt, rate = rate)
        return MetallicSphericalGalaxy(sg, metal), acptrt
    else
        return MetallicSphericalGalaxy(simulate_galaxy(p, samplesize, x0 = x0, αt = αt, rate = rate), metal)
    end
end

"""
Use a grid rejection sampler instead of a standard rejection sampler
"""
function simulate_galaxy_grid(p::SFWParameters{G}, samplesize::T; x0::G = 10.0^-8, αt::G = 1.0, rate::Bool = false, ngrid::Integer = 3) where {G <: AbstractFloat, T <: Integer}

    vmax0 = escape_velocity(x0, p)

    function sfwopt(x::Array{G, 1})
        return -(x[1] * x[1] * x[3] * profile_density(x[1], x[2], x[3], p, αt = αt))
    end

    ##grid bounds
    xgrid = linspace(0.0, p.rlim  / p.rs, ngrid + 1)
    vrgrid = linspace(0.0, vmax0, ngrid + 1)
    vtgrid = linspace(0.0, vmax0, ngrid + 1)

    ##midpoints
    xmids = xgrid[1:ngrid] + 0.5 * step(xgrid)
    vrmids = vrgrid[1:ngrid] + 0.5 * step(vrgrid)
    vtmids = vtgrid[1:ngrid] + 0.5 * step(vtgrid)

    ##optimization variables
    startpoints = Array{Float64}(ngrid^3, 3)
    lowers = Array{Float64}(ngrid^3, 3)
    uppers = Array{Float64}(ngrid^3, 3)
    fmaxes = Array{Float64}(ngrid^3)
    ii = 1
    for xx in 1:ngrid, rr in 1:ngrid, tt in 1:ngrid
        startpoints[ii, :] = [xmids[xx], vrmids[rr], vtmids[tt]]
        lowers[ii, :] = [xgrid[xx], vrgrid[rr], vtgrid[tt]]
        uppers[ii, :] = [xgrid[xx + 1], vrgrid[rr + 1], vtgrid[tt + 1]]
        minProb = Optim.optimize(sfwopt, startpoints[ii, :], lowers[ii, :], uppers[ii, :], Fminbox{Optim.NelderMead}())
        fmaxes[ii] = max(0.0,  -1.1 * minProb.minimum)
        ii += 1
    end
    w = StatsBase.Weights(fmaxes)
    #http://julianlsolvers.github.io/Optim.jl/latest/user/minimization/#box-minimization

    sampledr = Array{Float64}(samplesize)
    sampledvr = Array{Float64}(samplesize)
    sampledvt = Array{Float64}(samplesize)
    if rate
        tested = 0
    end
    accepted = 0
    while accepted < samplesize
        idx = StatsBase.sample(w)
        x = lowers[ii, 1] + step(xgrid) * rand()
        vr = lowers[ii, 2] + step(vrgrid) * rand()
        vt = lowers[ii, 3] + step(vtgrid) * rand()
        u = rand()
        sampleDensity = x * x * vt * profile_density(x, vr, vt, p, αt = αt)
        if rate
            tested += 1
        end
        if sampleDensity >= (u * fmax)
            accepted += 1
            sampledr[accepted] = x * p.rs
            sampledvr[accepted] =  vr
            sampledvt[accepted] = vt
        end
    end
    if rate
        return SphericalGalaxy(sampledr, sampledvr, sampledvt, samplesize), accepted / tested
    end
    return SphericalGalaxy(sampledr, sampledvr, sampledvt, samplesize)
end


text:
type DumbFeatureExtractor <: AbstractFeatureExtractor
end

let

    roadway = gen_straight_roadway(3, 1000.0, lane_width=1.0)
    rec = SceneRecord(1, 0.1, 5)
    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2( 0.0,0.0,0.0), roadway, 10.0), VehicleDef(1, AgentClass.CAR, 5.0, 2.0)),
            Vehicle(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(2, AgentClass.CAR, 5.0, 2.0)),
        ]))

    ext = DumbFeatureExtractor()
    @test rec_length(ext) == 1
    @test_throws Exception length(ext)
    @test_throws Exception pull_features!(ext, Array(Float64, 0), rec, roadway, 1)

    ext = FeatureExtractor(AbstractFeature[SPEED, POSFT], 2)
    @test rec_length(ext) == 2
    @test length(ext) == 2
    @test pull_features!(ext, [NaN,NaN], rec, roadway, 1) == [10.0, 0.0]

    subset_ext = SubsetExtractor(ext, [1])
    @test rec_length(subset_ext) == 2
    @test length(subset_ext) == 1
    @test pull_features!(subset_ext, [NaN], rec, roadway, 1) == [10.0]

    subset_ext.subset[1] = 2
    @test pull_features!(subset_ext, [NaN], rec, roadway, 1) == [0.0]

    stan_ext = StandardizingExtractor(ext, [1.0, 2.0], [1.0, 4.0])
    @test rec_length(stan_ext) == 2
    @test length(stan_ext) == 2
    @test pull_features!(stan_ext, [NaN,NaN], rec, roadway, 1) == [(10.0-1.0)/1.0, (0.0-2.0)/4.0]
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

path = "examples/examples_files/fourbar.urdf"
mech = Mechanism(path)

initializeConstraints!(mech, newtonIter=1000)

steps = Base.OneTo(1000)
storage = Storage{Float64}(steps,4)

storage = simulate!(mech, storage, record = true)
visualize(mech, storage)


text:
"""
    function agentArguments_(abm::Agent)

Returns the list of vectors arguments to be declared in the model function arguments.
"""
function agentArguments_(abm::Agent)
    l = [:t,:N]

    #Floats
    if length(abm.declaredSymbols["Local"])>0
        push!(l,:loc_)
    end
    if length(abm.declaredSymbols["Identity"])>0
        push!(l,:id_)
    end
    if length(abm.declaredSymbols["Global"])>0
        push!(l,:glob_)
    end
    if length(abm.declaredSymbols["GlobalArray"])>0
        for i in abm.declaredSymbols["GlobalArray"]
            push!(l,i)
        end
    end

    return l
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using JudiLing
using CSV
using Test
using SparseArrays
using DataFrames

@testset "cholesky" begin
    C = [1 1 1 1 0 0 0 0; 1 0 1 0 1 0 1 0; 0 0 0 0 1 1 1 1]
    S = [1 0 1 0; 1 1 0 0; 0 0 1 1]

    F = JudiLing.make_transform_matrix(C, S)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    JudiLing.make_transform_matrix(C, S, method = :multiplicative)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    C = sparse(C)
    S = Matrix(S)

    JudiLing.make_transform_matrix(C, S)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    JudiLing.make_transform_matrix(S, C)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    JudiLing.make_transform_matrix(C, S, method = :multiplicative)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    JudiLing.make_transform_matrix(S, C, method = :multiplicative)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    C = sparse(C)
    S = sparse(S)

    JudiLing.make_transform_matrix(C, S)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    JudiLing.make_transform_matrix(C, S, method = :multiplicative)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    fac_C = JudiLing.make_transform_fac(C)

    JudiLing.make_transform_matrix(fac_C, C, S)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    C = Matrix(C)
    fac_C = JudiLing.make_transform_fac(C)

    JudiLing.make_transform_matrix(fac_C, C, S)
    Shat = C * F
    @test -0.05 < sum(Shat-S) < 0.05

    JudiLing.make_transform_matrix(fac_C, C, S,
        output_format = :auto, sparse_ratio = 0.9)

    JudiLing.make_transform_matrix(fac_C, C, S,
        output_format = :sparse)
end


text:
include("low_volumefraction.jl")

using JLD
using Plots
plotly()

strfs = ["f_slow1","f_slow2","f_fast1","f_fast2"];
f = strfs[4]
species = load("data/$f.jld")["species"]
ωs = load("data/$f.jld")["ωs"]
medium = load("data/$f.jld")["medium"]

kTs_arr = [
    [ sqrt(wavenumber_low_volumefraction(ω,medium, sps)) for ω in ωs]
for sps in [[species[1]],[species[2]],species]];

speed_arr = [ ωs./real(kTs) for kTs in kTs_arr]
atten_arr = imag(kTs_arr)

labs = ["void" "stone" "mix"];
p1 = plot(ωs./real(medium.c),speed_arr, xlabel="k", ylabel="wave speed", labels=labs,);
p2 = plot(ωs./real(medium.c),atten_arr, xlabel="k", ylabel="attenuation", labels=labs,);
plot(p1,p2)

max_stars_repo_name: arturgower/EffectiveWaves


text:
@testset "sql" begin

using DataFrames

sess = SparkSession()

Spark.set_log_level(Spark.context(sess), "ERROR")

# testing IO
ds = read_json(sess, joinpath(@__DIR__, "people.json"))
@test count(ds) == 2

mktempdir() do dir
    parquet_file = joinpath(dir, "people.parquet")

    write_parquet(ds, parquet_file)
    ds2 = read_parquet(sess, parquet_file)
    write_json(ds2, joinpath(dir, "people2.json"))

    @test count(ds2) == 2
end

# testing sql
Spark.create_temp_view(ds, "people")
ds_all = sql(sess, "SELECT * from people")

@test Spark.as_named_tuple(Spark.head(ds_all)) == (age = 32, name = "Peter")

# test conversion to DataFrames
df = DataFrame(ds_all)
@test nrow(df) == 2
@test ncol(df) == 2

close(sess)

end


text:
module PanelMatrices

using StaticNumbers
using StaticArrays
using LinearAlgebra
#using UnsafeArrays # maybe we don't need this, but just use the same idea.

export PanelMatrix

const panel_stride_error = ErrorException("Panel stride is less than minimum required.")
const too_small_error = ErrorException("PanelMatrix size must either be static or greater than the panel size.")

# The inner PanelMatrix constructor checks the size of the data vector.
# Use @inbounds to bypass size and stride check.
# (Mutating the size of .data may lead to memory corroption. Don't do that!)
# We force the size to be static when smaller than or equal to one panel.
# All subsequent code can assume that if size[i] isa Int then that Int
# is not smaller than panel_size[i] and therefore that pat_first[i] and
# pad_last[i] never apply to the same panel.
# pad_last is computed from pad_first, size and panel_size. Only its type can be specified.
# n_panels is computed from pad_first, size and panel_size. It will be Static if all of those are.

"""
A `PanelMatrix` is a panel-major matrix, or a view into a panel-major matrix.

The size of a PanelMatrix can either be Int or StaticInteger. In either case,
the  matrix cannot be resized. (Using Int requires fewer methods to be compiled
when matrices of varying sizes are used.)

If either the width or height of a PanelMatrix is smaller than the default panel
size, then it is recommended to use a StaticInteger for that element of the size,
and to reduce the panel size so that it does not exceed the matrix size.

A PanelMatrix contains a `.data` field which holds the actual data. This must
be an array that maps to continuous memory storage. (For efficiency, panels may
be accessed using hand-coded assembly that makes this assumtion.)

Care should be taken if mutating the `.data` array. Changing its size
(e.g. using `push!` or `pop!`) may lead to data corruption.

Memory is allocated to cover the entire panel, even in cases where only parts
of a panel is covered by the matrix.
"""
struct PanelMatrix{T,D,S1,S2,P1,P2,C,F1,F2,L1,L2,R,N1,N2} <: AbstractMatrix{T}
    data::D # The actual data <: AbstractVector{T}
    size::Tuple{S1, S2} # (rows, columns) in the matrix
    panel_size::Tuple{P1, P2} # (rows, columns) in each panel
    panel_class::Val{C}  # :CM = column major, :RM = row major
    pad_first::Tuple{F1, F2} # (rows, columns) of padding before the first element of first tile
    pad_last::Tuple{L1, L2} # (rows, columns) of padding after the last element until last tile ends
    panel_stride::R # number of panels per column in data
    n_panels::Tuple{N1, N2} # (rows, columns) of panels in the matrix

    function PanelMatrix{T,D}(data::D, size::Tuple{S1, S2}, panel_size::Tuple{P1, P2},
           panel_class::Val{C}, pad_first::Tuple{F1, F2}, pad_last_types::Tuple{<:Union{Nothing,Type}, <:Union{Nothing,Type}},
           panel_stride::R) where {T, D<:AbstractVector{T}, S1<:Integer, S2<:Integer, P1<:Integer, P2<:Integer, C, F1<:Integer, F2<:Integer, R<:Integer}

        n_panels = maybe_static.((f,s,p) -> cld.(f+s,p), pad_first, size, panel_size)
        @boundscheck begin
            panel_stride * panel_size[1] >= size[1] + pad_first[1] || throw(panel_stride_error)
            last_panel = n_panels[1] + (n_panels[2]-1)*panel_stride
            checkbounds(data, Base.OneTo(last_panel*prod(panel_size)))
        end
        any(.!(isa.(StaticInteger, typeof.(size))) .& (size .<= panel_size)) && throw(too_small_error)
        pad_last_types = something.(pad_last_types, typeof.(size))
        pad_last = convert.(pad_last_types, panel_size .* n_panels .- pad_first .- size)
        (L1, L2) = typeof.(pad_last)
        (N1, N2) = typeof.(n_panels)
        new{T,D,S1,S2,P1,P2,C,F1,F2,L1,L2,R,N1,N2}(data, size, panel_size, panel_class, pad_first, pad_last, panel_stride, n_panels)
    end
end

# TODO: We could use keyword args, but first benchmark without and check that
#       there's no performance issue.

"""
    PanelMatrix{T}(undef, size, panel_size, panel_class, pad_first, pad_last_types)

Create an uninitialized `PanelMatrix` for element type `T`, with the given `size`.
"""
function PanelMatrix{T}(
       ::UndefInitializer,
       size::Tuple{<:Integer, <:Integer},
       panel_size::Tuple{<:Integer, <:Integer} = static.(min.(size, default_panelsize(T))),
       panel_class::Val = default_panelclass(T),
       pad_first::Tuple{<:Integer, <:Integer} = static.((0, 0)),
       pad_last_types::Tuple = (nothing, nothing)
       ) where {T}
    panel_stride = cld(size[1]+pad_first[1], panel_size[1])
    n_panels = cld.(pad_first .+ size, panel_size)
    last_panel = n_panels[1] + (n_panels[2]-1)*panel_stride
    data = Vector{T}(undef, last_panel*prod(panel_size))
    PanelMatrix{T, Vector{T}}(data, size, panel_size, panel_class, pad_first, pad_last_types, panel_stride)
end

"""
    PanelMatrix(x, panel_size, panel_class, pad_first, pad_last_types)

Create a `PanelMatrix` from the matrix `x` (making a copy).
"""
function PanelMatrix(x::AbstractMatrix,
        panel_size::Tuple{<:Integer, <:Integer} = static.(min.(size(x), default_panelsize(eltype(x)))),
        panel_class::Val = default_panelclass(T),
        pad_first::Tuple{<:Integer, <:Integer} = static.((0, 0)),
        pad_last_types::Tuple = (nothing, nothing)
        )
    z = PanelMatrix{eltype(x)}(undef, size(x), panel_size, panel_class, pad_first, pad_last_types)
    z .= x
    #  TODO: fill padding with zeros ?
    return z
end

default_panelsize(T) = static.((4, 4))

default_panelclass(T) = Val(:CM)

include("Panels.jl")

@inline function Base.getindex(x::PanelMatrix, i, j)
    @boundscheck checkbounds(x, i, j)
    (p, q) = divrem.((i, j) .+ x.pad_first .- 1, x.panel_size)
    @inbounds full_panel_view(x, p[1]+1, q[1]+1)[p[2]+1, q[2]+1]
end

@inline function Base.setindex!(x::PanelMatrix, v, i, j)
    @boundscheck checkbounds(x, i, j)
    (p, q) = divrem.((i, j) .+ x.pad_first .- 1, x.panel_size)
    @inbounds full_panel_view(x, p[1]+1, q[1]+1)[p[2]+1, q[2]+1] = v
end

Base.size(x::PanelMatrix) = x.size
Base.eltype(x::PanelMatrix) = eltype(x.data)

# TODO: Base.similar

include("matmul.jl")

end # module


text:
using CSV
using DataFrames
using LinearAlgebra
using Plots

function ansatz(x,α)
	x = norm(x)
	return exp(-α*x)
end

function probability_density(x,α)
	# note that it is not normalised becuase it does not matter!
	return ansatz(x,α)^2
end

function E_local(x,α)
	x = norm(x)
	return -1/x - (α/2)*(α - 2/x)
end

function metropolis(N,α)
	L = 3/sqrt(2*α)
	x = rand(3)*2*L .- L
	E = 0
	E2 = 0
	Eln_avg = 0
	ln_avg = 0
	rejection_ratio = 0
	#start the loop
	for i in 1:N
		x_t = x + 0.2*(rand(3)*2*L .- L)
		if probability_density(x_t,α) >= probability_density(x,α)
			x = x_t
		else
			rand_number = rand()
			if rand_number < probability_density(x_t,α)/probability_density(x,α)
				x = x_t
			else
				rejection_ratio += 1/N
			end
		end

		E += E_local(x,α)/N
		E2 += E_local(x,α)^2/N
		Eln_avg += -1 * E_local(x,α) * (norm(x)/N)
        ln_avg += -1 * norm(x)/N

	end
	return E,E2,Eln_avg,ln_avg,rejection_ratio
end

function iterate_alpha()
	α = 0.2
	alpha_iter = 30
	N = 500
	random_walkers = 200
	γ = 0.5
	energy = Vector{Float64}(undef,alpha_iter)
	alpha = Vector{Float64}(undef,alpha_iter)
	variance = Vector{Float64}(undef,alpha_iter)
	for i in 1:alpha_iter
		E = 0
		E2 = 0
		dE_dalpha = 0
		Eln = 0
		ln = 0
		rejection_ratio = 0
		for j in 1:random_walkers
			E, E2, Eln, ln, rejections = metropolis(N,α)
			E += E/random_walkers
			E2 += E2/random_walkers
			Eln += Eln/random_walkers
			ln += ln/random_walkers
			# rejection_ratio += rejections_met/random_walkers
		end

		# Define the next α
		dE_dalpha = 2*(Eln - E*ln)
		α = α - γ*dE_dalpha
		# α = α + 0.05

		# Update the arrays for plotting
		energy[i] = E
		alpha[i] = α
		variance[i] = E2 - E^2

	end
	return energy,alpha,variance
end

energy,alpha,variance = iterate_alpha()
CSV.write(datadir("sims/hydrogen_atom.csv"),(Energy = energy, Variance = variance, Alpha = alpha))

energy_plot = plot()
plot!(energy_plot,energy,
	label = false,
	title = "Energy expenctation Plot",
	linecolor = :blue,
	ylabel = "<E>",
	xlabel = "timestep"
 )

safesave(plotsdir("energy_hydrogen_atom.png"),energy_plot)

variance_plot = plot()
plot!(variance_plot,variance,
	title = "Variance in energy",
	label = false,
	linecolor = :blue,
	ylabel = "Var(E)",
	xlabel = "timestep"
     )
safesave(plotsdir("variance_hydrogen_atom.png"),variance_plot)

alpha_plot = plot()
plot!(alpha_plot,alpha,
	title = "Variation Parameter",
	label = false,
	linecolor = :blue,
	ylabel = "α",
	xlabel = "timestep"
)

safesave(plotsdir("alpha_hydrogen_atom.png"),alpha_plot)

alpha_vs_energy_plot = plot()
scatter!(alpha_vs_energy_plot,alpha,
	energy,
	title = "Energy vs α",
	xlabel = "α",
	ylabel = "<E>",
	label = false,
	markercolor = :blue
)
safesave(plotsdir("alpha_vs_energy_hydrogen_atom.png"),alpha_vs_energy_plot)

max_stars_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_issues_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_forks_repo_name: Ved-Mahajan/Variational-Monto-Carlo


text:
# General structure for reachability methods

# Performs layer-by-layer propagation
# It is called by all solvers under reachability
# TODO: also called by ReluVal and FastLin, so move to general utils (or to network.jl)
function forward_network(solver, nnet::Network, input::AbstractPolytope)
    reach = input
    for layer in nnet.layers
        reach = forward_layer(solver, layer, reach)
    end
    return reach
end

# Checks whether the reachable set belongs to the output constraint
# It is called by all solvers under reachability
# Note vertices_list is not defined for HPolytope: to be defined
function check_inclusion(reach::Vector{<:AbstractPolytope}, output)
    for poly in reach
        issubset(poly, output) || return ReachabilityResult(:violated, reach)
    end
    return ReachabilityResult(:holds, similar(reach, 0))
end

function check_inclusion(reach::P, output) where P<:AbstractPolytope
    if issubset(reach, output)
        return ReachabilityResult(:holds, P[])
    end
    return ReachabilityResult(:violated, [reach])
end

# return a vector so that append! is consistent with the relu forward_partition
forward_partition(act::Id, input::HPolytope) = [input]

function forward_partition(act::ReLU, input::HPolytope)
    n = dim(input)
    output = Vector{HPolytope}(undef, 0)
    C, d = tosimplehrep(input)
    dh = [d; zeros(n)]
    for h in 0:(2^n)-1
        P = getP(h, n)
        Ch = [C; I - 2P]
        input_h = HPolytope(Ch, dh)
        if !isempty(input_h)
            push!(output, linear_map(Matrix{Float64}(P), input_h))
        end
    end
    return output
end

function getP(h::Int64, n::Int64)
    str = string(h, pad = n, base = 2)
    vec = Vector{Int64}(undef, n)
    for i in 1:n
        vec[i] = ifelse(str[i] == '1', 1, 0)
    end
    return Diagonal(vec)
end


text:
module Net

export sctenifoldnet, tenrnet, manialn, drgenes, pcnet, tensordecomp, tensordecomp2

const NCOMP1,NCOMP2=3,5
const NLAYERS,NCELLS=10,500

using Statistics, LinearAlgebra, Distributions, Random
using MultipleTesting, Random, SparseArrays, TSVD, TensorToolbox
# import KrylovKit

include("support.jl")
include("pcnet.jl")                   # include the contents of other files in the module
include("tenrnet.jl")
include("tensordecomp.jl")
include("manialn.jl")
include("drgenes.jl")

function sctenifoldnet(X::AbstractMatrix{T}, Y::AbstractMatrix{T}; donorm::Bool=true) where T<:Real
    Z0,_=tenrnet(X,donorm=donorm)
    Z1,_=tenrnet(Y,donorm=donorm)
    Z0=0.5*(Z0+Z0')
    Z1=0.5*(Z1+Z1')
    d=manialn(Z0,Z1)
    fc,p,adjp=drgenes(d)
    return d,fc,p,adjp
end

end

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
module DiskData

using Lazy

typealias AVector AbstractVector

include("iter/iter.jl")
include("iter/split.jl")
include("vectors/chunked.jl")
include("vectors/cached.jl")
include("vectors/disk.jl")
include("vectors/bigvector.jl")

end # module


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

"""
The [Anholt Wind Farm](https://en.wikipedia.org/wiki/Anholt_Offshore_Wind_Farm)

"""

# Load pkgs
using Plots
using Unitful
using MultiStateSystems

# Pkg const
const _MSS = MultiStateSystems

# Settings for a specific analysis
cbl = true
mss = true
wtr = true
number_of_clusters = 8

# Include the stochastic models for the wind turbines and cables
include(joinpath(_MSS.BASE_DIR,"examples/anholt_wind_farm/elements/wind_turbine.jl"))
include(joinpath(_MSS.BASE_DIR,"examples/anholt_wind_farm/elements/sea_cable.jl"))

# Feeder 1 - Nodes 30, 62:65, 76, 109:111
ntwᶠ¹ = Network()
add_user!(ntwᶠ¹, node = 1)
add_sources!(ntwᶠ¹, node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ¹, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(5,7),(3,8),(8,9),(9,10)],
                             std = [stdᶜᵇˡ[(8.500u"km",500u"mm^2")],stdᶜᵇˡ[(1.678u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.809u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(1.464u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")]]) :
     add_components!(ntwᶠ¹, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(5,7),(3,8),(8,9),(9,10)]) ;

# Feeder 2 - Nodes 23:29, 42, 61
ntwᶠ² = Network()
add_user!(ntwᶠ², node = 1)
add_sources!(ntwᶠ², node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ², edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)],
                             std = [stdᶜᵇˡ[(4.500u"km",500u"mm^2")],stdᶜᵇˡ[(0.603u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",240u"mm^2")],stdᶜᵇˡ[(0.603u"km",240u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",240u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.809u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.809u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ², edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)]) ;

# Feeder 3 - Nodes 22, 59, 60 102:108,
ntwᶠ³ = Network()
add_user!(ntwᶠ³, node = 1)
add_sources!(ntwᶠ³, node = 2:11, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ³, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(5,10),(10,11)],
                             std = [stdᶜᵇˡ[(4.000u"km",500u"mm^2")],stdᶜᵇˡ[(2.100u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(0.691u"km",500u"mm^2")],stdᶜᵇˡ[(1.678u"km",240u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ³, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(5,10),(10,11)]) ;

# Feeder 4 - Nodes 19:21, 56:58, 99:101
ntwᶠ⁴ = Network()
add_user!(ntwᶠ⁴, node = 1)
add_sources!(ntwᶠ⁴, node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ⁴, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(3,7),(7,8),(2,9),(9,10)],
                             std = [stdᶜᵇˡ[(2.000u"km",500u"mm^2")],stdᶜᵇˡ[(2.100u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(2.100u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.873u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.873u"km",150u"mm^2")],stdᶜᵇˡ[(0.691u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(1.272u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ⁴, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(3,7),(7,8),(2,9),(9,10)]) ;

# Feeder 5 - Nodes 17,18, 41, 54,55, 75, 96:98
ntwᶠ⁵ = Network()
add_user!(ntwᶠ⁵, node = 1)
add_sources!(ntwᶠ⁵, node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ⁵, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(6,8),(4,9),(2,10)],
                             std = [stdᶜᵇˡ[(2.000u"km",500u"mm^2")],stdᶜᵇˡ[(1.464u"km",240u"mm^2")],
                                    stdᶜᵇˡ[(1.095u"km",240u"mm^2")],stdᶜᵇˡ[(1.095u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(1.095u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.873u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.691u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ⁵, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(6,8),(4,9),(2,10)]) ;

# Feeder 6 - Nodes 15,16, 40, 52,53, 74, 92:95
ntwᶠ⁶ = Network()
add_user!(ntwᶠ⁶, node = 1)
add_sources!(ntwᶠ⁶, node = 2:11, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ⁶, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(6,9),(4,10),(2,11)],
                             std = [stdᶜᵇˡ[(2.000u"km",500u"mm^2")],stdᶜᵇˡ[(1.272u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(1.095u"km",240u"mm^2")],stdᶜᵇˡ[(1.272u"km",240u"mm^2")],
                                    stdᶜᵇˡ[(1.464u"km",150u"mm^2")],stdᶜᵇˡ[(0.691u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.809u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.873u"km",150u"mm^2")],stdᶜᵇˡ[(0.691u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ⁶, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(6,9),(4,10),(2,11)]) ;

# Feeder 7 - Nodes 14, 39, 51, 73, 83,84, 88:91
ntwᶠ⁷ = Network()
add_user!(ntwᶠ⁷, node = 1)
add_sources!(ntwᶠ⁷, node = 2:11, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ⁷, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10),(6,11)],
                             std = [stdᶜᵇˡ[(3.000u"km",500u"mm^2")],stdᶜᵇˡ[(1.272u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(1.272u"km",500u"mm^2")],stdᶜᵇˡ[(1.464u"km",240u"mm^2")],
                                    stdᶜᵇˡ[(1.095u"km",240u"mm^2")],stdᶜᵇˡ[(1.272u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.691u"km",150u"mm^2")],stdᶜᵇˡ[(0.691u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.691u"km",150u"mm^2")],stdᶜᵇˡ[(0.809u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ⁷, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10),(6,11)]) ;

# Feeder 8 - Nodes 13, 50, 72, 80:82, 85:87
ntwᶠ⁸ = Network()
add_user!(ntwᶠ⁸, node = 1)
add_sources!(ntwᶠ⁸, node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ⁸, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)],
                             std = [stdᶜᵇˡ[(3.000u"km",500u"mm^2")],stdᶜᵇˡ[(2.100u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(1.464u"km",240u"mm^2")],stdᶜᵇˡ[(1.272u"km",240u"mm^2")],
                                    stdᶜᵇˡ[(0.809u"km",240u"mm^2")],stdᶜᵇˡ[(0.809u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(1.464u"km",150u"mm^2")],stdᶜᵇˡ[(0.809u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.691u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ⁸, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)]) ;

# Feeder 9 - Nodes 12, 38, 49, 69:71, 77:79
ntwᶠ⁹ = Network()
add_user!(ntwᶠ⁹, node = 1)
add_sources!(ntwᶠ⁹, node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ⁹, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)],
                             std = [stdᶜᵇˡ[(4.000u"km",500u"mm^2")],stdᶜᵇˡ[(1.464u"km",500u"mm^2")],
                                    stdᶜᵇˡ[(1.272u"km",240u"mm^2")],stdᶜᵇˡ[(1.678u"km",240u"mm^2")],
                                    stdᶜᵇˡ[(0.873u"km",240u"mm^2")],stdᶜᵇˡ[(0.873u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.952u"km",150u"mm^2")],stdᶜᵇˡ[(0.809u"km",150u"mm^2")],
                                    stdᶜᵇˡ[(0.809u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ⁹, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)]) ;

# Feeder 10 - Nodes 11, 37, 45:48, 66:68
ntwᶠ¹⁰ = Network()
add_user!(ntwᶠ¹⁰, node = 1)
add_sources!(ntwᶠ¹⁰, node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ¹⁰, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(6,8),(8,9),(9,10)],
                              std = [stdᶜᵇˡ[(4.500u"km",500u"mm^2")],stdᶜᵇˡ[(1.678u"km",500u"mm^2")],
                                     stdᶜᵇˡ[(1.464u"km",240u"mm^2")],stdᶜᵇˡ[(0.952u"km",240u"mm^2")],
                                     stdᶜᵇˡ[(0.952u"km",240u"mm^2")],stdᶜᵇˡ[(0.952u"km",150u"mm^2")],
                                     stdᶜᵇˡ[(0.873u"km",150u"mm^2")],stdᶜᵇˡ[(0.809u"km",150u"mm^2")],
                                     stdᶜᵇˡ[(0.809u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ¹⁰, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(6,8),(8,9),(9,10)]) ;

# Feeder 11 - Nodes 10, 31:36, 43,44
ntwᶠ¹¹ = Network()
add_user!(ntwᶠ¹¹, node = 1)
add_sources!(ntwᶠ¹¹, node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ¹¹, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(6,9),(9,10)],
                              std = [stdᶜᵇˡ[(5.000u"km",500u"mm^2")],stdᶜᵇˡ[(1.678u"km",500u"mm^2")],
                                     stdᶜᵇˡ[(0.873u"km",240u"mm^2")],stdᶜᵇˡ[(0.873u"km",240u"mm^2")],
                                     stdᶜᵇˡ[(0.873u"km",240u"mm^2")],stdᶜᵇˡ[(0.873u"km",150u"mm^2")],
                                     stdᶜᵇˡ[(0.809u"km",150u"mm^2")],stdᶜᵇˡ[(0.691u"km",150u"mm^2")],
                                     stdᶜᵇˡ[(0.809u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ¹¹, edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(6,9),(9,10)]) ;

# Feeder 12 - Nodes
ntwᶠ¹² = Network()
add_user!(ntwᶠ¹², node = 1)
add_sources!(ntwᶠ¹², node = 2:10, ntw = (ntwʷᵗ,1))
cbl ? add_components!(ntwᶠ¹², edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)],
                              std = [stdᶜᵇˡ[(5.750u"km",500u"mm^2")],stdᶜᵇˡ[(0.603u"km",500u"mm^2")],
                                     stdᶜᵇˡ[(0.603u"km",240u"mm^2")],stdᶜᵇˡ[(0.603u"km",240u"mm^2")],
                                     stdᶜᵇˡ[(0.603u"km",240u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                     stdᶜᵇˡ[(0.603u"km",150u"mm^2")],stdᶜᵇˡ[(0.603u"km",150u"mm^2")],
                                     stdᶜᵇˡ[(0.603u"km",150u"mm^2")]]) :
      add_components!(ntwᶠ¹², edge = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,10)]) ;

# Feeders 1-6
ntwᶠ¹⁻⁶ = Network()
add_user!(ntwᶠ¹⁻⁶, node = 1)
add_sources!(ntwᶠ¹⁻⁶, node = 1, ntw = [(ntwᶠ¹,1),(ntwᶠ²,1),(ntwᶠ³,1),(ntwᶠ⁴,1),(ntwᶠ⁵,1),(ntwᶠ⁶,1)])

# Feeders 7-12
ntwᶠ⁷⁻¹² = Network()
add_user!(ntwᶠ⁷⁻¹², node = 1)
add_sources!(ntwᶠ⁷⁻¹², node = 1, ntw = [(ntwᶠ⁷,1),(ntwᶠ⁸,1),(ntwᶠ⁹,1),(ntwᶠ¹⁰,1),(ntwᶠ¹¹,1),(ntwᶠ¹²,1)])

# Overall network
ntw = Network()
add_user!(ntw, node = 1, ind = [:GRA,:EENS])
add_sources!(ntw, node = 1, ntw = [(ntwᶠ¹⁻⁶,1),(ntwᶠ⁷⁻¹²,1)])

# Solve the problem
@time solve!(ntw)

# Plot the probability distribution (p>=1e-5)
ugf = ntw.usr[1][:ugf]
scatter(ustrip.(ugf.val[ugf.prb.>=1e-5]),ugf.prb[ugf.prb.>=1e-5],legend=false)


text:
using GeneratorArrays
using Test

struct Foo0 end
struct Foo1
    x
end
struct Foo2
    x
    y
end

@testset "GeneratorArrays" begin

    @testset "iterators" begin
        let i0 = 1:0.5:5, i1 = [1;2;3]; i2 = [Float16(1) 2; 3 4; 5 6]
        @testset "arrays $it" for (it,T,N) in ((i0,Float64,1),
                                               (i1,Int,1),
                                               (i2,Float16,2))
            @test array(it) === it
            gen = (x for x in it)
            @test length(axes(array(gen))) == N
            agen = array(gen)
            @test axes(agen) == axes(it)
            @test size(agen) == size(it)
            @test length(agen) == length(it)
            @test iterate(agen) == iterate(it)
            @test agen == collect(it)
            @test agen[3] == it[3]

            @test Base.IteratorSize(typeof(agen)) == Base.HasShape{N}()
            @test Base.IteratorEltype(typeof(agen)) == Base.HasEltype()
        end

        @test array("abγ") == ['a', 'b', 'γ']
        end
    end

    @testset "generators" begin
        a = [1.0; 2; 3]
        let g1 = (Float16(sin(x)) for x in 1:5), g3 = (x*y for x in a, y in 1:5, z = 1:2)
        @testset "generator dim $N" for (it,T,N,s) in ((g1,Float16,1,(5,)),
                                                    (g3,Float64,3,(3,5,2)))
            @test array(it).gen === it
            agen = array(it)
            @test length(axes(agen)) == N
            @test size(agen) == s
            @test length(agen) == prod(s)
            @test collect(agen) == collect(it)
            @test agen == collect(it)
            @test agen[1] == iterate(it)[1]
        end
        agen = array(g1)
        @test agen[3] == Float16(sin(3.0))
        agen = array(g3)
        @test agen[1] == 1
        @test agen[1,2] == 2
        @test agen[1,2,2] == 2
        @test_throws BoundsError agen[1,2,3]
        @test agen[1,2,1,1] == 2
        @test_throws BoundsError agen[]
        @test_throws BoundsError agen[1,1,1,17]
        @test agen[1:2] == a[1:2]
        @test agen[:,2,1] == a .* 2
        @test agen[:,:,1] == a * (1:5)'
        @test_throws DomainError array((x*y for x in a for y in a))
        end
    end

    @testset "recursive" begin
        let g0 = (x^2 for x in -2:2), g1 = (x*y for x in 1:2, y in 0:10)

            g2 = ((s, t) for s in g0, t in g1)
            a2 = array(g2)
            @test size(a2) == (5,2,11)
            @test eltype(a2) == Tuple{Int,Int}
            @test a2 == collect(g2)
        end
    end

    let constructors = ((Foo0, 0), (Foo1, 1), (Foo2, 2)),
        iterators = ( (), (1,), (1,2), Int[], [1], [1,2])

        @testset "constructors with $F - $b" for b in iterators, (F, narg) in constructors
            FF(i) = i == narg ? F : Base.Bottom
            narg == 0 && @test eltype(array(F() for x in b)) == FF(0)
            narg == 1 && @test eltype(array(F(x) for x in b)) == FF(1)
            narg == 2 && @test eltype(array(F(x,x) for x in b)) == FF(2)
        end
    end

    let iterators = ( ((), Union{}), ((1,), Tuple{Int,Int}))
        @testset "product iterators $b" for (b, res)  in iterators
            g = ((x,y) for x in b, y in 1:2)
            @test eltype(array(g)) == res
        end
    end
end


text:
function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.NUM)

    p[C.S] = 2.0
    p[C.kE] = 0.4
    p[C.kM] = 1.0
    p[C.kCD] = 0.03
    p[C.kCDS] = 0.45
    p[C.kR] = 0.18
    p[C.kRE] = 180
    p[C.kb] = 0.003
    p[C.KS] = 0.5
    p[C.kCE] = 0.35
    p[C.dM] = 0.7
    p[C.dE] = 0.25
    p[C.dCD] = 1.5
    p[C.dCE] = 1.5
    p[C.dR] = 0.06
    p[C.dRP] = 0.06
    p[C.dRE] = 0.03
    p[C.kP1] = 18.0
    p[C.kP2] = 18.0
    p[C.kDP] = 3.6
    p[C.KM] = 0.15
    p[C.KE] = 0.15
    p[C.KCD] = 0.92
    p[C.KCE] = 0.92
    p[C.KRP] = 0.01

    return p
end

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
"""
The NoiseModel supertype corresponds to different methods
    of simulating noise coming from the electronics chain
    components.

It is not simply a string in the ElecChain struct because
    it must exist as a separate instance for the stp->pss step
    that is independent from the elec simulation, but has to know
    what type of noise model is used to deal with fano noise.
"""
abstract type NoiseModel end

"""
The NoiseFromSim model means simulating noise from scratch
    starting from fano noise of the germanium crystal
    and ending with noise coming from electronics components.
"""
@with_kw struct NoiseFromSim <: NoiseModel
    "Noise σ in keV for Gaussian noise simulation"
    noise_σ::typeof(1.0*energy_unit) = 0
end

"""
    NoiseFromSim(sim_conf)

PropDict -> NoiseFromSim

Construct a NoiseFromSim struct based on given simulation configuration <sim_conf>
"""
function NoiseFromSim(sim_conf::PropDict)
    @info "//\\//\\// Noise simulated from scratch (fano, preamp noise)"
    noise_σ = haskey(sim_conf.setup, :preamp) ? T(sim_conf.setup.preamp.noise_sigma)u"keV" : 0u"keV"

    NoiseFromSim(noise_σ)
end

"""
The NoiseFromData model means instead of simulating given
    preamp noise, and also gain and offset, these parameters will be
    accounted for or inferred from real data baselines
    contained in the <baseline_catalog>
"""
struct NoiseFromData <: NoiseModel
    "A table of data baselines in RDWaveform format"
    baseline_catalog::Table
end

"""
    NoiseFromData(sim_conf)

PropDict -> NoiseFromData

Construct a NoiseFromData struct based on given simulation configuration <sim_conf>
"""
function NoiseFromData(sim_conf::PropDict)
    @info "//\\//\\// Noise levels and offset added via slapping the baseline from data on top of the waveform"
    # construct table of baselines based on given raw data hdf5 file
    baseline_table = baseline_catalog(sim_conf.noise_data)
    NoiseFromData(baseline_table)
end

"""
    NoiseModel(sim_config)

PropDict -> <NoiseModel>

Constuct a NoiseModel supertype instance based on simulation settings
    given in <sim_config>
Type of <NoiseModel> depends on <sim_config> settings.
"""
function NoiseModel(sim_config::PropDict)
    if haskey(sim_config, :noise_data)
        NoiseFromData(sim_config)
    else
        NoiseFromSim(sim_config)
    end
end

# -------------------------------------------------------------------

"""
    fano_noise(events, det_meta, env, ::NoiseFromSim)

Table, PropDict, Environment -> Table

Calculate fano noise level based on the detector specification provided in
    LEGEND metadata <det_meta> and environment settings provided in <env>,
    and add it to given <events>
"""
function fano_noise(events::Table, det_meta::PropDict, env::Environment, ::NoiseFromSim)
    println("Adding fano noise")
    ssd_conf = ssd_config(det_meta, env)
    simulation = Simulation(SolidStateDetector{T}(ssd_conf))
    det_material = simulation.detector.semiconductors[1].material
    add_fano_noise(events, det_material.E_ionisation, det_material.f_fano)
end

"""
    fano_noise(events, ::PropDict, ::Environment, ::NoiseFromData)

Table -> Table

Do nothing since we do not need to simulate fano noise separately when
    using data baselines to account for noise levels.
"""
function fano_noise(events::Table, ::PropDict, ::Environment, ::NoiseFromData)
    println("Not adding fano noise because using noise levels from data")
    events
end

"""
    simulate_noise(wf, preamp)

RDWaveform, PreAmp -> RDWaveform

Simulate effects of the preamplifier <preamp> on the given waveform <wf>.
"""
function simulate_noise(wf::RDWaveform, preamp::PreAmp)
    # wf values are in eV (without u"eV" units attached), noise sigma is in keV
    noise_σ = ustrip(uconvert(u"eV", preamp.noise_σ))
    gaussian_noise_dist = Normal(T(0), T(noise_σ))
    RDWaveform(wf.time, wf.value .+ rand!(gaussian_noise_dist, similar(wf.value)))
end

"""
    simulate_noise(wf, noise_model)

RDWaveform, NoiseFromData -> RDWaveform

Simulate noise and offset by picking a random baseline from the baseline catalog
    contained in <noise_model> and slapping it on top of the given waveform <wf>.
"""
function simulate_noise(wf::RDWaveform, noise_model::NoiseFromData)
    # NOTE: has to be done after FADC i.e. after ElecChain
    baseline = rand(Tables.getcolumn(Tables.columns(noise_model.baseline_catalog), :waveform))
    # extend to match the wf
    baseline_long = extend_baseline(baseline, wf)
    # slap on top of the waveform as promised
    RDWaveform(wf.time, wf.value .+ baseline_long.value)
end


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

mutable struct ConnectivityInformation <: SwaggerModel
    hops::Any # spec type: Union{ Nothing, Vector{ConnectivityHop} } # spec name: hops
    connectionStatus::Any # spec type: Union{ Nothing, String } # spec name: connectionStatus
    avgLatencyInMs::Any # spec type: Union{ Nothing, Int32 } # spec name: avgLatencyInMs
    minLatencyInMs::Any # spec type: Union{ Nothing, Int32 } # spec name: minLatencyInMs
    maxLatencyInMs::Any # spec type: Union{ Nothing, Int32 } # spec name: maxLatencyInMs
    probesSent::Any # spec type: Union{ Nothing, Int32 } # spec name: probesSent
    probesFailed::Any # spec type: Union{ Nothing, Int32 } # spec name: probesFailed

    function ConnectivityInformation(;hops=nothing, connectionStatus=nothing, avgLatencyInMs=nothing, minLatencyInMs=nothing, maxLatencyInMs=nothing, probesSent=nothing, probesFailed=nothing)
        o = new()
        validate_property(ConnectivityInformation, Symbol("hops"), hops)
        setfield!(o, Symbol("hops"), hops)
        validate_property(ConnectivityInformation, Symbol("connectionStatus"), connectionStatus)
        setfield!(o, Symbol("connectionStatus"), connectionStatus)
        validate_property(ConnectivityInformation, Symbol("avgLatencyInMs"), avgLatencyInMs)
        setfield!(o, Symbol("avgLatencyInMs"), avgLatencyInMs)
        validate_property(ConnectivityInformation, Symbol("minLatencyInMs"), minLatencyInMs)
        setfield!(o, Symbol("minLatencyInMs"), minLatencyInMs)
        validate_property(ConnectivityInformation, Symbol("maxLatencyInMs"), maxLatencyInMs)
        setfield!(o, Symbol("maxLatencyInMs"), maxLatencyInMs)
        validate_property(ConnectivityInformation, Symbol("probesSent"), probesSent)
        setfield!(o, Symbol("probesSent"), probesSent)
        validate_property(ConnectivityInformation, Symbol("probesFailed"), probesFailed)
        setfield!(o, Symbol("probesFailed"), probesFailed)
        o
    end
end # type ConnectivityInformation

const _property_map_ConnectivityInformation = Dict{Symbol,Symbol}(Symbol("hops")=>Symbol("hops"), Symbol("connectionStatus")=>Symbol("connectionStatus"), Symbol("avgLatencyInMs")=>Symbol("avgLatencyInMs"), Symbol("minLatencyInMs")=>Symbol("minLatencyInMs"), Symbol("maxLatencyInMs")=>Symbol("maxLatencyInMs"), Symbol("probesSent")=>Symbol("probesSent"), Symbol("probesFailed")=>Symbol("probesFailed"))
const _property_types_ConnectivityInformation = Dict{Symbol,String}(Symbol("hops")=>"Vector{ConnectivityHop}", Symbol("connectionStatus")=>"String", Symbol("avgLatencyInMs")=>"Int32", Symbol("minLatencyInMs")=>"Int32", Symbol("maxLatencyInMs")=>"Int32", Symbol("probesSent")=>"Int32", Symbol("probesFailed")=>"Int32")
Base.propertynames(::Type{ ConnectivityInformation }) = collect(keys(_property_map_ConnectivityInformation))
Swagger.property_type(::Type{ ConnectivityInformation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectivityInformation[name]))}
Swagger.field_name(::Type{ ConnectivityInformation }, property_name::Symbol) =  _property_map_ConnectivityInformation[property_name]

const _allowed_ConnectivityInformation_connectionStatus = ["Unknown", "Connected", "Disconnected", "Degraded"]

function check_required(o::ConnectivityInformation)
    true
end

function validate_property(::Type{ ConnectivityInformation }, name::Symbol, val)
    if name === Symbol("connectionStatus")
        Swagger.validate_param(name, "ConnectivityInformation", :enum, val, _allowed_ConnectivityInformation_connectionStatus)
    end
end


text:
using GLVisualize, GLAbstraction, GLWindow, GeometryTypes, Reactive, ColorTypes, ModernGL

w, r = Screen();
@async r();
cs = cubecamera(w);

cos_01{T}(i::T) = (cos(i)+T(1.5)) / T(2.5)
sin_01{T}(i::T) = (sin(i)+T(1.5)) / T(2.5)
glClearColor(1,1,1,1)

cap(a, lower=0f0, upper=1f0) = min(max(a, lower), upper)

color_update(k::Float32) = RGBA{U8}[RGBA{U8}(
		cos_01((i/2f0-sin(k)))*sin_01((j/2f0+cos(k))) * abs(sin(k)),
		cos_01(j/2f0)*sin_01(i/2f0) * abs(cos(k)),
		0.01f0,
		cap((cos_01((i/2f0-sin(k)))*sin_01((j/2f0+cos(k))) * abs(sin(k))) +
		(cos_01(j/2f0)*sin_01(i/2f0) * abs(cos(k))))
	) for i=0f0:0.1f0:30f0, j=0f0:0.1f0:2f0
]
robj = visualize(lift(color_update, bounce(0f0:0.01f0:100f0)));
view(robj);
bb = robj.boundingbox.value;

view(visualize("random access"), model=translationmatrix(-bb.maximum-bb.minimum));


text:
cd("/Volumes/SSD Hans/Github/MacroFall2020/AdvMacro/Assignment5")
# mkpath("Figures")
using Plots
# using LateXStrings # Pkg.add("LaTeXStrings") # https://github.com/stevengj/LaTeXStrings.jl
using Dierckx # Pkg.add("Dierckx") # https://github.com/kbarbary/Dierckx.jl
using Interpolations # Pkg.add("Interpolations") # https://github.com/JuliaMath/Interpolations.jl
using ForwardDiff # Pkg.add("ForwardDiff") # https://github.com/JuliaDiff/ForwardDiff.jl
using Optim # Pkg.add("Optim") # https://julianlsolvers.github.io/Optim.jl/stable/
    using Optim: converged, maximum, maximizer, minimizer, iterations
using Roots # Pkg.add("Roots") # https://github.com/JuliaMath/Roots.jl
using Parameters # Pkg.add("Parameters") # https://github.com/mauro3/Parameters.jl
using Distributions #Pkg.add("Distributions")
using QuadGK # Pkg.add("QuadGK") # https://juliamath.github.io/QuadGK.jl/latest/
using LinearAlgebra
using Random
using Statistics
using Latexify
Pkg.add("StatsPlots")
using StatsPlots
# Call Scaled Interpolation Functions
    include("../Scaled_Interpolation_Functions.jl")

# Set random seed
Random.seed!(0021);

# Paramters
    # Generate structure for parameters using Parameters module
    # Set default values for parameters
    @with_kw struct Par
        # Model Parameters
        ρ::Float64 = 0.9    ; # AR(1) persistence
        σ::Float64 = 2.0    ; # Standard deviation of inovations η
        n_draw::Int64 = 10000 ; # Number of draws for simulated markov chain
    end
# Allocate parameters to object p for future calling
p = Par()

# Function to generate markov chain z using monte carlo simulation
function Simulate_MC(p::Par)
    @unpack ρ,σ,n_draw = p
    # Generate inovations η
    η = rand(Normal(0.0,σ),n_draw)
    # initialize process for z
    z = zeros(n_draw)
    # Compute z according to AR(1) process
    z[1] = η[1]
    for i in 2:n_draw
        z[i] = ρ*z[i-1] + η[i]
    end
    return z
end

# Function to distretize AR(1) markov process with Tauchen (1986)
function Tauchen86(N,p::Par,Ω::Any=3)
    @unpack ρ,σ = p
    # INPUTS:
        # ρ: persistence of unerlying AR(1) process where z' = ρz+η
        # σ: Std dev of inovation η in AR(1) process where η∼N(0,σ^2)
        # N: Size of grid for discrete process
        # Ω: Denotes boundaries of the grid as the number of standard deviations below and above the mean (lb=-Ωσ,ub=Ωσ)
    # OUTPUT:
        # z: All possible values of discretized AR(1) process, equally spaced grid of size N
        # Π: Matrix of transition probabilities
        # PDF_z: Stationary PDF of z
    #--------------------------------------------------------------------------
     z = collect(range(-Ω*σ,Ω*σ;length=N))
     Δz= 2*Ω*σ/(N-1)
     # initialize and compute probability transition matrix
     Π = zeros(N,N)
     F(x) = cdf.(Normal(),x)
     function prob_ij(zzp,zz,z::Vector)
         if zzp == z[1]
             return F((zzp+Δz/2-ρ*zz)/σ)
         elseif zzp == z[N]
             return 1-F((zzp-Δz/2-ρ*zz)/σ)
         else
             return F((zzp+Δz/2-ρ*zz)/σ)-F((zzp-Δz/2-ρ*zz)/σ)
         end
     end
     for i in 1:N
         for j in 1:N
            Π[i,j] = prob_ij(z[j],z[i],z)
        end
     end
     # Stationary distribution
     PDF_z = real(eigvecs(Π')[:,end]); PDF_z = PDF_z/sum(PDF_z)
     return (z,Π,PDF_z)
end

# Function to distretize AR(1) markov process with Rouwenhorst (1995)
function Rouwenhorst95(N,p::Par)
    @unpack ρ,σ=p
    # INPUTS:
        # ρ: persistence of unerlying AR(1) process where z' = ρz+η
        # σ: Std dev of inovation η in AR(1) process where η∼N(0,σ^2)
        # N: Size of grid for discrete process
    # OUTPUT:
        # z: All possible values of discretized AR(1) process, equally spaced grid of size N
        # Π: Matrix of transition probabilities
        # PDF_z: Stationary PDF of z
    #---------------------------------------------------------------------------
    Π = zeros(N,N)
    Π_Nm = zeros(N-1,N-1)
    P = (1+ρ)/2
    ϕ = σ*(sqrt((N-1)/(1-ρ^2)))
    z = range(-ϕ,ϕ;length=N)
    if N==2
        Π = [P 1-P;1-P P]
    else
        Π_Nm = Rouwenhorst95(N-1,p)[2]
        o = zeros(N-1)
        Π = P*[Π_Nm o; o' 0] + (1-P)*[o Π_Nm; 0 o'] + (1-P)*[o' 0; Π_Nm o] + P*[0 o';o Π_Nm]
        Π = Π./repeat(sum(Π,dims=2),1,N)
    end
    PDF_z = pdf.(Binomial(N-1,0.5),(0:N-1))
    return (z,Π,PDF_z)
end

# Function to simulate a discrete markov process
function Simulate_discrete(N,z,PDF_s,Π,p::Par)
    @unpack ρ,σ,n_draw = p
    # INPUT:
        # n_draw: Number of draws for the simulation
        # N     : Number of possible values for the discrete process
        # z     : Vector of possible values for the discrete process
        # PDF_s : Stationary PDF (used to draw initial point)
        # Π     : Transition probability matrix
    # OUTPUT:
        # z_MC  : Simulated draws
    # Initialization
    z_ind = zeros(Int64,n_draw) # index values of z
    z_MC = zeros(n_draw)
    # Initial value
    z_ind[1] = rand(Categorical(PDF_s))
    z_MC[1] = z[z_ind[1]]
    # Simulation
    for i in 2:n_draw
        z_ind[i] = rand(Categorical(Π[z_ind[i-1],:]))
        z_MC[i] = z[z_ind[i]]
    end
    return z_MC
end

# Function that calculate numerical moments out of a given array
function num_moments(z)
    # First four moments
    num_mean = mean(z)
    num_var = var(z)
    num_skew = skewness(z)
    num_kurt = kurtosis(z)
    # Autocorrelation going for t-1 to t-4
    autocorr = zeros(4)
    for i in 1:4
        autocorr[i] = cor(z[(i+1):end],z[1:(end-i)])
    end
    return (num_mean,num_var,num_skew,num_kurt,autocorr)
end

                    ### a) ###
# Simulate markov chain
z = Simulate_MC(p)

                    ### b) ###
N = [5,15]
# Tauchen (1986)
ztauc_MC = zeros(p.n_draw,2)
for i in 1:2
    (ztauc,Π_tauc,PDF_ztauc) = Tauchen86(N[i],p)
    ztauc_MC[:,i] = Simulate_discrete(N[i],ztauc,PDF_ztauc,Π_tauc,p)
end
# Rouwenhorst (1995)
zrouw_MC = zeros(p.n_draw,2)
for i in 1:2
    (zrouw,Π_rouw,PDF_zrouw) = Rouwenhorst95(N[i],p)
    zrouw_MC[:,i] = Simulate_discrete(N[i],zrouw,PDF_zrouw,Π_rouw,p)
end

                    ### c) ###
# Report moments for simulations, Ns=10000
tmean,rmean,tvar,rvar,tskew,rskew,tkur,rkur = [zeros(2) for i in 1:8]
mcmean,mcvar,mcskew,mckur = [zeros(1) for i in 1:4]
tautocorr,rautocorr = [zeros(4,2) for i in 1:2]
mcautocorr = zeros(4)
for i in 1:2
    (tmean[i],tvar[i],tskew[i],tkur[i],tautocorr[:,i]) = num_moments(ztauc_MC[:,i])
    (rmean[i],rvar[i],rskew[i],rkur[i],rautocorr[:,i]) = num_moments(zrouw_MC[:,i])
end
(mcmean,mcvar,mcskew,mckur,mcautocorr) = num_moments(z)
r(x) = round(x,digits=3)
Moments_Mat = ["" "Data" "Tauchen" "Tauchen" "" "Rouwenhorst" "Rouwenhorst";
               "" " " "N=5" "N=15" "" "N=5" "N=15" ;
               "mean" r(mcmean) r(tmean[1]) r(tmean[2]) "" r(rmean[1]) r(rmean[2]);
               "var"  r(mcvar)  r(tvar[1]) r(tvar[2]) "" r(rvar[1]) r(rvar[2]);
               "skew"  r(mcskew)  r(tskew[1]) r(tskew[2]) "" r(rskew[1]) r(rskew[2]);
               "kurt" r(mckur) r(tkur[1]) r(tkur[2]) "" r(rkur[1]) r(rkur[2]);
               "acorr" r(mcautocorr[1]) r(tautocorr[1,1]) r(tautocorr[1,2]) "" r(rautocorr[1,1]) r(rautocorr[1,2])]

println("\n Moments from simulated paths: \n")
display(Moments_Mat)

copy_to_clipboard(true)
mdtable(Moments_Mat)

#-------------------------------------------------------
# Plot paths for simulations, Ns=100
gr()
plot(title="Tauchen Simulation Paths",legend=:outerright,foreground_color_legend = nothing,background_color_legend = nothing)
plot!(1:100,z[1:100], linewidth=3,label="Simulation")
plot!(1:100,ztauc_MC[1:100,1] ,linewidth=3,label="N=5",linestyle=(:dot) )
plot!(1:100,ztauc_MC[1:100,2],linewidth=3,label="N=15",linestyle=(:dash))
xlabel!("Time")
savefig("./Figures/Tauchen")

gr()
plot(title="Rouwenhorst Simulation Paths",legend=:outerright,foreground_color_legend = nothing,background_color_legend = nothing)
plot!(1:100,R_s_5[1:100] ,linewidth=3,label="N=5" )
plot!(1:100,R_s_11[1:100],linewidth=3,label="N=11",linestyle=(:dash))
plot!(1:100,R_s_21[1:100],linewidth=3,label="N=21",linestyle=(:dot))
xlabel!("Time")
savefig("./Figures/Rouwenhorst")])

density(z, label="Simulation", linewidth=3,)
density!(ztauc_MC[:,1], label="N=5", linewidth=3,)
density!(ztauc_MC[:,2], label="N=15", linewidth=3,)
title!("Tauchen")
savefig("./Figures/Tauchen")

density(z, label="Simulation", linewidth=3,)
density!(zrouw_MC[:,1], label="N=5", linewidth=3,)
density!(zrouw_MC[:,2], label="N=15", linewidth=3,)
title!("Rouwenhorst")
savefig("./Figures/Rouwen")

max_stars_repo_name: py-r-hans/AdvMacro

max_issues_repo_name: py-r-hans/AdvMacro

max_forks_repo_name: py-r-hans/AdvMacro


text:
using TestImages

# load mandril from testimages and save it to disk
img = testimage("mandril_color");
save("mandril_color.png", img);