text:
@testset "recv_connect.jl" begin
    @testset begin receiver = TCPReceiver(4000)
        @test window_size(receiver) == 4000
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        segment_arrives(receiver; with_syn=true, result=SegmentArrives_Result_OK)
        @test ackno(receiver).val == 1
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
    end

    @testset begin receiver = TCPReceiver(5435)
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
        @test window_size(receiver) == 5435

        segment_arrives(receiver; seqno=893475, result=SegmentArrives_Result_NOT_SYN)
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
    end

    @testset begin receiver = TCPReceiver(5435)
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
        @test window_size(receiver) == 5435

        segment_arrives(receiver; with_fin=true, with_ack=0, seqno=893475, result=SegmentArrives_Result_NOT_SYN)
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
    end

    @testset begin receiver = TCPReceiver(5435)
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
        @test window_size(receiver) == 5435

        segment_arrives(receiver; with_fin=true, with_ack=0, seqno=893475, result=SegmentArrives_Result_NOT_SYN)
        @test ackno(receiver) === nothing
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; with_syn=true, seqno=89347598, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(89347599)
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
    end

    @testset begin receiver = TCPReceiver(4000)
        segment_arrives(receiver; with_syn=true, seqno=5, with_fin=true, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(7)
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 0
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
module SortedVectors

export SortedVector

using ArgCheck: @argcheck
using DocStringExtensions: SIGNATURES
using Lazy: @forward
using Parameters: @unpack

"""
Flag for indicating that

1. sorting should be verified,
2. the argument vector will not be modified later.

Not exported.
"""
struct CheckSorted end

"""
Flag for indicating that

1. the argument should be assumed to be sorted, and **this should not be checked**,
2. the argument vector will not be modified later.

Not exported.
"""
struct AssumeSorted end

struct SortedVector{T, F, V <: AbstractVector{T}} <: AbstractVector{T}
    "comparison function for sorting"
    lt::F
    "sorted contents"
    sorted_contents::V
    function SortedVector(::AssumeSorted, lt::F,
                          sorted_contents::V) where {F, V <: AbstractVector}
        new{eltype(sorted_contents), F, V}(lt, sorted_contents)
    end
end

function SortedVector(::CheckSorted, lt, sorted_contents::AbstractVector)
    @argcheck issorted(sorted_contents; lt = lt)
    SortedVector(AssumeSorted(), lt, sorted_contents)
end

"""
    SortedVector([lt=isless], xs)

Sort `xs` by `lt` (which defaults to `isless`) and wrap in a SortedVector. For reverse
sorting, use `!lt`.

    SortedVector(SortedVectors.CheckSorted(), lt, sorted_contents)

Checks that the vector is sorted, throws an `ArgumentError` if it isn't. This is a
relatively cheap operation if the vector is supposed to be sorted but this should be
checked. Caller should `copy` the `sorted_contents` if they are mutable and may be modified.

    SortedVector(SortedVectors.AssumeSorted(), lt, sorted_contents)

Unchecked, unsafe constructor. Use only if you are certain that `sorted_contents` is sorted
according to `lt`, otherwise results are undefined. `copy` the `sorted_contents` if they are
mutable and may be modified.
"""
function SortedVector(lt, xs::AbstractVector)
    SortedVector(AssumeSorted(), lt, sort(xs; lt = lt))
end

SortedVector(xs::AbstractVector) = SortedVector(isless, xs)

Base.parent(sv::SortedVector) = sv.sorted_contents

####
#### array interface
####

@forward SortedVector.sorted_contents (Base.size, Base.getindex, Base.length, Base.axes)

Base.IndexStyle(::Type{<:SortedVector}) = Base.IndexLinear()

function Base.setindex!(sv::SortedVector, x, i::Integer)
    @unpack lt, sorted_contents = sv
    a, b = firstindex(sorted_contents), lastindex(sorted_contents)
    a < i ≤ b && @argcheck lt(sorted_contents[i-1], x)
    a ≤ i < b && @argcheck lt(x, sorted_contents[i+1])
    sorted_contents[i] = x
end

####
#### cut
####

Base.searchsortedfirst(sv::SortedVector, x) =
    searchsortedfirst(sv.sorted_contents, x; lt = sv.lt)

"""
$(SIGNATURES)

Return `i` such that `breaks[i] < x ≤ breaks[i + 1]`, where `<` is the sorting from the
second argument.

When `x ≤ breaks[1]`, `left` (0 by default) is used, and `right` (default: length) for `x >
breaks[end]`.
"""
function cut(x, breaks::SortedVector; left = 0, right = length(breaks))
    ix = searchsortedfirst(breaks, x) - 1
    ix == 0 ? left : (ix == length(breaks) ? right : ix)
end

end # module

ext: jl


text:
function make_script()
	sample_size = collect(10000:10000:120000)
	for samples in sample_size
		open("hi.sh", "w") do file
    		write(file, "#!/bin/bash")
		end
	end
end


text:
mutable struct P_filt{T<:Real}
	fs::Float64
	freq::Vector{Float64}
	p_conv::Conv.P_conv{T,1,1,1}
end

function P_filt(nt::Int=128, filt_func=x->butterworth(x, order=2, fmin=0.0, fmax=1.0); fs=1)
	# create a plan for conv
	p_conv=Conv.P_conv(gsize=[nt], dsize=[nt], ssize=[nt],)

	freq=rfftfreq(p_conv.np2, fs)
	filt=filt_func(freq)
	for i in eachindex(p_conv.sfreq)
		p_conv.sfreq[i] = filt[i]
	end

	A_mul_B!(p_conv.spad, p_conv.sifftp, p_conv.sfreq)
	Conv.pad_truncate!(p_conv.s, p_conv.spad, p_conv.slags[1],
		    p_conv.slags[2], p_conv.np2, -1)

	fs=Float64(fs)
	return P_filt(fs, Vector(freq), p_conv)

end

function filter!(xx, x, pa::P_filt)

	for ii in 1:size(xx,2)
		xv=view(x, :, ii)
		xxv=view(xx, :, ii)
		Conv.mod!(pa.p_conv, :d, d=xxv, g=xv)
	end
end

"""
A bandpass butterworth filter using `fmin` and `fmax` in the frequency domain.
Return either a zero-phase or minimum-phase filter using `attrib`.
More info on minimum-phase filtering: http://www.katjaas.nl/minimumphase/minimumphase.html.
! then it is converted to minimum phase filter in the cepstral domain
! more info: http://www.katjaas.nl/minimumphase/minimumphase.html
! positive quefrencies correspond to minimum phase component of the signal
! negetive quefrencies correspond to maximum phase component of the signal
! signal = minimum phase [convolved with] maximum phase

# Arguments

* ``

# Keyword Arguments

* `order::Int64` : order of the butterworth filter
* `attrib::Symbol`
  * `=:zp` means zero phase
  * `=:mp` means zero phase
* `fmin::Float64`
  * `=0.` means a lowpass filter at fmax
* `fmax::Float64`
  * `=0.` means a highpass filter at fmin

"""
function butterworth(fvec; order::Int64=2, attrib::Symbol=:zp, fmin=0.0, fmax=0.0)
# Author : Pawan Bharadwaj
#          p.b.pisupati@tudelft.nl
# August 2017, imported to Julia from FORTRAN90
	(fmin < 0.0) && error("fmin cannot be .lt. zero")
	(fmax < 0.0) && error("fmax cannot be .lt. zero")
	(fmax ≠ 0.0) && (fmax <= fmin) && error("fmax .le. fmin")

	if(fmin == 0.0)
		# lowpass filter
		F  = complex.((1.0 + (fvec./fmax).^(2.0*order)).^(-1.0), 0.0)
	elseif(fmax == 0.0)
		# highpass filter
		F  = complex.((1.0 + (fmin./fvec).^(2.0*order)).^(-1.0), 0.0)
	elseif((fmin ≠ 0.0) & (fmax ≠ 0.0))
		# bandpass filter
		x1 = sqrt(fmin*fmax) / (fmax-fmin) * (fvec./sqrt(fmin*fmax) + sqrt(fmax*fmin)./fvec);
		F  = complex.((1.0 + (x1).^(2.0*order)).^(-1.0), 0.0)
	end
	# DC component is always forced zero
	F[1] = complex(0.0, 0.0);

	# conversion to minimum phase
	if(attrib == :mp)
	# to prevent log(0)
		damp = 1e-20 * maximum(abs.(F))
		# logarithm
		X = log.(complex.(abs.(F) + damp, 0.0))
		# to cepstral domain - IFFT
		ifft!(X)
		# only real part output
		X = complex.(real.(X), 0.0);
		# scaling
		# X = X / complex(real(npow2), 0.0)

		# positive cepstrum x 2
		X[2 : npow2/2 + 1] *= complex(2.0, 0.0)
		# remove negative quefrencies
		X[npow2/2 + 2 : npow2] = complex(0.0, 0.0)

		# FFT
		fft!(X)

		# exponential
		F = exp.(X)

		F /= complex(maximum(abs.(F)), 0.0)
	end

	return F

end # butterworth


text:
#----- API -----#
trace(results::SimulationResults) = length(results.trace) > 0 ? results.trace : error("No trace in optimization results. To get a trace, run simulation with store_trace = true.")
prices(results::SimulationResults) = results.p̄
#---------------------------------------------------------------------------------------------------

max_stars_repo_name: IvanJericevich/FinancialModels

max_issues_repo_name: IvanJericevich/FinancialModels

max_forks_repo_name: IvanJericevich/FinancialModels


text:
#
# Address, Person types for deep match test
#

type Address
    street::AbstractString
    city::AbstractString
    zip::AbstractString
end

type Person
    firstname::AbstractString
    lastname::AbstractString
    address::Address
end

#
# Untyped lambda calculus definitions
#

abstract Term

immutable Var <: Term
    name::AbstractString
end

immutable Fun <: Term
    arg::AbstractString
    body::Term
end

immutable App <: Term
    f::Term
    v::Term
end


text:
using WK3

function main()
    # option flags
    rstflag = "no"

    # initialization
    numbeats = 30; # total number of cardiac cycles
    ensemblesize = 50;
    if rstflag == "no"
        y0 = [[rand(Distributions.Normal(10,1))*WK3.mmHgToPa/WK3.Ps;
            rand(Distributions.Normal(135,10))*WK3.cm3Tom3/WK3.Vs;0*WK3.cm3Tom3/WK3.Qs] for i=1:ensemblesize]; # random
        # y0 = [[10*WK3.mmHgToPa/WK3.Ps;135*WK3.cm3Tom3/WK3.Vs;0*WK3.cm3Tom3/WK3.Qs] for i=1:ensemblesize]; # uniform
        nvar = [3 for i=1:ensemblesize];
        systems = pmap((a1)->WK3.CVSystem(a1),nvar);
    elseif rstflag == "yes"
        fnames = ["system_$i.mat" for i=1:ensemblesize];
        vars = [MAT.matread(fnames[i]) for i=1:ensemblesize];
        nvar = [3 for i=1:ensemblesize];
        rst = [rstflag for i=1:ensemblesize];
        systems = pmap((a1,a2,a3)->WK3.CVSystem(a1,a2,a3),nvar,vars,rst);
        y0 = [[vars[1]["system"]["Pa"][end];1;vars[1]["system"]["Q"][end]] for i=1:ensemblesize];
    end
    t0 = [0/WK3.ts for i=1:ensemblesize];
    sparams = [systems[i].sparams for i=1:ensemblesize];
    mparams = [systems[i].mparams for i=1:ensemblesize];
    tm = linspace(0.,systems[1].mparams.th[1],1e2);
    g1 = zeros(length(tm));
    g2 = zeros(length(tm));
    h1 = zeros(length(tm));
    h2 = zeros(length(tm));
    to = Vector{Float64}[];
    yo = Vector{Vector{Float64}}[];
    ke = zeros(ensemblesize);
    append!(yo,[[zeros(1)]]);
    append!(to,[zeros(1)]);

    # data assimilation setup
    # load patient data
    filename = "patient.mat";
    tf,qdata,vdata = patdatainterp(systems[1].mparams.th[1],filename);

    # allocators for ensemble augmented state, measurements
    nparams = 9;
    X = [zeros(nvar[1]) for i in (1:length(systems))];
    θ = [zeros(nparams) for i in (1:length(systems))];
    Y = [zeros(2) for i in (1:length(systems))];

    # parameter scalings
    θs = zeros(nparams);
    for i = 1:length(systems)
        θs[1] += mparams[i].Zc;
        θs[2] += mparams[i].R;
        θs[3] += mparams[i].C;
        θs[4] += mparams[i].V0;
        θs[5] += mparams[i].t1;
        θs[6] += mparams[i].t2;
        θs[7] += mparams[i].m1;
        θs[8] += mparams[i].m2;
        θs[9] += mparams[i].Emax;
    end
    θs /= ensemblesize;

    # allocators for state, forecast measurement mean
    xhat = zeros(nvar[1]);
    θhat = zeros(nparams);
    yhat = zeros(2);

    # allocators for measurement replicates
    yi = [zeros(2) for i in 1:length(systems)];

    # output variables for ensemble avg. state, meas., params. w/ times
    tout = Float64[];
    xout = Vector{Float64}[];
    xoutv = Vector{Float64}[];
    yout = Vector{Float64}[];
    youtv = Vector{Float64}[];
    innov = Vector{Float64}[];
    ioutv = Vector{Float64}[];
    θout = Vector{Float64}[];
    θoutv = Vector{Float64}[];
    Pθout = Vector{Float64}[];
    Pθoutv = Vector{Float64}[];
    Pxout = Vector{Float64}[];
    Pxoutv = Vector{Float64}[];
    lbx = Vector{Float64}[];
    ubx = Vector{Float64}[];
    lbxv = Vector{Float64}[];
    ubxv = Vector{Float64}[];
    Pvout = Float64[];
    numassims = 0;
    println("Number of assimilations: $(length(tf))")

    # distributions for state, measurement, parameters
    error = WK3.Errors(nparams);

    # allocators for calculating normalized RMSE ratio
    r1dot = 0;
    r2dot = zeros(ensemblesize);

    # RTPS allocators
    p = 0.5; # relaxation amount, ∃ [0.5, 0.95]
    c = zeros(length(systems));
    σb = zeros(length(xhat));
    σa = zeros(length(xhat));
    σtb = zeros(length(θhat));
    σta = zeros(length(θhat));

    # solver loop
    tic()
    for n in 1:numbeats
        for nn in 1:length(tf)+1
            # randomize initial ventricular volume if starting new cycle
            if nn == 1 && rstflag == "no"
                for i=1:ensemblesize
                    y0[i][2] = rand(Distributions.Normal(1,0.1)); # random
                    # y0[i][2] = 1.; # uniform (for debugging)
                end
            end
            # normal time integration
            if nn <= length(tf)
                tfcur = [(n-1)*0.8/WK3.ts + tf[nn]/WK3.ts for i=1:ensemblesize];
            elseif nn > length(tf)
                tfcur = [n*0.8/WK3.ts for i=1:ensemblesize];
            end
            # calculate elastance scaling
            for i = 1:ensemblesize
                g1 .= (tm./mparams[i].t1).^mparams[i].m1;
                g2 .= (tm./mparams[i].t2).^mparams[i].m2;
                h1 .= g1./(1.+g1);
                h2 .= 1./(1.+g2);
                ke[i] .= (mparams[i].Emax.-mparams[i].Emin)./maximum(h1.*h2);
            end
            # collect outputs
            # println("Time integration:")
            soln = pmap((a1,a2,a3,a4,a5,a6)->WK3.odeint(a1,a2,a3,a4,a5,a6),y0,t0,tfcur,sparams,mparams,ke);
            to = [soln[i][1] for i in 1:length(soln)];
            yo = [soln[i][2] for i in 1:length(soln)];
            for i=1:ensemblesize
                append!(systems[i].t,to[i])
            end

            # diagnostic variables
            for i=1:ensemblesize
                g1 = (tm/systems[i].mparams.t1).^systems[i].mparams.m1;
                g2 = (tm/systems[i].mparams.t2).^systems[i].mparams.m2;
                h1 = g1./(1+g1);
                h2 = 1./(1+g2);
                k = (systems[i].mparams.Emax-systems[i].mparams.Emin)/maximum(h1.*h2);
                for j = 1:length(to[i])
                    g1t = (mod(to[i][j]*WK3.ts,systems[i].mparams.th[end])/systems[i].mparams.t1).^systems[i].mparams.m1;
                    g2t = (mod(to[i][j]*WK3.ts,systems[i].mparams.th[end])/systems[i].mparams.t2).^systems[i].mparams.m2;
                    h1t = g1t/(1+g1t);
                    h2t = 1/(1+g2t);
                    push!(systems[i].E,k*h1t*h2t+systems[i].mparams.Emin);
                    push!(systems[i].Pv,systems[i].E[end]*(yo[i][j][2]*WK3.Vs-systems[i].mparams.V0)/WK3.Ps);
                end
            end

            # data assimilation
            if nn <= length(tf)
                println("Assimilating data at t = $(tfcur[1])")

                # get last time step taken, use as time step for assimilation
                if tfcur[1] != t0[1] # ignore if measurement is at t0
                    # warn("Assimilating data. Fixing time step.")
                    for i in 1:length(soln)
                        sparams[i].h0 = to[i][end] - to[i][end-1];
                    end
                end

                # non-dimensional measurement
                y = [vdata[nn]/WK3.Vs,qdata[nn]/WK3.Qs];
                # y = [pdata[nn]/WK3.Ps,vdata[nn]/WK3.Vs,qdata[nn]/WK3.Qs];

                # generate measurement replicates
                for i = 1:length(systems)
                    # yi[i][1] = y[1] + rand(Distributions.Normal(0,error.odev[1]/WK3.Ps));
                    # yi[i][2] = y[2] + rand(Distributions.Normal(0,error.odev[2]/WK3.Vs));
                    # yi[i][3] = y[3] + rand(Distributions.Normal(0,error.odev[3]/WK3.Qs));
                    yi[i][1] = y[1] + rand(Distributions.Normal(0,error.odev[2]/WK3.Vs));
                    yi[i][2] = y[2] + rand(Distributions.Normal(0,error.odev[3]/WK3.Qs));
                end

                # forecast parameters and means
                for i = 1:length(systems)
                    θ[i][1] = mparams[i].Zc;
                    θ[i][2] = mparams[i].R;
                    θ[i][3] = mparams[i].C;
                    θ[i][4] = mparams[i].V0;
                    θ[i][5] = mparams[i].t1;
                    θ[i][6] = mparams[i].t2;
                    θ[i][7] = mparams[i].m1;
                    θ[i][8] = mparams[i].m2;
                    θ[i][9] = mparams[i].Emax;
                end

                # forecast parameters (dimensional)
                # println("Forecast parameters:")
                θ = [WK3.paramwalk!(error,θs,θ[i]) for i in 1:length(systems)];

                # non-dimensionalize parameters
                for i = 1:length(systems)
                    θ[i] = θ[i]./θs;
                end

                # RTPS prior standard deviation
                for i in 1:length(θ[1])
                    for j in 1:length(systems)
                        c[j] = θ[j][i];
                    end
                    σtb[i] = std(c;corrected=true);
                end

                # parameters back into model
                for i = 1:length(soln)
                    # mparams[i].Zc = θ[i][1]*θs[1];
                    # mparams[i].R = θ[i][2]*θs[2];
                    # mparams[i].C = θ[i][3]*θs[3];
                    # mparams[i].V0 = θ[i][4]*θs[4];
                    # mparams[i].t1 = θ[i][5]*θs[5];
                    # mparams[i].t2 = θ[i][6]*θs[6];
                    # mparams[i].m1 = θ[i][7]*θs[7];
                    # mparams[i].m2 = θ[i][8]*θs[8];
                    # mparams[i].Emax = θ[i][9]*θs[9];
                    if mparams[i].Emax <= mparams[i].Emin
                        mparams[i].Emin = mparams[i].Emax - 1e5;
                    end
                end

                # calculate elastance scaling
                for i = 1:ensemblesize
                    g1 .= (tm./mparams[i].t1).^mparams[i].m1;
                    g2 .= (tm./mparams[i].t2).^mparams[i].m2;
                    h1 .= g1./(1.+g1);
                    h2 .= 1./(1.+g2);
                    ke[i] .= (mparams[i].Emax.-mparams[i].Emin)./maximum(h1.*h2);
                end

                # forecast state w/ forecast parameters (single time step)
                soln = pmap((a1,a2,a3,a4,a5,a6)->WK3.odeint(a1,a2,a3,a4,a5,a6),y0,t0,tfcur,sparams,mparams,ke);
                to = [soln[i][1] for i in 1:length(soln)];
                yoa = [soln[i][2] for i in 1:length(soln)];

                # vector of forecast state vectors
                X = [yoa[i][end][1:end] for i in 1:length(soln)];

                # vector of forecast measurements
                Y = [[yoa[i][end][2],yoa[i][end][3]] for i in 1:length(soln)]; # V, flowrate

                # forecast mean state, parameters, measurement
                xhat = mean(X);
                yhat = mean(Y);
                θhat = mean(θ);
                # println("Normalized ̂x, first forecast: $xhat")
                # println("Normalized ̂y, first forecast: $yhat")
                # println("Normalized ̂θ, first forecast: $θhat")

                # forecast params./meas. cross covariance, meas. covariance
                Pty = zeros(length(θhat),length(yhat))
                Pyy = zeros(length(yhat),length(yhat))
                for i = 1:length(soln)
                    Pty += *((θ[i] .- θhat),(Y[i] .- yhat)');
                    Pyy += *((Y[i] .- yhat),(Y[i] .- yhat)');
                end
                Pty ./= (ensemblesize);
                Pyy ./= (ensemblesize);

                # add meas. noise to meas. covariance (allows invertibility)
                Pyy += diagm([(error.odev[2]/WK3.Vs)^2,(error.odev[3]/WK3.Qs)^2],0);

                # parameter Kalman gain
                K = Pty*inv(Pyy);
                # println("Parameter Kalman gain:")
                # display(K)

                # parameter analysis step
                for i = 1:length(soln)
                    # println("ith measurement replicate: $(yi[i])")
                    # println("ith forecast measurement: $(Y[i])")
                    θ[i][:] += K*(yi[i] .- Y[i]);
                end

                # println("Normalized θ after analysis: $(mean(θ))")

                # RTPS parameter covariance inflation
                for i in 1:length(θ[1])
                    for j in 1:length(systems)
                        c[j] = θ[j][i];
                    end
                    σta[i] = std(c;corrected=true);
                end
                θhat = mean(θ);
                for i = 1:length(soln)
                    θ[i] .= θ[i] .+ p.*((σtb.-σta)./σta).*(θ[i].-θhat);
                end

                # println("Normalized θ after RTPS: $(mean(θ))")

                # # compute FD Jacobian of observation vector w.r.t. parameters (per realization)
                # yvec = [[Y[i]] for i in 1:length(soln)];
                # J = pmap((a1,a2,a3,a4,a5,a6)->WK3.fdjac(a1,a2,a3,a4,a5,a6),y0,t0,tfcur,sparams,mparams,yvec);
                #
                # # non-dimensionalize Jacobian
                # for i in 1:length(soln) # i indexes ensemble members
                #     for j in 1:size(J[1],1) # j indexes observations
                #         J[i][j,:] .*= θs;
                #     end
                # end
                #
                # # compute dimensionless FIM (per realization)
                # I = [J[i]'*((WK3.Qs/error.odev[1])^2)*J[i] for i in 1:length(soln)]
                # println("Minimum magnitude of I[1]: $(minimum(abs.(I[1])))")
                # println("Maximum magnitude of I[1]: $(maximum(abs.(I[1])))")
                # println("Minimum rank of FIM: $(minimum(rank.(I)))")

                # analysis parameters back into ensemble members
                for i = 1:length(soln)
                    # mparams[i].Zc = θ[i][1]*θs[1];
                    # mparams[i].R = θ[i][2]*θs[2];
                    # mparams[i].C = θ[i][3]*θs[3];
                    # mparams[i].V0 = θ[i][4]*θs[4];
                    # mparams[i].t1 = θ[i][5]*θs[5];
                    # mparams[i].t2 = θ[i][6]*θs[6];
                    # mparams[i].m1 = θ[i][7]*θs[7];
                    # mparams[i].m2 = θ[i][8]*θs[8];
                    # mparams[i].Emax = θ[i][9]*θs[9];
                    # if mparams[i].Emax <= mparams[i].Emin
                    #     mparams[i].Emin = mparams[i].Emax - 1e5;
                    # end
                    systems[i].mparams = mparams[i];
                    if mparams[i].Zc <= error.lb[1]
                        println("Warning: analysis Zc below lower bound for member $i.
                            Setting to lower bound of $(error.lb[1]).")
                        mparams[i].Zc = error.lb[1];
                        θ[i][1] = error.lb[1]/θs[1];
                    end
                    if mparams[i].R <= error.lb[2]
                        println("Warning: analysis R below lower bound for member $i.
                            Setting to lower bound of $(error.lb[2]).")
                        mparams[i].R = error.lb[2];
                        θ[i][2] = error.lb[2]/θs[2];
                    end
                    if mparams[i].C <= error.lb[3]
                        println("Warning: analysis C below lower bound for member $i.
                            Setting to lower bound of $(error.lb[3]).")
                        mparams[i].C = error.lb[3];
                        θ[i][3] = error.lb[3]/θs[3];
                    end
                    if mparams[i].V0 <= error.lb[4]
                        println("Warning: analysis V0 below lower bound for member $i.
                            Setting to lower bound of $(error.lb[4]).")
                        mparams[i].V0 = error.lb[4];
                        θ[i][4] = error.lb[4]/θs[4];
                    end
                    if mparams[i].t1 <= error.lb[5]
                        println("Warning: analysis τ1 below lower bound for member $i.
                            Setting to lower bound of $(error.lb[5]).")
                        mparams[i].t1 = error.lb[5];
                        θ[i][5] = error.lb[5]/θs[5];
                    end
                    if mparams[i].t2 <= error.lb[6]
                        println("Warning: analysis τ2 below lower bound for member $i.
                            Setting to lower bound of $(error.lb[6]).")
                        mparams[i].t2 = error.lb[6];
                        θ[i][6] = error.lb[6]/θs[6];
                    end
                    if mparams[i].m1 <= error.lb[7]
                        println("Warning: analysis m1 below lower bound for member $i.
                            Setting to lower bound of $(error.lb[7]).")
                        mparams[i].m1 = error.lb[7];
                        θ[i][7] = error.lb[7]/θs[7];
                    end
                    if mparams[i].m2 <= error.lb[8]
                        println("Warning: analysis m2 below lower bound for member $i.
                            Setting to lower bound of $(error.lb[8]).")
                        mparams[i].m2 = error.lb[8];
                        θ[i][8] = error.lb[8]/θs[8];
                    end
                    if mparams[i].Emax <= error.lb[9]
                        println("Warning: analysis Emax below lower bound for member $i.
                            Setting to lower bound of $(error.lb[9]).")
                        mparams[i].Emax = error.lb[9];
                        θ[i][9] = error.lb[9]/θs[9];
                    end
                end

                # recalculate and output parameter averages (non-dimensional)
                θhat = mean(θ);
                append!(θout,[θhat])
                # println("Normalized analysis parameter averages: $θhat")

                # analysis parameter variance (for post-processing)
                Ptt = zeros(length(θhat))
                for i = 1:length(soln)
                    Ptt += (θ[i] .- θhat).^2;
                end
                Ptt ./= ensemblesize;
                append!(Pθout,[Ptt])

                # calculate elastance scaling
                for i = 1:ensemblesize
                    g1 .= (tm./mparams[i].t1).^mparams[i].m1;
                    g2 .= (tm./mparams[i].t2).^mparams[i].m2;
                    h1 .= g1./(1.+g1);
                    h2 .= 1./(1.+g2);
                    ke[i] .= (mparams[i].Emax.-mparams[i].Emin)./maximum(h1.*h2);
                end

                # corrected forecast w/ analysis parameters
                soln = pmap((a1,a2,a3,a4,a5,a6)->WK3.odeint(a1,a2,a3,a4,a5,a6),y0,t0,tfcur,sparams,mparams,ke);
                to = [soln[i][1] for i in 1:length(soln)];
                yoa = [soln[i][2] for i in 1:length(soln)];

                # vector of forecast state vectors
                X = [yoa[i][end][1:end] for i in 1:length(soln)];

                # for i in 1:ensemblesize
                #     println("Forecast values, member $i: $(X[i][1:end])")
                # end

                # vector of forecast measurements
                # Y = [yoa[i][end][1:end] for i in 1:length(soln)]; # whole state
                Y = [[yoa[i][end][2],yoa[i][end][3]] for i in 1:length(soln)]; # V, flowrate

                # prior standard deviation
                for i in 1:length(X[1])
                    for j in 1:length(systems)
                        c[j] = X[j][i];
                    end
                    σb[i] = std(c;corrected=true);
                end

                # second forecast mean state, measurement
                xhat = mean(X);
                yhat = mean(Y);
                # println("Normalized ̂x, second forecast: $xhat")
                # println("Normalized ̂y, second forecast: $yhat")

                # output second forecast measurement
                append!(yout,[yhat])
                append!(innov,[yhat-y])

                # # forecast state/meas. cross covariance, meas. covariance
                # Pxy = zeros(length(xhat),length(yhat))
                # Pyy = zeros(length(yhat),length(yhat))
                # for i = 1:length(soln)
                #     Pxy += *((X[i] .- xhat),(Y[i] .- yhat)');
                #     Pyy += *((Y[i] .- yhat),(Y[i] .- yhat)');
                # end
                # Pxy ./= (ensemblesize);
                # Pyy ./= (ensemblesize);
                #
                # # add meas. noise to meas. covariance (allows invertibility)
                # # Pyy += diagm((error.odev[1]/WK3.Qs)^2.*ones(length(yhat)),0);
                # Pyy += diagm([(error.odev[2]/WK3.Vs)^2,(error.odev[3]/WK3.Qs)^2],0);
                # # Pyy += diagm([(error.odev[1]/WK3.Ps)^2,(error.odev[2]/WK3.Vs)^2,(error.odev[3]/WK3.Qs)^2],0);
                #
                # # println("Added noise:")
                # # display(diagm([(error.odev[1]/WK3.Ps)^2,(error.odev[2]/WK3.Qs)^2],0))
                # # println("State cross-covariance:")
                # # display(Pxy)
                # # println("Measurement covariance:")
                # # display(Pyy)
                #
                # # Kalman gain
                # K = Pxy*inv(Pyy);
                # # println("State Kalman gain:")
                # # display(K)
                #
                # # analysis step, NRR tracking
                # for i = 1:length(soln)
                #     X[i][:] += K*(yi[i] .- Y[i]);
                #     r2dot[i] += dot((Y[i].-yi[i]),(Y[i].-yi[i]));
                # end
                # r1dot += sqrt.(dot((yhat.-y),(yhat.-y)));
                #
                # # for i = 1:ensemblesize
                # #     println("Analysis values, member $i: $(X[i])")
                # # end
                #
                # # RTPS multiplicative covariance inflation
                # for i in 1:length(X[1])
                #     for j in 1:length(systems)
                #         c[j] = X[j][i];
                #     end
                #     σa[i] = std(c;corrected=true);
                # end
                # # println("Prior standard deviation: $σb")
                # # println("Posterior standard deviation: $σa")
                # xhat = mean(X);
                # # println("Prior standard deviation: $(σb)")
                # # println("Posterior standard deviation: $(σa)")
                # # println("Normalized difference: $((σb-σa)./σa)")
                # for i = 1:length(soln)
                #     for j = 1:length(X[i])
                #         if σa[j] != 0
                #             X[i][j] .= X[i][j] .+ p.*((σb[j].-σa[j])./σa[j]).*(X[i][j].-xhat[j]);
                #         end
                #     end
                # end
                #
                # # for i = 1:ensemblesize
                # #     println("Analysis values after RTPS, member $i: $(X[i])")
                # # end
                #
                # # output ensemble average state
                # xhat = mean(X);
                append!(xout,[xhat])

                # state variance (for post-processing)
                Pxx = zeros(length(xhat))
                for i = 1:length(soln)
                    Pxx += (X[i] .- xhat).^2;
                end
                Pxx ./= ensemblesize;
                append!(Pxout,[Pxx])

                # state 2-sigma quantiles (for post-processing)
                for i = 1:length(xhat)
                    xq = Float64[];
                    for j = 1:length(soln)
                        push!(xq,X[j][i])
                    end
                    q = quantile(xq,[0.025,0.975]);
                    if i == 1
                        append!(lbx,[ones(1)*q[1]])
                        append!(ubx,[ones(1)*q[2]])
                    else
                        push!(lbx[end],q[1])
                        push!(ubx[end],q[2])
                    end
                end

                # output ensemble average ventricular pressure
                ecur = 0;
                for i = 1:ensemblesize
                    g1 .= (tm./mparams[i].t1).^mparams[i].m1;
                    g2 .= (tm./mparams[i].t2).^mparams[i].m2;
                    h1 .= g1./(1.+g1);
                    h2 .= 1./(1.+g2);
                    ke[i] .= (mparams[i].Emax.-mparams[i].Emin)./maximum(h1.*h2);
                    ei,~ = WK3.elastancefn(tfcur[1]*WK3.ts,mparams[i],ke[i])
                    ecur += ei;
                    systems[i].E[end] = ei;
                    systems[i].Pv[end] = systems[i].E[end].*(X[i][2]*WK3.Vs-systems[i].mparams.V0)/WK3.Ps;
                end
                ecur /= ensemblesize;
                push!(Pvout,ecur.*(xhat[2].*WK3.Vs .- mparams[1].V0)./WK3.Ps)

                # analysis back into ensemble members
                # reshape output to vectors of individual state variables' time series
                for k = 1:ensemblesize # k indexes ensemble members
                    yo[k][end][1:end] = X[k][1:end]; # replace last time step w/ analysis
                    for i = 1:length(yo[k][1]) # i indexes state variables
                        for j = 1:length(yo[k]) # j indexes time steps
                            if i == 1
                                push!(systems[k].Pa,yo[k][j][i])
                            elseif i == 2
                                push!(systems[k].V,yo[k][j][i])
                            elseif i == 3
                                push!(systems[k].Q,yo[k][j][i])
                            end
                            if j == length(yo[k])
                                y0[k][1] = yo[k][j][1];
                                # reset initial volume if ending cardiac cycle
                                if tfcur[1]*WK3.ts == systems[1].mparams.th[end] && nn == length(tf);
                                    y0[k][2] = 1;
                                else
                                    y0[k][2] = yo[k][j][2];
                                end
                                y0[k][3] = yo[k][j][3];
                            end
                        end
                    end
                end

                # reset initial step size
                for i in 1:length(soln)
                    sparams[i].h0 = sparams[i].h0nom;
                end
            elseif nn == length(tf)+1 # make sure we get last set of outputs w/out measurement data
                # reshape output to vectors of individual state variables' time series
                for k = 1:ensemblesize # k indexes ensemble members
                    for i = 1:length(yo[k][1]) # i indexes state variables
                        for j = 1:length(yo[k]) # j indexes time steps
                            if i == 1
                                push!(systems[k].Pa,yo[k][j][i])
                            elseif i == 2
                                push!(systems[k].V,yo[k][j][i])
                            elseif i == 3
                                push!(systems[k].Q,yo[k][j][i])
                            end
                            if j == length(yo[k])
                                y0[k][1] = yo[k][j][1];
                                if tfcur[1]*WK3.ts == systems[1].mparams.th[end];
                                    y0[k][2] = 1;
                                else
                                    y0[k][2] = yo[k][j][2];
                                end
                                y0[k][3] = yo[k][j][3];
                            end
                        end
                    end
                end
            end

            # setup initial time for normal forward time integration
            t0 = [tfcur[1] for i in 1:length(soln)];

            # output measurement time
            if nn <= length(tf)
                push!(tout,t0[1])
            end
        end
    end
    toc()

    # reshape ensemble averages into vectors of individual time series
    for i in 1:length(xout[1]) # i indexes state variables
        xo = Float64[];
        Pxo = Float64[];
        lb = Float64[];
        ub = Float64[];
        for j in 1:length(xout) # j indexes time steps
            push!(xo,xout[j][i])
            push!(Pxo,Pxout[j][i])
            push!(lb,lbx[j][i])
            push!(ub,ubx[j][i])
        end
        append!(xoutv,[xo])
        append!(Pxoutv,[Pxo])
        append!(lbxv,[lb])
        append!(ubxv,[ub])
    end

    for i in 1:length(θout[1]) # i indexes state variables
        θo = Float64[];
        Pθo = Float64[];
        for j in 1:length(θout) # j indexes time steps
            push!(θo,θout[j][i])
            push!(Pθo,Pθout[j][i])
        end
        append!(θoutv,[θo])
        append!(Pθoutv,[Pθo])
    end

    for i in 1:length(yout[1])
        yo = Float64[];
        io = Float64[];
        for j in 1:length(yout)
            push!(yo,yout[j][i])
            push!(io,innov[j][i])
        end
        append!(youtv,[yo])
        append!(ioutv,[io])
    end

    # normalized RMSE ratio (optimal ensemble yields NRR ~ 1)
    r2dot = sqrt.(1/systems[1].t[end]*r2dot);
    R1 = 1/systems[1].t[end]*r1dot;
    R2 = sum(r2dot);
    Ra = R1/R2;
    ERa = sqrt.((ensemblesize)/(2*(ensemblesize)));
    println("Normalized RMSE ratio: $(Ra/ERa)")

    # save for post-processing
    vnames = ["t" "x" "Px" "theta" "Pt" "Pv" "lb" "ub"];
    for i in 1:length(vnames)
        file = MAT.matopen("$(vnames[i])_p.mat","w");
        if vnames[i] == "t"
            write(file,"t",tout)
        elseif vnames[i] == "x"
            write(file,"x",xoutv)
        elseif vnames[i] == "Px"
            write(file,"Px",Pxoutv)
        elseif vnames[i] == "theta"
            write(file,"theta",θoutv)
        elseif vnames[i] == "Pt"
            write(file,"Pt",Pθoutv)
        elseif vnames[i] == "Pv"
            write(file,"Pv",Pvout)
        elseif vnames[i] == "lb"
            write(file,"lb",lbxv)
        elseif vnames[i] == "ub"
            write(file,"ub",ubxv)
        end
        close(file)
    end

    # output
    return systems,tout,xoutv,youtv,ioutv,θoutv,Pθoutv,Pxoutv,Pvout,lbxv,ubxv
end


text:
using OptimizationProblems

probs = filter(name -> name != :OptimizationProblems
               && name != :penalty2
               && name != :penalty3
               , names(OptimizationProblems))

mpb_probs = (MathProgNLPModel(eval(p)(n_max),  name=string(p) )   for p in probs)


text:
using RelevanceStacktrace
using Revise

includet("../src/UnicodePlotsSimple.jl")
using .UnicodePlotsSimple

distribution([0.3,0.2334,0.1,0.8], height=3, values=[0.3,0.2334,0.1,0.8])
vv = abs.(randn(100))
distribution(vv, height=3, values=vv)

rnds = randn(1000).+2
histogram(rnds, width=70, height=5, title="Random num generation histrogram", printstat=true)

xticks = ["a","b","c","h","p"]
yticks = ["a","b","c","h","p"]

data = [randn(5,5),randn(5,5),randn(5,5)]
desnitry = [rand(1:150,5,5),rand(1:1350,5,5),rand(1:6,5,5)]

heatmap(data, title=["first", "sec", "3rd"], xticks=xticks, yticks=yticks, reversey=true)
heatmap(data, format=Val(:ONLYNUMBER), title=["first", "sec", "3rd"], xticks=xticks, yticks=yticks, reversey=true)
densitymap(data, desnitry, title=["first", "sec", "3rd"], xticks=xticks, yticks=yticks, reversey=true)
densitymap(data, desnitry, format=Val(:NUMBER), title=["first", "sec", "3rd"], xticks=xticks, yticks=yticks, reversey=true)

data = [randn(5,5),randn(5,5),randn(5,5)]
desnitry = [rand(1:150,5,5),rand(1:350,5,5),rand(1:6,5,5)]
densitymap(data, desnitry, format=Val(:DISTRIBUTION), title=["first", "sec", "3rd"], xticks=xticks, yticks=yticks, reversey=true)

#%%

using InteractiveUtils
densitymap(data, desnitry,
["first", "s2ec", "323c"], [1,2,2,2,3], [1,2,2,5,3], false, Val(:NUMBER))
@code_warntype densitymap(data, desnitry,
["first", "s2ec", "323c"], [1,2,2,2,3], [1,2,2,5,3], false, Val(:NUMBER))

#%%
heatmap(data,["first", "s2ec", "323c"], [1,2,2,2,3], [1,2,2,5,3], false)
# @code_warntype heatmap(data,["first", "s2ec", "323c"], [1,2,2,2,3], [1,2,2,5,3], false)

#%%
@code_warntype distribution(Float32[0.3,0.2334,0.1,0.8], 3, Float32[0.3,0.2334,0.1,0.8])

#%%
@code_warntype histogram(randn(1000).+2, 80, 5, "Random num generation histrogram", true, false, "","","")

#%%

vv = abs.(randn(120))
distribution(vv, height=3, values=vv)


text:
function recursivecombat!(player1, player2, display = false)

    haswinner = isempty(player1) || isempty(player2)

    round  = 1

    knowdecks = []

    while !haswinner
        if display
            println("Round $round ($(length(player1)) vs $(length(player2)) - cache $(length(knowdecks)))")
            round += 1
        end

        deck = (copy(player1), copy(player2))

        if deck in knowdecks
            return 1
        else
            push!(knowdecks, deck)
        end

        turn!(player1, player2)
        haswinner = isempty(player1) || isempty(player2)
    end

    return isempty(player1) ? 2 : 1
end

function turn!(player1, player2)
    value1 = popfirst!(player1)
    value2 = popfirst!(player2)

    isplayer1win = false

    if value1 > length(player1) || value2 > length(player2)
        isplayer1win = value1 > value2
    else
        subplayer1 = player1[1:value1]
        subplayer2 = player2[1:value2]

        isplayer1win = recursivecombat!(subplayer1, subplayer2) == 1
    end

    isplayer1win ? append!(player1, [value1, value2]) : append!(player2, [value2, value1])
end


text:
using Logging
using MOCNeutronTransport
using Test

# Turn off log info, so we don't clutter up the test output
Logging.disable_logging(Logging.Info)
const Floats = [Float32, Float64, BigFloat]
include("setup/setup_geometry.jl")
include("setup/setup_mesh.jl")
tests = ["quadrature/Quadrature",
         "geometry/Geometry",
         "mesh/Mesh",
         "raytracing/Raytracing"
        ]
for t in tests
  include("$(t).jl")
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# Load packages

using Vann
using PyPlot
using DataFrames
using ExcelReaders
using JLD
using CSV

# Settings

file_exper  = "C:/Work/Studies/vann/experiments.xlsx"
path_inputs = "C:/Work/Studies/vann/data/norway"
path_save   = "C:/Work/Studies/vann/calib_res"

df_exper = readxlsheet(DataFrame, file_exper, "Calib")

# Run one experiment

function run_single_experiment(path_inputs, path_save, df_exper, i_exper)

  # Folder for saving results

  name_experiment = "experiment_" * string(i_exper)

  path_save = joinpath(path_save, name_experiment)

  mkpath(path_save * "/calib_txt")
  mkpath(path_save * "/calib_png")
  mkpath(path_save * "/valid_txt")
  mkpath(path_save * "/valid_png")
  mkpath(path_save * "/param_snow")
  mkpath(path_save * "/param_hydro")
  mkpath(path_save * "/model_data")

  # Get settings

  opt = Dict()

  opt["epot_choice"]  = parse(df_exper[:Epot][i_exper])
  opt["snow_choice"]  = parse(df_exper[:Snow][i_exper])
  opt["hydro_choice"] = parse(df_exper[:Hydro][i_exper])
  opt["force_states"] = parse(df_exper[:force_states][i_exper])

  opt["calib_start"] = df_exper[:calib_start][i_exper]
  opt["calib_stop"]  = df_exper[:calib_stop][i_exper]
  opt["valid_start"] = df_exper[:valid_start][i_exper]
  opt["valid_stop"]  = df_exper[:valid_stop][i_exper]

  opt["tstep"] = df_exper[:tstep][i_exper]

  opt["warmup"] = round(Int64, df_exper[:warmup][i_exper])

  opt["path_save"]   = path_save
  opt["path_inputs"] = path_inputs

  # Run for all stations

  calib_all_stations(opt)

end

# Calibrate all stations for one experiment

function calib_all_stations(opt)

  # Empty dataframes for summary statistics

  df_calib = DataFrame(Station = String[], NSE = Float64[], KGE = Float64[])
  df_valid = DataFrame(Station = String[], NSE = Float64[], KGE = Float64[])

  # Loop over all watersheds

  dir_all = readdir(opt["path_inputs"])

  for dir_cur in dir_all

    # Only run calibration for watersheds without glaciers

    df_meta = CSV.read("$(opt["path_inputs"])/$dir_cur/metadata.txt", delim = ";")

    perc_glacier = df_meta[:perc_glacier]

    if perc_glacier[1].value < 0.5

      print("Calibration for data in $(dir_cur)\n")

      # Run for calibration period

      # Load data

      try
        date, tair, prec, q_obs, frac = load_data("$(opt["path_inputs"])/$dir_cur")
      catch
        info("Unable to read files in directory $(dir_cur)\n")
        continue
      end

      # Crop data

      date, tair, prec, q_obs = crop_data(date, tair, prec, q_obs, opt["calib_start"], opt["calib_stop"])

      # Compute potential evapotranspiration

      epot = eval(Expr(:call, opt["epot_choice"], date))

      # Precipitation correction step

      if sum(~isnan(q_obs)) > (3*365)

        ikeep = ~isnan(q_obs)

        prec_tmp = sum(prec .* repmat(frac, 1, size(prec,2)), 1)

        prec_sum = sum(prec_tmp[ikeep])
        q_sum = sum(q_obs[ikeep])
        epot_sum = sum(epot[ikeep])

        pcorr = (q_sum + 0.5 * epot_sum) / prec_sum

        prec = pcorr * prec

      else

        warn("Not enough runoff data for calibration (see folder $dir_cur)")
        continue

      end

      # Initilize model

      st_snow = eval(Expr(:call, opt["snow_choice"], opt["tstep"], date[1], frac))
      st_hydro = eval(Expr(:call, opt["hydro_choice"], opt["tstep"], date[1]))

      # Run calibration

      param_snow, param_hydro = run_model_calib(st_snow, st_hydro, tair, prec, epot, q_obs;
                                                force_states = opt["force_states"], warmup = opt["warmup"])

      println("Snow model parameters: $param_snow")
      println("Hydro model parameters: $param_hydro")

      # Reinitilize model

      st_snow = eval(Expr(:call, opt["snow_choice"], opt["tstep"], date[1], param_snow, frac))
      st_hydro = eval(Expr(:call, opt["hydro_choice"], opt["tstep"], date[1], param_hydro))

      # Run model with best parameter set

      q_sim, st_snow, st_hydro = run_model(st_snow, st_hydro, tair, prec, epot; return_all = true)

      # Store results in data frame

      q_obs = round(q_obs, 2)
      q_sim = round(q_sim, 2)

      df_res = DataFrame(date = Dates.format(date,"yyyy-mm-dd"), q_sim = q_sim, q_obs = q_obs)

      # Save results to txt file

      file_save = dir_cur[1:end-5]

      writetable(string(opt["path_save"], "/calib_txt/", file_save, "_station.txt"), df_res, quotemark = '"', separator = '\t')

      # Plot results

      ioff()

      file_name = string(opt["path_save"], "/calib_png/", file_save, "_hydro.png")

      plot_sim(st_hydro; q_obs = q_obs, file_name = file_name)

      file_name = string(opt["path_save"], "/calib_png/", file_save, "_snow.png")

      plot_sim(st_snow; file_name = file_name)

      # Compute summary statistics

      station = dir_cur[1:end-5]
      kge_res = kge(q_sim[opt["warmup"]:end], q_obs[opt["warmup"]:end])
      nse_res = nse(q_sim[opt["warmup"]:end], q_obs[opt["warmup"]:end])

      push!(df_calib, [station nse_res kge_res])

      # Run for validation period

      # Reload data

      date, tair, prec, q_obs, frac = load_data("$(opt["path_inputs"])/$dir_cur")

      # Crop data

      date, tair, prec, q_obs = crop_data(date, tair, prec, q_obs, opt["valid_start"], opt["valid_stop"])

      # Compute potential evapotranspiration

      epot = eval(Expr(:call, opt["epot_choice"], date))

      # Reinitilize model

      st_snow = eval(Expr(:call, opt["snow_choice"], opt["tstep"], date[1], param_snow, frac))
      st_hydro = eval(Expr(:call, opt["hydro_choice"], opt["tstep"], date[1], param_hydro))

      # Precipitation correction step

      prec = pcorr * prec

      # Run model with best parameter set

      q_sim, st_snow, st_hydro = run_model(st_snow, st_hydro, tair, prec, epot; return_all = true)

      # Store results in data frame

      q_obs = round(q_obs, 2)
      q_sim = round(q_sim, 2)

      df_res = DataFrame(date = Dates.format(date,"yyyy-mm-dd"), q_sim = q_sim, q_obs = q_obs)

      # Save results to txt file

      file_save = dir_cur[1:end-5]

      writetable(string(opt["path_save"], "/valid_txt/", file_save, "_station.txt"), df_res, quotemark = '"', separator = '\t')

      # Plot results

      ioff()

      file_name = string(opt["path_save"], "/valid_png/", file_save, "_hydro.png")

      plot_sim(st_hydro; q_obs = q_obs, file_name = file_name)

      file_name = string(opt["path_save"], "/valid_png/", file_save, "_snow.png")

      plot_sim(st_snow; file_name = file_name)

      # Save parameter values

      writedlm(opt["path_save"] * "/param_snow/" * file_save * "_param_snow.txt", param_snow)
      writedlm(opt["path_save"] * "/param_hydro/" * file_save * "_param_hydro.txt", param_hydro)

      st_snow = eval(Expr(:call, opt["snow_choice"], opt["tstep"], date[1], param_snow, frac))
      st_hydro = eval(Expr(:call, opt["hydro_choice"], opt["tstep"], date[1], param_hydro))

      jldopen(opt["path_save"] * "/model_data/" * file_save * "_modeldata.jld", "w") do file
        addrequire(file, Vann)
        write(file, "st_snow", st_snow)
        write(file, "st_hydro", st_hydro)
      end

      # Compute summary statistics

      station = dir_cur[1:end-5]
      kge_res = kge(q_sim[opt["warmup"]:end], q_obs[opt["warmup"]:end])
      nse_res = nse(q_sim[opt["warmup"]:end], q_obs[opt["warmup"]:end])

      push!(df_valid, [station nse_res kge_res])

    end

  end

  writetable(string(opt["path_save"], "/summary_calib_period.txt"), df_calib, quotemark = '"', separator = '\t')
  writetable(string(opt["path_save"], "/summary_valid_period.txt"), df_valid, quotemark = '"', separator = '\t')

end

# Run all experiments

for i_exper = 1:size(df_exper, 1)

  run_single_experiment(path_inputs, path_save, df_exper, i_exper)

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
import LightGraphs.SimpleGraphs.SimpleEdge
@testset "SimpleEdge" begin
    e = SimpleEdge(1,2)
    re = SimpleEdge(2,1)

    for s in [0x01, UInt16(1), 1]
        T = typeof(s)
        d = s+one(T)
        p = Pair(s, d)

        ep1 = SimpleEdge(p)
        ep2 = SimpleEdge{UInt8}(p)
        ep3 = SimpleEdge{Int16}(p)

        t1 = (s, d)
        t2 = (s, d, "foo")

        @test src(ep1) == src(ep2) == src(ep3) == s
        @test dst(ep1) == dst(ep2) == dst(ep3) == s + one(T)

        @test eltype(p) == typeof(s)
        @test SimpleEdge(p) == e
        @test SimpleEdge(t1) == SimpleEdge(t2) == e
        @test SimpleEdge(t1) == SimpleEdge{UInt8}(t1) == SimpleEdge{Int16}(t1)
        @test SimpleEdge{Int64}(ep1) == e

        @test Pair(e) == p
        @test Tuple(e) == t1
        @test reverse(ep1) == re
        @test sprint(show, ep1) == "Edge 1 => 2"
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
#=
how to compute mic with two variable
x = linespace(0,1,10) |> Array
y = sin(x)
f = MINE(method=:mic, alpha=0.6, c=15)
# formula = MINE(method=[:mic,:mas], alpha=0.6, c=15)
s = f(x, y)
=#

function MINE(fn::Symbol=:mic;
              alpha::Float64=0.6, c::Float64=15.0)

    @assert 0 <= alpha <= 1
    @assert c > 0
    return score = problem_score(alpha, c)
    @show score
#    return (x,y)-> (score(x,y) |> eval(CAPI, fn))
    return (x,y)-> (score(x,y) |> CAPI.mic)
end

function MINE(fns::Array{Symbol,1}=[:mic];
              alpha::Float64=0.6, c::Float64=15.0)

    @assert 0 <= alpha <= 1
    @assert c > 0

    (x,y) -> map(fn->MINE(fn, alpha=alpha, c=c)(x, y), fns)
end

function problem_score(alpha, c)
    param  = Parameter(alpha, c)
    pparam = convert(Ptr{Parameter}, pointer_from_objref(param))

    function problem_score(x::Array{Float64,1}, y::Array{Float64,1})
        prob  = Problem(x, y)
        pprob = convert(Ptr{Problem}, pointer_from_objref(prob))
        CAPI.compute_score(pprob, pparam)
    end
end

function unsafe_load_score(pscore::Ptr{Score})
    pscore
end

function main()
    f = MINE(:mic)
    f([1,2.0], [2.0, 4])
end


text:
function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.EqualTo)
    set_lower_bound(info, set.value)
    set_upper_bound(info, set.value)
end

function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.GreaterThan)
    set_lower_bound(info, set.lower)
end

function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.LessThan)
    set_upper_bound(info, set.upper)
end

function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.Interval)
    set_lower_bound(info, set.lower)
    set_upper_bound(info, set.upper)
end

MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{<:Bounds}) = true

function MOIU.load_constraint(model::Optimizer, ci::CI, f::SV, set::Bounds)
    variable_info = find_variable_info(model, f.variable)
    set_bounds(variable_info, set)
    return
end

MOI.supports_constraint(::Optimizer, ::Type{<:Union{SAF, SQF}}, ::Type{<:Bounds}) = true

function MOIU.load_constraint(model::Optimizer, ci::CI, f::Union{SAF, SQF}, set::Bounds)
    check_variable_indices(model, f)
    constraint_info = model.inner.constraint_info[ci.value]
    constraint_info.expression = to_expr(f)
    set_bounds(constraint_info, set)
    return
end

MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{MOI.ZeroOne}) = true
MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{MOI.Integer}) = true

function MOIU.load_constraint(model::Optimizer, ci::CI, f::SV, set::Union{MOI.ZeroOne, MOI.Integer})
    vi = f.variable
    check_variable_indices(model, vi)
    variable_info = model.inner.variable_info[vi.value]
    if set isa MOI.ZeroOne
        variable_info.category = :Bin
    elseif set isa MOI.Integer
        variable_info.category = :Int
    else
        error()
    end
    return
end

# MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

# function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_block_data::MOI.NLPBlockData)
#     model.nlp_block_data = nlp_block_data
#     return
# end


text:
x=1:10
y=rand(10)
p1 = plot(x,y) # Make a line plot
p2 = scatter(x,y) # Make a scatter plot
p3 = plot(x,y,xlabel="This one is labelled",lw=3,title="Subtitle")
p4 = histogram(x,y) # Four histograms each with 10 points? Why not!
plot(p1,p2,p3,p4,layout=(2,2),legend=false)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module julia_sets
using PyPlot
export gen_jset,show_jset

function gen_jset{T<:Real,U<:Real}(R::Function,x::Array{T,1},y::Array{T,1},n_iter::Int64,escape_tol::U)

  A = zeros(length(x),length(y));
  for i=1:length(x)
    for j=1:length(y)
      z = Complex(x[i],y[j])
      for k = 1:n_iter
        z = R(z)
        if abs(z) > escape_tol
          A[i,j] = k
          break;
        end
      end
      if A[i,j]==0
        A[i,j] = escape_tol+1;
        end
    end
  end

  return A
end

function show_jset{T<:Real}(A::Array{T,2})
 matshow(A)
end

end


text:
import Base.show

export loo_compare, ModelComparison

"""
    ModelComparison

A struct containing the results of model comparison.

# Fields

- `pointwise::KeyedArray`: A `KeyedArray` of pointwise estimates. See [`PsisLoo`]@ref.
  - `estimates::KeyedArray`: A table containing the results of model comparison, with the
    following columns --
      + `cv_elpd`: The difference in total leave-one-out cross validation scores
        between models.
      + `cv_avg`: The difference in average LOO-CV scorees between models.
      + `weight`: A set of Akaike-like weights assigned to each model, which can be used in
        pseudo-Bayesian model averaging.
  - `std_err::NamedTuple`: A named tuple containing the standard error of `cv_elpd`. Note
    that these estimators (incorrectly) assume all folds are independent, despite their
    substantial overlap, which creates a downward biased estimator. LOO-CV differences are
    *not* asymptotically normal, so these standard errors cannot be used to
    calculate a confidence interval.
  - `gmpd::NamedTuple`: The geometric mean of the predictive distribution. It equals the
    geometric mean of the probability assigned to each data point by the model, that is,
    `exp(cv_avg)`. This measure is only meaningful for classifiers (variables with discrete
    outcomes). We can think of it as measuring how often the model was right: A model that
    always predicts incorrectly will have a GMPD of 0, while a model that always predicts
    correctly will have a GMPD of 1. However, the GMPD gives a model "Partial points"
    between 0 and 1 whenever the model assigns a probability other than 0 or 1 to the
    outcome that actually happened.

See also: [`PsisLoo`](@ref)
"""
struct ModelComparison{RealType<:Real, N}
    pointwise::KeyedArray{RealType, 3, <:NamedDimsArray, <:Any}
    estimates::KeyedArray{RealType, 2, <:NamedDimsArray, <:Any}
    std_err::NamedTuple{<:Any, Tuple{Vararg{RealType, N}}}
    gmpd::NamedTuple{<:Any, Tuple{Vararg{RealType, N}}}
end

"""
    function loo_compare(
        cv_results...;
        sort_models::Bool=true,
        best_to_worst::Bool=true,
        [, model_names::Tuple{Symbol}]
    ) -> ModelComparison

Construct a model comparison table from several [`PsisLoo`](@ref) objects.

# Arguments
  - `cv_results`: One or more [`PsisLoo`](@ref) objects to be compared. Alternatively,
    a tuple or named tuple of `PsisLoo` objects can be passed. If a named tuple is passed,
    these names will be used to label each model.
  - `model_names`: A vector or tuple of strings or symbols used to identify models. If
    none, models are numbered using the order of the arguments.
  - `sort_models`: Sort models by total score.
  - `high_to_low`: Sort models from best to worst score. If `false`, reverse the order.

See also: [`ModelComparison`](@ref), [`PsisLoo`](@ref), [`psis_loo`](@ref)
"""
function loo_compare(
    cv_results::AbstractVector{<:PsisLoo};
    model_names::S=[Symbol("model_$i") for i in 1:length(cv_results)],
    sort_models::Bool=true,
    high_to_low::Bool=true,
) where S <: Base.AbstractVecOrTuple{<:Union{AbstractString, Symbol}}

    model_names = [Symbol(model_names[i]) for i in 1:length(model_names)]  # array version
    n_models, data_size = _get_dims(cv_results, model_names)

    @views if sort_models
        sorting_by = x -> x.estimates(:cv_elpd, :total)
        order = sortperm(cv_results; by=sorting_by, rev=high_to_low)
        cv_results = cv_results[order]
        model_names = model_names[order]
    end

    # Extract relevant values from PsisLoo objects
    @views begin
        cv_elpd = [cv_results[i].estimates(:cv_elpd, :total) for i in 1:n_models]
        pointwise = [cv_results[i].pointwise for i in 1:n_models]
    end

    # Construct pointwise differences
    pointwise_diffs = cat(pointwise...; dims=3)
    pointwise_diffs = KeyedArray(
        pointwise_diffs;
        data=1:data_size,
        statistic=[:cv_elpd, :naive_lpd, :p_eff, :mcse, :pareto_k],
        model=model_names,
    )
    # Subtract the effective number of params and elpd ests; leave mcse+pareto_k the same
    @views base_case = pointwise_diffs[:, 1:3, 1]
    @views @. pointwise_diffs[:, 1:3, Not(1)] = pointwise_diffs[:, 1:3, Not(1)] - base_case
    pointwise_diffs[:, 1:3, 1] .= 0

    name_tuple = ntuple(i -> model_names[i], n_models)  # convert to tuple

    variance = var(pointwise_diffs(data=:, statistic=:cv_elpd, model=:); dims=:data)
    se_total = @. sqrt(data_size * variance)
    se_total = ntuple(i -> se_total[i], n_models)
    se_total = NamedTuple{name_tuple}(se_total)

    log_norm = logsumexp(cv_elpd)
    weights = @. exp(cv_elpd - log_norm)

    gmpd = @. exp(cv_elpd / data_size)
    gmpd = NamedTuple{name_tuple}(gmpd)

    @. cv_elpd = cv_elpd - cv_elpd[1]
    avg_elpd = cv_elpd ./ data_size
    total_diffs = KeyedArray(
        hcat(cv_elpd, avg_elpd, weights);
        model=model_names,
        statistic=[:cv_elpd, :cv_avg, :weight],
    )

    return ModelComparison(pointwise_diffs, total_diffs, se_total, gmpd)

end

function loo_compare(cv_results::NamedTuple; kwargs...)
    @nospecialize cv_results
    return loo_compare(values(cv_results); model_names=keys(cv_results), kwargs...)
end

function loo_compare(cv_results::Tuple; kwargs...)
    @nospecialize cv_results
    @views cv_results = [cv_results[i] for i in 1:length(cv_results)]
    return loo_compare(cv_results; kwargs...)
end

function loo_compare(cv_results::PsisLoo...; kwargs...)
    @nospecialize cv_results
    @views cv_results = [cv_results[i] for i in 1:length(cv_results)]
    return loo_compare(cv_results; kwargs...)
end

function Base.show(io::IO, ::MIME"text/plain", model_comparison::ModelComparison)
    estimates = model_comparison.estimates
    return pretty_table(
        estimates;
        compact_printing=false,
        header=estimates.statistic,
        row_names=estimates.model,
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end

function _get_dims(cv_results::AbstractVector{<:PsisLoo}, model_names::Vector{Symbol})

    n_models = length(cv_results)
    if n_models != length(model_names)
        throw(ArgumentError("Number of model names does not match size of `cv_results`."))
    end

    data_sizes = [cv_results[i].psis_object.data_size for i in 1:n_models]
    data_size = data_sizes[1]
    if any(data_sizes .≠ data_size)
        throw(
            ArgumentError(
                "Number of data points differs across `PsisLoo` objects, implying models " *
                "were not trained on the same sample. LOO-CV cannot be used."
            )
        )
    end

    return n_models, data_size
end

max_stars_repo_name: ParadaCarleton/JuLOOa

max_issues_repo_name: ParadaCarleton/JuLOOa

max_forks_repo_name: ParadaCarleton/JuLOOa


text:
"""
    ncvar = NCDatasets.ancillaryvariables(ncv::NCDatasets.CFVariable,modifier)

Return the first ancillary variables from the NetCDF variable `ncv` with the
standard name modifier `modifier`. It can be used for example to access
related variable like status flags.
"""
function ancillaryvariables(ncv::NCDatasets.CFVariable,modifier)
    ds = Dataset(ncv.var.ncid,ncv.var.isdefmode)
    varname = name(ncv)

    if !haskey(ncv.attrib,"ancillary_variables")
        return nothing
    end

    ancillary_variables = split(ncv.attrib["ancillary_variables"])

    for j = 1:length(ancillary_variables)
        ncv_ancillary = ds[ancillary_variables[j]]
        if occursin(modifier,ncv_ancillary.attrib["standard_name"])
            @debug ancillary_variables[j]
            return ncv_ancillary
        end
    end

    # nothing found
    return nothing
end

"""
    data = NCDatasets.filter(ncv, indices...; accepted_status_flags = nothing)

Load and filter observations by replacing all variables without an acepted status
flag to `missing`. It is used the attribute `ancillary_variables` to identify
the status flag.

```
# da["data"] is 2D matrix
good_data = NCDatasets.filter(ds["data"],:,:, accepted_status_flags = ["good_data","probably_good_data"])
```

"""
function filter(ncv, indices...; accepted_status_flags = nothing)
#function filter_(ncv, indices...)
#    accepted_status_flags = ("good_value", "probably_good_value")
    data = ncv[indices...];

    if (accepted_status_flags != nothing)
        ncv_ancillary = ancillaryvariables(ncv,"status_flag");
        if ncv_ancillary == nothing
            error("no variable with the attribute status_flag as standard_name among $(ancillary_variables) found")
        end

        flag_values = ncv_ancillary.attrib["flag_values"]
        flag_meanings = ncv_ancillary.attrib["flag_meanings"]::String
        if typeof(flag_meanings) <: AbstractString
            flag_meanings = split(flag_meanings)
        end

        accepted_status_flag_values = zeros(eltype(flag_values),length(accepted_status_flags))
        for i = 1:length(accepted_status_flags)
            tmp = findfirst(accepted_status_flags[i] .== flag_meanings)

            if tmp == nothing
                error("cannot recognise flag $(accepted_status_flags[i])")
            end
            accepted_status_flag_values[i] = flag_values[tmp]
        end
        #@debug accepted_status_flag_values

        dataflag = ncv_ancillary.var[indices...];
        for i in eachindex(data)
            good = false;
            for accepted_status_flag_value in accepted_status_flag_values
                good = good || (dataflag[i] .== accepted_status_flag_value)
            end
            if !good
                #@show i,dataflag[i]
                data[i] = missing
            end
        end
    end

    return data
end


text:
module Plots
    using Makie
    using JSServe
    using Colors
    using Meshing
    using MeshIO
    using Markdown
    using GeometryBasics
    using DataFrames
    using FileIO
    import ..AllenNeuropixels as AN

    function htmlexport(fig, file::String="plot.html")
        open(file, "w") do io
                println(io, """
                <html>
                    <head>
                    </head>
                    <body>
                """)
                show(io, MIME"text/html"(), Page(exportable=true, offline=true))
                show(io, MIME"text/html"(), fig)
                println(io, """
                    </body>
                </html>
                """)
        end
    end
    exporthtml = htmlexport

    include("./ReferenceAtlas.jl")
    include("./LFP.jl")
end


text:
using StatsBase: mode, mean

abstract type TreeNode end

mutable struct ClassificationNode <: TreeNode
    depth::Int
    isleaf::Bool
    regionid::Union{Nothing, Int}
    prediction::Union{Nothing, Int}
    featureindex::Union{Nothing, Int}
    thresh::Union{Nothing, Float64}
    leftnode::Union{Nothing, TreeNode}
    rightnode::Union{Nothing, TreeNode}

    function ClassificationNode(depth)
        new(depth, false, nothing, nothing, nothing, nothing, nothing, nothing)
    end

end

struct TreeStoppingRule
    maxdepth::Int
    minsamplesplit::Int
end

"""
Tree node for regression task.
"""
mutable struct RegressionNode <: TreeNode
    depth::Int
    isleaf::Bool
    regionid::Union{Nothing, Int}
    prediction::Union{Nothing, Float64}
    featureindex::Union{Nothing, Int}
    thresh::Union{Nothing, Float64}
    leftnode::Union{Nothing, TreeNode}
    rightnode::Union{Nothing, TreeNode}

    function RegressionNode(depth)
        new(depth, false, nothing, nothing, nothing, nothing, nothing, nothing)
    end

end

abstract type DecisionTree end

struct ClassificationTree <: DecisionTree
    rootnode::TreeNode
    lossfunc::Function
    stoppingrule::TreeStoppingRule
    numregions::Int

    function ClassificationTree(X, y, lossfunc, stoppingrule)
        rootnode, lastregionid = fittree!(ClassificationNode(1), X, y, lossfunc, stoppingrule)
        new(rootnode, lossfunc, stoppingrule, lastregionid)
    end
end

struct RegressionTree <: DecisionTree
    rootnode::TreeNode
    lossfunc::Function
    stoppingrule::TreeStoppingRule
    numregions::Int

    function RegressionTree(X, y, lossfunc, stoppingrule)
        rootnode, lastregionid = fittree!(RegressionNode(1), X, y, lossfunc, stoppingrule)
        new(rootnode, lossfunc, stoppingrule, lastregionid)
    end
end

function RegressionTree(X, y, stoppingrule)
    return RegressionTree(X, y, squaredloss, stoppingrule)
end

function worthsplitting(y, leftindices, currentdepth, stoppingrule)
    rightindices = .!leftindices
    if (currentdepth > stoppingrule.maxdepth ||
        !checkminsamplesplit(sum(leftindices), sum(rightindices),
                             stoppingrule.minsamplesplit))
        return false
    else
        return true
    end
end

function checkminsamplesplit(leftnum, rightnum, minsamplesplit)
    return leftnum >= minsamplesplit && rightnum >= minsamplesplit
end

function findbestsplit(X, y, lossfunc, minsamplesplit)
    # no split case, loss value is evluated on all observations
    bestfeature = 1
    bestthresh = -Inf
    lossval = lossfunc(y)

    for (featureindex, featurearray) in enumerate(eachcol(X))
        sortedindices = sortperm(featurearray)
        sortedX = X[sortedindices, :]
        sortedy = y[sortedindices]
        for i in 1:(length(sortedy) - 1)
            templossval = lossfunc(sortedy[1:i], sortedy[i+1:end])
            if (templossval < lossval && checkminsamplesplit(i, length(y)-i, minsamplesplit))
                lossval = templossval
                bestfeature = featureindex
                bestthresh = (sortedX[i, featureindex] + sortedX[i+1, featureindex])/2
            end
        end
    end

    return bestfeature, bestthresh
end

function split(X, y, lossfunc, minsamplesplit)
    featureindex, thresh = findbestsplit(X, y, lossfunc, minsamplesplit)
    leftindices = X[:, featureindex] .< thresh
    return featureindex, thresh, leftindices
end

function _buildbranch(node::N, X, y, indices, depth) where N<:TreeNode
    selectedX = X[indices, :]
    selectedy = y[indices]
    node = N(depth)
    return selectedX, selectedy, node
end

function predict!(node::ClassificationNode, y)
    node.prediction = mode(y)
    nothing
end

function predict!(node::RegressionNode, y)
    node.prediction = mean(y)
    nothing
end

function _createyhat(model::RegressionTree, X)
    zeros(Float64, size(X)[1])
end

function _createyhat(model::ClassificationTree, X)
    zeros(Int, size(X)[1])
end

function predict(model::DecisionTree, X)
    yhat = _createyhat(model, X)
    for (i, x) in enumerate(eachrow(X))
        node = model.rootnode
        while !node.isleaf
            if x[node.featureindex] < node.thresh
                node = node.leftnode
            else
                node = node.rightnode
            end
        end
        yhat[i] = node.prediction
    end
    return yhat
end

function fittree!(node::TreeNode, X, y, lossfunc, stoppingrule)
    function fit!(node::TreeNode, X, y)
        featureindex, thresh, leftindices = split(X, y, lossfunc, stoppingrule.minsamplesplit)
        nextdepth = node.depth + 1

        if !worthsplitting(y, leftindices, nextdepth, stoppingrule)
            predict!(node, y)
            regionid += 1
            node.regionid = regionid
            node.isleaf = true
            return node
        else
            node.featureindex = featureindex
            node.thresh = thresh
            rightindices = .!leftindices

            leftX, lefty, leftnode = _buildbranch(node, X, y, leftindices, nextdepth)
            node.leftnode = fit!(leftnode, leftX, lefty)

            rightX, righty, rightnode = _buildbranch(node, X, y, rightindices, nextdepth)
            node.rightnode = fit!(rightnode, rightX, righty)
            return node
        end
    end

    regionid = 0
    fittednode = fit!(node::TreeNode, X, y)
    return fittednode, regionid
end

# Loss functions
function squaredloss(ys...)
    totalloss = 0.0
    for y in ys
        yhat = mean(y)
        totalloss += sum((y.-yhat).^2)
    end
    return totalloss
end

function _giniindex(y)
    num_ones = float(sum(y.==1))
    total_num = length(y)
    p = num_ones / total_num
    return 2*p*(1-p)
end

function giniindex(ys...)
    totalval = 0.0
    totalnum = sum([length(y) for y in ys])
    for y in ys
        n = float(length(y))
        totalval += n / totalnum * _giniindex(y)
    end
    return totalval
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
"""
Valuation of a product for given market data
"""
function price end

# product generic
function price(product:: AbstractProduct, container:: AbstractMarketDataContainer, valuation_date:: Date)
    market_data = get_specific_market_data(product, container)
    return price(product, market_data, valuation_date)
end

# product specific because of different market data requirements
function price(product:: FXForward, market_data:: SpecificFXForwardMarketData, valuation_date:: Date)
    fx_forward_curves = market_data.fx_forward_curves
    discount_curves = market_data.discount_curves
    return price(product, fx_forward_curves, discount_curves, valuation_date; flip_fx_direction=market_data.flip_fx_direction)
end

function price(product:: FXForward, fx_forward_curves:: FXForwardTimeSeries, discount_curves:: DiscountTimeSeries, valuation_date:: Date; flip_fx_direction=false)
    fx_forward_curve = fx_forward_curves[valuation_date]
    discount_curve = discount_curves[valuation_date]
    return price(product, fx_forward_curve, discount_curve; flip_fx_direction=flip_fx_direction)
end

function price(product:: FXForward, fx_forward_curve:: Curve, discount_curve:: Curve; flip_fx_direction=false)
    payoff = product.notional_base_currency * product.forward_fx_rate
    interpolated_fx_forward_rate = interpolate(product.tenor, fx_forward_curve)
    par_payoff = flip_fx_direction ? product.notional_base_currency / interpolated_fx_forward_rate : product.notional_base_currency * interpolated_fx_forward_rate
    pv = (payoff - par_payoff) * interpolate(product.tenor, discount_curve)
    return pv
end

"""
Prices a complete time series of a product
"""
function price_time_series(product:: AbstractProduct, container:: AbstractMarketDataContainer)
    market_data = get_specific_market_data(product, container)
    return price_time_series(product, market_data)
end

function price_time_series(product:: FXForward, market_data:: SpecificFXForwardMarketData)
    fx_forward_curves = values(market_data.fx_forward_curves.curves)
    discount_curves = values(market_data.discount_curves.curves)
    pvs = [price(product, fx_curve, discount_curve; flip_fx_direction=market_data.flip_fx_direction) for (fx_curve, discount_curve) in zip(fx_forward_curves, discount_curves)]
    return pvs
end

max_issues_repo_name: lungben/MarketRiskPrototype

max_forks_repo_name: lungben/MarketRiskPrototype


text:
using Thebes, Luxor, Random, Rotations

function main()
    Drawing(800, 800, "test1.svg")
    origin()
    background("grey95")
    setopacity(0.7)
    setlinejoin("bevel")
    eyepoint(Point3D(500, 500, 100))
    uppoint(Point3D(0, 0, 20))
    centerpoint(Point3D(0, 0, 10))
    tiles = Tiler(800, 800, 6, 6)
    fontsize(20)
    for (pos, n) in tiles
        @layer begin
            r = RotXYZ(0, 0, 0)
            translate(pos)
            axes3D(150)
            object = make(Cube)
            scaleby!(object, 50, 50, 50)
            moveby!(object, 0, 0, 0)
            rotateby!(object, 0, n/5, 0)
            sortfaces!(object)
            pin(object)
            sethue("black")
            text(string(n), O  + (0, 30))
            fontsize(50)
            r *= RotXYZ(0, n/5, 0)
            text3D("$n", halign=:center, object.vertices[1], rotation = r)
        end
    end
    finish()
    preview()
end

main()


text:
using RobustFactorizations, ProximalOperators
using Test
using Random, Statistics, LinearAlgebra

@testset "RobustFactorizations.jl" begin

    @testset "Robust PCA" begin
        @info "Testing Robust PCA"

        @testset "Nonneg" begin
            @info "Testing Nonneg"

            D = [0.462911    0.365901  0.00204357    0.692873    0.935861;
            0.0446199    0.108606   0.0664309   0.0736707    0.264429;
            0.320581    0.287788    0.073133    0.188872    0.526404;
            0.356266    0.197536 0.000718338    0.513795    0.370094;
            0.677814    0.011651    0.818047   0.0457694    0.471477]

            A = [0.462911   0.365901  0.00204356   0.345428   0.623104;
            0.0446199  0.108606  0.0429271    0.0736707  0.183814;
            0.320581   0.203777  0.073133     0.188872   0.472217;
            0.30725    0.197536  0.000717701  0.201626   0.370094;
            0.234245   0.011651  0.103622     0.0457694  0.279032]

            E = [0.0        0.0        0.0        0.347445  0.312757 ;
            0.0        0.0        0.0235038  0.0       0.0806151;
            0.0        0.0840109  0.0        0.0       0.0541868;
            0.0490157  0.0        6.5061e-7  0.312169  0.0      ;
            0.443569   0.0        0.714425   0.0       0.192445]

            res = rpca(D, nonnegS=true, nonnegL=true, verbose=false)

            @test res.L ≈ A atol=1e-6
            @test res.S ≈ E atol=1e-6
            @test norm(D - (res.L + res.S))/norm(D) < 1e-6

            res = rpca(D, nonnegS=false, nonnegL=false, verbose=false)
            @test norm(D - (res.L + res.S))/norm(D) < 1e-6
        end

        @testset "Complex" begin
            @info "Testing Complex"
            u = randn(ComplexF64, 100)
            v = randn(ComplexF64, 20)
            E = randn(ComplexF64,100,20) .* 10 .* (rand.() .< 0.01)
            A = u*v'
            D = A .+ E
            res = rpca(D)
            @test sum(abs2, res.S-E)/sum(abs2, E) < 1e-5
            @test sum(abs2, res.L-A)/sum(abs2, A) < 1e-5
        end

        @testset "ElasticNet" begin
            @info "Testing ElasticNet"

            # Dense and sparse noise
            D = lowrank(100,10,3)
            E = 0.1randn(100,10) + sparserandn(100,10)
            Dn = D + E
            λ = 1/sqrt(100)
            res1 = rpca(Dn, verbose=false)
            res2 = rpca(Dn, verbose=false, proxD=SqrNormL2(0.5))
            @test norm(Dn - (res1.L + res1.S + res1.D))/norm(D) < 1e-7
            @show norm(D - res1.L) > norm(D - res2.L) # This test was not robust

            N = 1000; r=4;s=1e-1
            for N in [100,1000,10_000], r in [2,4,6]
                D = lowrank(N,10,r)
                for s in [1e-3, 1e-2, 1e-1]
                    @show N,r,s
                    λ = 1/sqrt(N)

                    # Only sparse noise
                    E = 10sparserandn(N,10)
                    Dn = D + E
                    res1 = rpca(Dn, verbose=false)
                    res2 = rpca(Dn, verbose=false, proxD=SqrNormL2(λ/(s)))
                    @test norm(D - res1.L) < norm(D - res2.L)

                    # Dense and large sparse noise
                    E .+= s*randn(N,10)
                    Dn = D + E
                    res1 = rpca(Dn, verbose=false)
                    res2 = rpca(Dn, verbose=false, proxD=SqrNormL2(λ/(s)))
                    @show norm(D - res1.L)
                    @show norm(D - res2.L)

                    # @test norm(D - Â1) > norm(D - Â2)
                end
            end
        end
    end

    @testset "Robust PCA FISTA" begin
        @info "Testing Robust PCA FISTA"

        @testset "Nonneg" begin
            @info "Testing Nonneg"

            D = [0.462911    0.365901  0.00204357    0.692873    0.935861;
            0.0446199    0.108606   0.0664309   0.0736707    0.264429;
            0.320581    0.287788    0.073133    0.188872    0.526404;
            0.356266    0.197536 0.000718338    0.513795    0.370094;
            0.677814    0.011651    0.818047   0.0457694    0.471477]

            A = [0.462911   0.365901  0.00204356   0.345428   0.623104;
            0.0446199  0.108606  0.0429271    0.0736707  0.183814;
            0.320581   0.203777  0.073133     0.188872   0.472217;
            0.30725    0.197536  0.000717701  0.201626   0.370094;
            0.234245   0.011651  0.103622     0.0457694  0.279032]

            E = [0.0        0.0        0.0        0.347445  0.312757 ;
            0.0        0.0        0.0235038  0.0       0.0806151;
            0.0        0.0840109  0.0        0.0       0.0541868;
            0.0490157  0.0        6.5061e-7  0.312169  0.0      ;
            0.443569   0.0        0.714425   0.0       0.192445]

            res = rpca_fista(D, nonnegS=true, nonnegL=true, verbose=false)

            @test norm(res.L-A)/norm(A) < 0.9
            @test norm(res.S-E)/norm(E) < 0.9
            @test norm(D - (res.L + res.S))/norm(D) < 0.9

            res = rpca_fista(D, nonnegS=false, nonnegL=false, verbose=false, iters=50000, λ=0.2)
            @test norm(res.L-A)/norm(A) < 0.9
            @test norm(res.S-E)/norm(E) < 0.9
        end

        @testset "Complex" begin
            @info "Testing Complex"
            u = randn(ComplexF64, 100)
            v = randn(ComplexF64, 20)
            E = randn(ComplexF64,100,20) .* 10 .* (rand.() .< 0.01)
            A = u*v'
            D = A .+ E
            res = rpca_fista(D)
            @test sum(abs2, res.S-E)/sum(abs2, E) < 1e-3
            @test sum(abs2, res.L-A)/sum(abs2, A) < 1e-3
        end
    end

    @testset "Robust PCA ADMM" begin
        @info "Testing Robust PCA ADMM"

        @testset "Nonneg" begin
            @info "Testing Nonneg"

            D = [0.462911    0.365901  0.00204357    0.692873    0.935861;
            0.0446199    0.108606   0.0664309   0.0736707    0.264429;
            0.320581    0.287788    0.073133    0.188872    0.526404;
            0.356266    0.197536 0.000718338    0.513795    0.370094;
            0.677814    0.011651    0.818047   0.0457694    0.471477]

            A = [0.462911   0.365901  0.00204356   0.345428   0.623104;
            0.0446199  0.108606  0.0429271    0.0736707  0.183814;
            0.320581   0.203777  0.073133     0.188872   0.472217;
            0.30725    0.197536  0.000717701  0.201626   0.370094;
            0.234245   0.011651  0.103622     0.0457694  0.279032]

            E = [0.0        0.0        0.0        0.347445  0.312757 ;
            0.0        0.0        0.0235038  0.0       0.0806151;
            0.0        0.0840109  0.0        0.0       0.0541868;
            0.0490157  0.0        6.5061e-7  0.312169  0.0      ;
            0.443569   0.0        0.714425   0.0       0.192445]

            res = rpca_admm(D, nonnegS=true, nonnegL=true, verbose=true)

            @test norm(res.L-A)/norm(A) < 0.2
            @test norm(res.S-E)/norm(E) < 0.2
            @test norm(D - (res.L + res.S))/norm(D) < 1e-8

        end

        @testset "Complex" begin
            @info "Testing Complex"
            u = randn(ComplexF64, 100)
            v = randn(ComplexF64, 20)
            E = randn(ComplexF64,100,20) .* 10 .* (rand.() .< 0.01)
            A = u*v'
            D = A .+ E
            res = rpca_admm(D,ρ=50, tol=1e-5, iters=5000)
            @test sum(abs2, res.S-E)/sum(abs2, E) < 0.01
            @test sum(abs2, res.L-A)/sum(abs2, A) < 0.01
        end
    end

end


text:
import ForwardDiff
import LegibleLambdas
import LegibleLambdas: LegibleLambda
import ACE: read_dict, write_dict

export λ, lambda

struct Lambda{TL}
   ll::TL
   exstr::String
end

Base.show(io::IO, t::Lambda) = print(io, "λ($(t.exstr))")

function λ(str::String)
   ex = Meta.parse(str)
   ll = LegibleLambda(ex, eval(ex))
   return Lambda(ll, str)
end

lambda(str::String) = λ(str)

(t::Lambda)(x) = evaluate(t, x)

evaluate(t::Lambda, x) = t.ll.λ(x)

evaluate_d(t::Lambda, x::SVector) =  ForwardDiff.gradient(t.ll.λ, x)

ACE.evaluate_ed(t::Lambda, x) =  evaluate(t, x), evaluate_d(t, x)

evaluate_dd(t::Lambda, x::SVector) = ForwardDiff.hessian(t.ll.λ, x)

evaluate_d(t::Lambda, x::Real) =     ForwardDiff.derivative(t.ll.λ, x)

evaluate_dd(t::Lambda, x::Real) =     ForwardDiff.derivative(y -> evaluate_d(t, y), x)

function frule_evaluate(t::Lambda, x::Real, dx::SVector)
      f = evaluate(t, x)
      df = evaluate_d(t, x)
      return f, df * dx
end

write_dict(t::Lambda)  = Dict(
         "__id__" => "ACE_Lambda",
         "exstr" => t.exstr
      )

read_dict(::Val{:ACE_Lambda}, D::Dict) = λ(D["exstr"])

import Base: ==

==(F1::Lambda, F2::Lambda) = (F1.exstr == F2.exstr)


text:
using Oceananigans.Grids: Center, Face

"""
Notes:
This file defines grid lengths, areas, and volumes for staggered structured grids.
Each "reference cell" is associated with an index i, j, k.
The "location" of each reference cell is roughly the geometric centroid of the reference cell.
On the staggered grid, there are 7 cells additional to the "reference cell"
that are staggered with respect to the reference cell in x, y, and/or z.

The staggering is denoted by the locations "Center" and "Face":
    - "Center" is shared with the reference cell;
    - "Face" lies between reference cell centers, roughly at the interface between
      reference cells.

The three-dimensional location of an object is defined by a 3-tuple of locations, and
denoted by a triplet of superscripts. For example, an object `φ` whose cell is located at
(Center, Center, Face) is denoted `φᶜᶜᶠ`. `ᶜᶜᶠ` is Centered in `x`, `Centered` in `y`, and on
reference cell interfaces in `z` (this is where the vertical velocity is located, for example).
The super script `ᵃ` denotes "any" location.

The operators in this file fall into three categories:

1. Operators needed for an algorithm valid on rectilinear grids with
   at most a stretched vertical dimension and regular horizontal dimensions.
2. Operators needed for an algorithm on a grid that is curvilinear in the horizontal
   at rectilinear (possibly stretched) in the vertical.
"""

#####
##### Grid lengths for horizontally-regular algorithms
#####

@inline Δx(i, j, k, grid::ARG) = grid.Δx
@inline Δy(i, j, k, grid::ARG) = grid.Δy

@inline ΔzC(i, j, k, grid::RegularRectilinearGrid) = grid.Δz
@inline ΔzC(i, j, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.Δzᵃᵃᶠ[k]

@inline ΔzF(i, j, k, grid::RegularRectilinearGrid) = grid.Δz
@inline ΔzF(i, j, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.Δzᵃᵃᶜ[k]

@inline Δzᵃᵃᶠ(i, j, k, grid::RegularRectilinearGrid) = grid.Δz
@inline Δzᵃᵃᶠ(i, j, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.Δzᵃᵃᶠ[k]

@inline Δzᵃᵃᶜ(i, j, k, grid::RegularRectilinearGrid) = grid.Δz
@inline Δzᵃᵃᶜ(i, j, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.Δzᵃᵃᶜ[k]

#####
##### "Spacings" in Flat directions for rectilinear grids.
##### Here we dispatch all spacings to `one`. This abuse of notation
##### makes volumes and areas correct. E.g., with `z` Flat we have `v = dx * dy * 1`.
#####
##### Note: Vertical metrics are specific to each rectilinear grid type; therefore
##### we must dispatch on Flat vertical directions for each grid independently.
#####

using Oceananigans.Grids: Flat

#####
##### Horizontal metrics for AbstractRectilinearGrid
#####

const XFlatARG = AbstractRectilinearGrid{<:Any, <:Flat}
const YFlatARG = AbstractRectilinearGrid{<:Any, <:Any, <:Flat}

@inline Δx(i, j, k, grid::XFlatARG) = one(eltype(grid))
@inline Δy(i, j, k, grid::YFlatARG) = one(eltype(grid))

@inline Δxᶜᶜᵃ(i, j, k, grid::XFlatARG) = one(eltype(grid))
@inline Δxᶜᶠᵃ(i, j, k, grid::XFlatARG) = one(eltype(grid))
@inline Δxᶠᶠᵃ(i, j, k, grid::XFlatARG) = one(eltype(grid))
@inline Δxᶠᶜᵃ(i, j, k, grid::XFlatARG) = one(eltype(grid))

@inline Δyᶜᶜᵃ(i, j, k, grid::YFlatARG) = one(eltype(grid))
@inline Δyᶠᶜᵃ(i, j, k, grid::YFlatARG) = one(eltype(grid))
@inline Δyᶜᶠᵃ(i, j, k, grid::YFlatARG) = one(eltype(grid))
@inline Δyᶠᶠᵃ(i, j, k, grid::YFlatARG) = one(eltype(grid))

#####
##### Vertical metrics for RegularRectilinearGrid
#####

const ZFlatRRG = RegularRectilinearGrid{<:Any, <:Any, <:Any, <:Flat}

@inline ΔzC(  i, j, k, grid::ZFlatRRG) = one(eltype(grid))
@inline ΔzF(  i, j, k, grid::ZFlatRRG) = one(eltype(grid))
@inline Δzᵃᵃᶠ(i, j, k, grid::ZFlatRRG) = one(eltype(grid))
@inline Δzᵃᵃᶜ(i, j, k, grid::ZFlatRRG) = one(eltype(grid))

#####
##### Vertical metrics for VerticallyStretchedRectilinearGrid
#####

const ZFlatVSRG = VerticallyStretchedRectilinearGrid{<:Any, <:Any, <:Any, <:Flat}

@inline ΔzC(  i, j, k, grid::ZFlatVSRG) = one(eltype(grid))
@inline ΔzF(  i, j, k, grid::ZFlatVSRG) = one(eltype(grid))
@inline Δzᵃᵃᶠ(i, j, k, grid::ZFlatVSRG) = one(eltype(grid))
@inline Δzᵃᵃᶜ(i, j, k, grid::ZFlatVSRG) = one(eltype(grid))

#####
##### Areas for horizontally-regular algorithms
#####

@inline Axᵃᵃᶜ(i, j, k, grid) = Δy(i, j, k, grid) * ΔzF(i, j, k, grid)
@inline Axᵃᵃᶠ(i, j, k, grid) = Δy(i, j, k, grid) * ΔzC(i, j, k, grid)

@inline Ayᵃᵃᶜ(i, j, k, grid) = Δx(i, j, k, grid) * ΔzF(i, j, k, grid)
@inline Ayᵃᵃᶠ(i, j, k, grid) = Δx(i, j, k, grid) * ΔzC(i, j, k, grid)

@inline Azᵃᵃᵃ(i, j, k, grid) = Δx(i, j, k, grid) * Δy(i, j, k, grid)

#####
##### Volumes for horizontally-regular algorithms
#####

@inline Vᵃᵃᶜ(i, j, k, grid) = Δx(i, j, k, grid) * Δy(i, j, k, grid) * ΔzF(i, j, k, grid)
@inline Vᵃᵃᶠ(i, j, k, grid) = Δx(i, j, k, grid) * Δy(i, j, k, grid) * ΔzC(i, j, k, grid)

#####
##### Grid lengths for horizontally-curvilinear, vertically-rectilinear algorithms
#####

@inline Δxᶜᶜᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid)
@inline Δxᶜᶠᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid)
@inline Δxᶠᶠᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid)
@inline Δxᶠᶜᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid)

@inline Δyᶜᶜᵃ(i, j, k, grid::ARG) = Δy(i, j, k, grid)
@inline Δyᶠᶜᵃ(i, j, k, grid::ARG) = Δy(i, j, k, grid)
@inline Δyᶜᶠᵃ(i, j, k, grid::ARG) = Δy(i, j, k, grid)
@inline Δyᶠᶠᵃ(i, j, k, grid::ARG) = Δy(i, j, k, grid)

#####
##### Areas for algorithms that generalize to horizontally-curvilinear, vertically-rectilinear grids
#####

@inline Azᶜᶜᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid) * Δy(i, j, k, grid)
@inline Azᶠᶠᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid) * Δy(i, j, k, grid)
@inline Azᶜᶠᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid) * Δy(i, j, k, grid)
@inline Azᶠᶜᵃ(i, j, k, grid::ARG) = Δx(i, j, k, grid) * Δy(i, j, k, grid)

#####
##### Areas for three-dimensionally curvilinear algorithms
#####

@inline Axᶜᶜᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δyᶜᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # c
@inline Axᶠᶜᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δyᶠᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # u
@inline Axᶠᶠᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δyᶠᶠᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # ζ
@inline Axᶜᶠᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δyᶜᶠᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # η
@inline Axᶠᶜᶠ(i, j, k, grid::Union{ARG, AHCG}) = Δyᶠᶜᵃ(i, j, k, grid) * Δzᵃᵃᶠ(i, j, k, grid) # η
@inline Axᶜᶜᶠ(i, j, k, grid::Union{ARG, AHCG}) = Δyᶜᶜᵃ(i, j, k, grid) * Δzᵃᵃᶠ(i, j, k, grid) # η

@inline Ayᶜᶜᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δxᶜᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # c
@inline Ayᶜᶠᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δxᶜᶠᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # v
@inline Ayᶠᶜᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δxᶠᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # v
@inline Ayᶠᶠᶜ(i, j, k, grid::Union{ARG, AHCG}) = Δxᶠᶠᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid) # ζ
@inline Ayᶜᶠᶠ(i, j, k, grid::Union{ARG, AHCG}) = Δxᶜᶠᵃ(i, j, k, grid) * Δzᵃᵃᶠ(i, j, k, grid) # ξ
@inline Ayᶜᶜᶠ(i, j, k, grid::Union{ARG, AHCG}) = Δxᶜᶜᵃ(i, j, k, grid) * Δzᵃᵃᶠ(i, j, k, grid) # ξ

#####
##### Volumes for three-dimensionally curvilinear algorithms
#####

@inline Vᶜᶜᶜ(i, j, k, grid::Union{ARG, AHCG}) = Azᶜᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid)
@inline Vᶠᶜᶜ(i, j, k, grid::Union{ARG, AHCG}) = Azᶠᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid)
@inline Vᶜᶠᶜ(i, j, k, grid::Union{ARG, AHCG}) = Azᶜᶠᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid)
@inline Vᶜᶜᶠ(i, j, k, grid::Union{ARG, AHCG}) = Azᶜᶜᵃ(i, j, k, grid) * Δzᵃᵃᶠ(i, j, k, grid)

#####
##### Temporary place for grid spacings and areas for RegularLatitudeLongitudeGrid
#####

@inline hack_cosd(φ) = cos(π * φ / 180)
@inline hack_sind(φ) = sin(π * φ / 180)

@inline Δxᶜᶠᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶠᵃ[j]) * deg2rad(grid.Δλ)
@inline Δxᶠᶜᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶜᵃ[j]) * deg2rad(grid.Δλ)
@inline Δxᶜᶜᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = Δxᶠᶜᵃ(i, j, k, grid)
@inline Δxᶠᶠᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = Δxᶜᶠᵃ(i, j, k, grid)

@inline Δyᶜᶠᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.radius * deg2rad(grid.Δφ)
@inline Δyᶠᶜᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = Δyᶜᶠᵃ(i, j, k, grid)
@inline Δyᶜᶜᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = Δyᶜᶠᵃ(i, j, k, grid)
@inline Δyᶠᶠᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = Δyᶜᶠᵃ(i, j, k, grid)

@inline Δzᵃᵃᶜ(i, j, k, grid::RegularLatitudeLongitudeGrid) = grid.Δz
@inline Δzᵃᵃᶠ(i, j, k, grid::RegularLatitudeLongitudeGrid) = grid.Δz

@inline Azᶜᶜᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))
@inline Azᶠᶠᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλ) * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶠᶜᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = Azᶜᶜᵃ(i, j, k, grid)
@inline Azᶜᶠᵃ(i, j, k, grid::RegularLatitudeLongitudeGrid) = Azᶠᶠᵃ(i, j, k, grid)

#####
##### Temporary place for grid spacings and areas for ConformalCubedSphereFaceGrid
#####

@inline Δxᶜᶜᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δxᶜᶜᵃ[i, j]
@inline Δxᶠᶜᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δxᶠᶜᵃ[i, j]
@inline Δxᶜᶠᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δxᶜᶠᵃ[i, j]
@inline Δxᶠᶠᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δxᶠᶠᵃ[i, j]

@inline Δyᶜᶜᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δyᶜᶜᵃ[i, j]
@inline Δyᶠᶜᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δyᶠᶜᵃ[i, j]
@inline Δyᶜᶠᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δyᶜᶠᵃ[i, j]
@inline Δyᶠᶠᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Δyᶠᶠᵃ[i, j]

@inline Δzᵃᵃᶜ(i, j, k, grid::ConformalCubedSphereFaceGrid) = grid.Δz
@inline Δzᵃᵃᶠ(i, j, k, grid::ConformalCubedSphereFaceGrid) = grid.Δz

@inline Azᶜᶜᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Azᶜᶜᵃ[i, j]
@inline Azᶠᶜᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Azᶠᶜᵃ[i, j]
@inline Azᶜᶠᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Azᶜᶠᵃ[i, j]
@inline Azᶠᶠᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.Azᶠᶠᵃ[i, j]

#####
##### Generic functions for specified locations
#####
##### For example, Δx(i, j, k, Face, Center, LZ) is equivalent to = Δxᶠᶜᵃ(i, j, k, grid).
#####
##### We also use the function "volume" rather than `V`.
#####

location_code_xy(LX, LY) = Symbol(interpolation_code(LX), interpolation_code(LY), :ᵃ)
location_code(LX, LY, LZ) = Symbol(interpolation_code(LX), interpolation_code(LY), interpolation_code(LZ))

for LX in (:Center, :Face)
    for LY in (:Center, :Face)
        LXe = @eval $LX
        LYe = @eval $LY

        Ax_function = Symbol(:Ax, location_code(LXe, LYe, Center()))
        Ay_function = Symbol(:Ay, location_code(LXe, LYe, Center()))
        Az_function = Symbol(:Az, location_code_xy(LXe, LYe))

        @eval begin
            Az(i, j, k, grid, ::$LX, ::$LY, LZ) = $Az_function(i, j, k, grid)
            Ax(i, j, k, grid, ::$LX, ::$LY, ::Center) = $Ax_function(i, j, k, grid)
            Ay(i, j, k, grid, ::$LX, ::$LY, ::Center) = $Ay_function(i, j, k, grid)
        end
    end
end

Ax(i, j, k, grid, ::Face, ::Center, ::Face) = Axᶠᶜᶠ(i, j, k, grid)
Ay(i, j, k, grid, ::Center, ::Face, ::Face) = Ayᶜᶠᶠ(i, j, k, grid)

volume(i, j, k, grid, ::Center, ::Center, ::Center) = Vᶜᶜᶜ(i, j, k, grid)
volume(i, j, k, grid, ::Face,   ::Center, ::Center) = Vᶠᶜᶜ(i, j, k, grid)
volume(i, j, k, grid, ::Center, ::Face,   ::Center) = Vᶜᶠᶜ(i, j, k, grid)
volume(i, j, k, grid, ::Center, ::Center, ::Face)   = Vᶜᶜᶠ(i, j, k, grid)


text:
function ADbackward!(net;debug=false,AccumulateGradient=false)

    #=   calculate the gradient of a net function
    (c) David Barber, University College London 2015
    =#

    if debug; println("Get gradient:"); end

    fill!(net.gradient[net.FunctionNode],1.0)
    for n in net.ancestors
        if net.node[n].takederivative
            if !AccumulateGradient
                fill!(net.gradient[n],0.0)
            end
            for c in net.relevantchildren[n]
                if debug
                    println("-----------------")
                    println("node $n, child $c: $(net.node[c].df)($(net.node[c].parents))")
                     for parIDX in net.parentIDX[c,n] # deals with the case that x-->f<--x (which happens with e.g. x+x)
                        println("parent index [$parIDX]")
                        # This is why gradients must always add up whatever is currently there (ie not replace).
                        @time net.node[c].df(parIDX,net.value[c],net.auxvalue[c],net.gradient[c],net.gradient[n],net.value[net.node[c].parents]...)
                    end
                else
                    for parIDX in net.parentIDX[c,n] # deals with the case that x-->f<--x (which happens with e.g. x+x)
                        net.node[c].df(parIDX,net.value[c],net.auxvalue[c],net.gradient[c],net.gradient[n],net.value[net.node[c].parents]...)
                    end
                end
            end
        end
    end

    if debug; println("done backward pass");  end

end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff