text:
"""
## shape_fun!

This subroutine produces derivatives of shape functions with respect
to local coordinates.

### Function
```julia
shape_fun!(fun, point, i)
```

### Arguments
```julia
* fun::Vector{Float64}       : Shape function (Updated)
* points::Matrix{Float64}    : Local coordinates of integration points
* i::Int                   : Integration point
```
"""
function shape_fun!(fun::Vector{Float64}, points::Matrix{Float64}, i::Int)
  #
  #   This subroutine computes the values of the shape functions.
  #   to local coordinates
  #

  ndim = size(points,2)
  nod = length(fun)
  if ndim == 1                                # 1.0 dimensional case
    xi=points[i,1]
    if nod == 1
      t1 = -1.0-xi
      t2 = 1.0-xi
      fun[1] = t2/2.0
      fun[2] = -t1/2.0
    elseif nod == 3
      t1=-1.0-xi
      t2=-xi
      t3=1.0-xi
      fun[1] = t2*t3/2.0
      fun[2] = -t1*t3
      fun[3] = t1*t2/2.0
    elseif nod == 4
      t1=-1.0-xi
      t2=-1.0/3.0-xi
      t3=1.0/3.0-xi
      t4=1.0-xi
      fun[1] = t2*t3*t4*9.0/16.0
      fun[2] = -t1*t3*t4*27.0/16.0
      fun[3] = t1*t2*t4*27.0/16.0
      fun[4] = -t1*t2*t3*9.0/16.0
    elseif nod == 5
      t1=-1.0 -xi
      t2=-0.5-xi
      t3=-xi
      t4=0.5-xi
      t5=1.0-xi
      fun[1] = t2*t3*t4*t5*2.0/3.0
      fun[2] = -t1*t3*t4*t5*8.0/3.0
      fun[3] = t1*t2*t4*t5*4.0
      fun[4] = -t1*t2*t3*t5*8.0/3.0
      fun[5] = t1*t2*t3*t4*2.0/3.0
    else
      println("Wrong number of nodes in shape_fun.")
    end
  elseif ndim == 2                                # 2.0 dimensional case
    c1=points[i,1]
    c2=points[i,2]
    c3=1.0-c1-c2
    xi=points[i,1]
    eta=points[i,2]
    etam=0.25*(1.0-eta)
    etap=0.25*(1.0+eta)
    xim=0.25*(1.0-xi)
    xip=0.25*(1.0+xi)
    if nod == 3
      fun[:] = [c1,c3,c2]
    elseif nod == 4
      fun[:] = [4.0*xim*etam,4.0*xim*etap,4.0*xip*etap,4.0*xip*etam]
    elseif nod == 5
      fun[:] = [4.0*xim*etam-0.25*(1.0-xi*xi)*(1.0-eta*eta),
            4.0*xim*etap-0.25*(1.0-xi*xi)*(1.0-eta*eta),
            4.0*xip*etap-0.25*(1.0-xi*xi)*(1.0-eta*eta),
            4.0*xip*etam-0.25*(1.0-xi*xi)*(1.0-eta*eta),
            (1.0-xi*xi)*(1.0-eta*eta)]
    elseif nod == 6
      fun[1] = (2.0*c1-1.0)*c1
      fun[2] = 4.0*c3*c1
      fun[3] = (2.0*c3-1.0)*c3
      fun[4] = 4.0*c2*c3
      fun[5] = (2.0*c2-1.0)*c2
      fun[6] = 4.0*c1*c2
    elseif nod == 8
      fun[:] = [4.0*etam*xim*(-xi-eta-1.0),32.0*etam*xim*etap,
            4.0*etap*xim*(-xi+eta-1.0),32.0*xim*xip*etap,
            4.0*etap*xip*(xi+eta-1.0), 32.0*etap*xip*etam,
            4.0*xip*etam*(xi-eta-1.0), 32.0*xim*xip*etam]
    elseif nod == 9
      etam=eta-1.0
      etap=eta+1.0
      xim=xi-1.0
      xip=xi+1.0
      fun[:] = [0.25*xi*xim*eta*etam,-0.5*xi*xim*etap*etam,
            0.25*xi*xim*eta*etap,-0.5*xip*xim*eta*etap,
            0.25*xi*xip*eta*etap,-0.5*xi*xip*etap*etam,
            0.25*xi*xip*eta*etam,-0.5*xip*xim*eta*etam,
            xip*xim*etap*etam]
    elseif nod == 10
      fun[1] =  ((3.0*c1-1.0)*(3.0*c1-2.0)*c1)/2.0
      fun[2] =  -(9.0*(3.0*c1-1.0)*(c1+c2-1.0)*c1)/2.0
      fun[3] =   (9.0*(3.0*c1+3.0*c2-2.0)*(c1+c2-1.0)*c1)/2.0
      fun[4] = -((3.0*c1+3.0*c2-1.0)*(3.0*c1+3.0*c2-2.0)*(c1+c2-1.0))/2.0
      fun[5] =   (9.0*(3.0*c1+3.0*c2-2.0)*(c1+c2-1.0)*c2)/2.0
      fun[6] =  -(9.0*(c1+c2-1.0)*(3.0*c2-1.0)*c2)/2.0
      fun[7] =  ((3.0*c2-1.0)*(3.0*c2-2.0)*c2)/2.0
      fun[8] =   (9.0*(3.0*c2-1.0)*c1*c2)/2.0
      fun[9] =   (9.0*(3.0*c1-1.0)*c1*c2)/2.0
      fun[10] = -27.0*((c2-1.0)+c1)*c1*c2
    elseif nod == 15
      t1=c1-0.25
      t2=c1-0.5
      t3=c1-0.75
      t4=c2-0.25
      t5=c2-0.5
      t6=c2-0.75
      t7=c3-0.25
      t8=c3-0.5
      t9=c3-0.75
      fun[1] = 32.0/3.0*c1*t1*t2*t3
      fun[2] = 128.0/3.0*c3*c1*t1*t2
      fun[3] = 64.0*c3*c1*t1*t7
      fun[4] = 128.0/3.0*c3*c1*t7*t8
      fun[5] = 32.0/3.0*c3*t7*t8*t9
      fun[6] = 128.0/3.0*c2*c3*t7*t8
      fun[7] = 64.0*c2*c3*t4*t7
      fun[8] = 128.0/3.0*c2*c3*t4*t5
      fun[9] = 32.0/3.0*c2*t4*t5*t6
      fun[10] = 128.0/3.0*c1*c2*t4*t5
      fun[11] = 64.0*c1*c2*t1*t4
      fun[12] = 128.0/3.0*c1*c2*t1*t2
      fun[13] = 128.0*c1*c2*t1*c3
      fun[15] = 128.0*c1*c2*c3*t4
      fun[14] = 128.0*c1*c2*c3*t7
    else
      println("Wrong number of nodes in shape_fun.")
    end
  elseif ndim == 3                            # 3.0 dimensional case
    xi=points[i,1]
    eta=points[i,2]
    zeta=points[i,3]
    etam=1.0-eta
    xim=1.0-xi
    zetam=1.0-zeta
    etap=eta+1.0
    xip=xi+1.0
    zetap=zeta+1.0
    if nod == 4
     fun[1] = xi
     fun[2] = eta
     fun[3] = zeta
     fun[4] = 1.0-fun[1]-fun[2]-fun[3]
    elseif nod == 8
     fun[:] = [0.125*xim*etam*zetam,0.125*xim*etam*zetap,
            0.125*xip*etam*zetap,0.125*xip*etam*zetam,
            0.125*xim*etap*zetam,0.125*xim*etap*zetap,
            0.125*xip*etap*zetap,0.125*xip*etap*zetam]
    elseif nod == 14                          # type 6 fin_el
     fun[1] = (xi*eta+xi*zeta+2.0*xi+eta*zeta+2.0*eta+2.0*zeta+2.0)*
       (xi-1.0)*(eta-1.0)*(zeta-1.0)/8.0
     fun[2] =-(xi*eta-xi*zeta+2.0*xi-eta*zeta+2.0*eta-2.0*zeta+2.0)*
       (xi-1.0)*(eta-1.0)*(zeta+1.0)/8.0
     fun[3] =-(xi*eta-xi*zeta+2.0*xi+eta*zeta-2.0*eta+2.0*zeta-2.0)*
       (xi+1.0)*(eta-1.0)*(zeta+1.0)/8.0
     fun[4] = (xi*eta+xi*zeta+2.0*xi-eta*zeta-2.0*eta-2.0*zeta-2.0)*
       (xi+1.0)*(eta-1.0)*(zeta-1.0)/8.0
     fun[5] =-(xi+1.0)*(xi-1.0)*(eta-1.0)*(zeta+1.0)*(zeta-1.0)/2.0
     fun[6] =-(xi-1.0)*(eta+1.0)*(eta-1.0)*(zeta+1.0)*(zeta-1.0)/2.0
     fun[7] = (xi+1.0)*(xi-1.0)*(eta+1.0)*(eta-1.0)*(zeta+1.0)/2.0
     fun[8] = (xi+1.0)*(eta+1.0)*(eta-1.0)*(zeta+1.0)*(zeta-1.0)/2.0
     fun[9] =-(xi+1.0)*(xi-1.0)*(eta+1.0)*(eta-1.0)*(zeta-1.0)/2.0
     fun[10] =  (xi*eta-xi*zeta-2.0*xi+eta*zeta+2.0*eta-2.0*zeta-2.0)*
       (xi-1.0)*(eta+1.0)*(zeta-1.0)/8.0
     fun[11] = -(xi*eta+xi*zeta-2.0*xi-eta*zeta+2.0*eta+2.0*zeta-2.0)*
       (xi-1.0)*(eta+1.0)*(zeta+1.0)/8.0
     fun[12] = -(xi*eta+xi*zeta-2.0*xi+eta*zeta-2.0*eta-2.0*zeta+2.0)*
       (xi+1.0)*(eta+1.0)*(zeta+1.0)/8.0
     fun[13] =  (xi*eta-xi*zeta-2.0*xi-eta*zeta-2.0*eta+2.0*zeta+2.0)*
       (xi+1.0)*(eta+1.0)*(zeta-1.0)/8.0
     fun[14] =  (xi+1.0)*(xi-1.0)*(eta+1.0)*(zeta+1.0)*(zeta-1.0)/2.0
    elseif nod == 20
      xii = [-1,-1,-1,0,1,1,1,0,-1,-1,1,1,-1,-1,-1,0,1,1,1,0]
      etai = [-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,1,1,1]
      zetai = [-1,0,1,1,1,0,-1,-1,-1,1,1,-1,-1,0,1,1,1,0,-1,-1]
      for l in 1:20
        xi0 = xi*xii[l]
        eta0 = eta*etai[l]
        zeta0 = zeta*zetai[l]
        if l == 4 || l == 8 || l == 16 || l == 20
          fun[l] = 0.25*(1.0-xi*xi)*(1.0+eta0)*(1.0+zeta0)
        elseif l >= 9 && l <= 12
          fun[l] = 0.25*(1.0+xi0)*(1.0-eta*eta)*(1.0+zeta0)
        elseif l == 2 || l == 6 || l == 14 || l == 18
          fun[l] = 0.25*(1.0+xi0)*(1.0+eta0)*(1.0-zeta*zeta)
        else
          fun[l] = 0.125*(1.0+xi0)*(1.0+eta0)*(1.0+zeta0)*(xi0+eta0+zeta0-2)
        end
      end
    else
      println("Wrong number of nodes in shape_fun.")
    end
  else
   println("Wrong number of dimensions in shape_fun.")
  end
end


text:
using SDE
using Base.Test
require(Pkg.dir("SDE","src", "Randm.jl"))
#include(Pkg.dir("SDE","src", "quad.jl"))
import SDE.eps2

function mc(M::Vector)
  m = mean(M)
  ste = std(M)/sqrt(length(M))
  m, ste
end
qu = sqrt(2)*erfinv(2*0.995-1)

macro repeat(ex, n)
    quote
       r1 = $ex
       y = zeros(size(r1)..., $n)
       s = prod(size(r1))
       y[1:s] = r1
       for i in 1:$n-1
        y[i*s+1:(i+1)*s] = $ex
       end
       y
    end
end

for d in [1:3]
    N = 5000
    xi = 0.3
	global T, t, v, h, P, B, A, beta, s
	println("Test dimension $d xi $xi")

	if (d== 1) srand(7)
	elseif (d== 2) srand(8)
	else srand(6)
	end
	B = xi.*Randm.randstable(d)
#	A = Randm.randposdef(d)
	Si = rand(d,d)
	A = Si*Si'
	beta = 0.5*randn(d)

	P = MvLinPro(B, beta, Si)
	P0 = MvAffPro(beta, Si)
	P00 = MvLinPro(eye(d)*1E-7, beta, Si)

	function aa(s)
	 A
	end

	T = 1.7
	t = 0.5
	v = zeros(d)

	println("Test SDE.H, SDE.K")

    H1 = SDE.H(t,T, P)
    H2 = SDE.varH(t, T, B, aa)
    K = SDE.K(t, T, P)
    H3 = expm((T-t)*B')*inv(K)*expm((T-t)*B)

    H0 = SDE.H(t, T, P0)
    K0 = SDE.K(t, T, P0)
    H00 = SDE.H(t, T, P00)
    @test norm(H0 - H00) < 1E-5*d^2

	h = T-t
	phi = expm(h*B)
	phim = expm(-h*B)
	K2 = P.lambda - phi*P.lambda*phi'
	H4 = phi'*inv(P.lambda - phi*P.lambda*phi')*phi

	@test norm(H1 - H2) < 1E-10
	@test norm(H1 - H3) < 1E-10
	@test norm(H1 - H4) < 1E-10

	@test norm(K-K2) < d*24*eps()

	function varr3(t, x, T, v, P)
	    h = T - t
		binv = inv(P.B)
		phi = expm(h*P.B)
		phim = expm(-h*P.B)
		vt =  phim*( v) + phim * binv * beta -binv*beta
		SDE.H(t, T, P) * (vt-x )

	end

	println("Test SDE.r")

	srand(5)
	x0 = x = randn(d)/d
	v = randn(d)
	r1 = SDE.r(t, x, T, v, P)
	r2 = SDE.varr(t, x, T, v, P)
	@test norm(r1 - r2) < 1E-10

	r3 = varr3(t, x, T, v, P)
	@test norm(r1 - r3) < 1E-10

	r10 = SDE.r(t, x, T, v, P0)
	r100 = SDE.r(t, x, T, v, P00)
	@test norm(r10 - r100) < 1E-5*d^2

	println("Test SDE.mu")

	mu = SDE.mu(t, x, T, P)
	mu0 = SDE.mu(t, x, T, P0)
	mu2 = SDE.varmu(t, x, T, P)
	@test norm(mu - mu2) < 1E-13

    println("Testing mean and variance of samplep")

    y = @repeat(SDE.samplep(t, x, T, P), N)
    mcy = mc(sum(y,1)[:])
	println(repr(mcy), sum(mu))
    @test abs(mcy[1]- sum(mu)) < qu*mcy[2]
    muhat = mean(y,2)
	Khat = cov((y.-x)')
	println(norm(Khat - K))
	@test norm(Khat-K) < 0.04*d

    y = @repeat(SDE.samplep(t, x, T, P0), N)
    mcy = mc(sum(y,1)[:])
	println(repr(mcy), sum(mu0))
    @test abs(mcy[1]- sum(mu0)) < qu*mcy[2]
    muhat = mean(y,2)
	Khat = cov((y.-x)')
	println(norm(Khat - K0))
	@test norm(Khat-K0) < 0.04*d

	#checks if transition density p integrates to 1 by monte carlo integration
	println("Test SDE.lp")

	function testp(h, x0, N, B, beta, Si)
	    local P, P0
		mu =  B*x0+beta
		P = MvLinPro(B, beta, Si)
		P0 = MvAffPro(mu, Si)
    	# monte carlo integration of the transition density p with proposals according to p0
		w1 = zeros(N)
		for n in 1:N
			y = SDE.samplep(0., x0, h, P0 )
			w1[n] =  exp(SDE.lp(0., x0, h,  y, P)-SDE.lp(0., x0, h, y, P0))
		end
		# monte carlo integration of the transition density p0 with proposals according to p
	    w2 = zeros(N)
		for n in 1:N
			y = SDE.samplep(0., x0, h, P )
			w2[n] =  exp(SDE.lp(0., x0, h,  y, P0)-SDE.lp(0., x0, h, y, P))
		end

		return mc(w1), mc(w2), w1, w2
	end

	srand(5)
	n1, n2, w1, w2 = testp(0.05, x, N, B, beta, Si)
	println(repr(n1), repr(n2))
	@test abs(n1[1]-1.) < qu*n1[2]
	@test abs(n2[1]-1.) < qu*n2[2]

	@test abs(SDE.lp(t, x,T, v, P) - SDE.varlp(t, x,T, v, (t,s) -> t-s, B, beta, s->A)) < 1E-10

    tmin = t/2
    tmax = T
    T = tmax-tmin

	s = SDE.soft(t, tmin, T)
	@test norm(t - SDE.tofs(s, tmin, T)) < 1E-13
	@test norm(t - (tmax - (T - s)^2/T)) < 1E-13

	v1 = SDE.Vs(s, T, v, P0)
	v2 = SDE.V(t, tmax, v, P0)
#    println("v10, v20", repr(v1), repr(v2))
	@test norm(v1-v2) < d^2*1E-13

	v1 = SDE.Vs(s, T, v, P)
	v2 = SDE.V(t, tmax, v, P)
#    println("v1, v2", repr(v1), repr(v2))
 	@test norm(v1-v2) < d^2*1E-13

    v1 = SDE.dotVs(s, T, v, P0)
    v2 = SDE.dotVs(s, T, v, P00)
    println ("dotv ", repr(v1) ,repr(v2))
    @test norm(v1 - v2) < 1E-5

	us =  SDE.uofx(s, x, T, v,  P)
	xt =  SDE.xofu(s, us, T, v,  P)
	@test norm(x - xt) < 1E-10

	r1 =  SDE.r(t, x, tmax, v, P)
 	r2 =  SDE.J(s, T, P) *us*T/(T-s)
	@test norm(r1 .- r2) < 1E-10

end

ext: jl


text:
function tchanelSolve(
    p0,
    T,
    IRScale,
    UVScale,
    massfun,
    Imlambdaflowfun,
    Relambdaflowfun,
    ϵ = 0.1;
    u0,
    progress = true,
    atol = 1e-14,
    rtol = 1e-14,
    adaptive = true,
    dtmax = 0.01,
    kwargs...,
)
    function flow_t(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = Imlambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
        du[2] = Relambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
    end
    alg = Tsit5()
    kspan = (UVScale, IRScale)
    prob = ODEProblem(flow_t, u0, kspan, p0)
    solve(
        prob,
        alg,
        progress = progress,
        atol = atol,
        rtol = rtol,
        adaptive = adaptive,
        dtmax = dtmax,
        kwargs...,
    )
end

function tchanelSolveParallel(
    T,
    IRScale,
    UVScale,
    massfun,
    Imlambdaflowfun,
    Relambdaflowfun, numberOfParameters,
    ϵ = 0.1;
    u0,
    atol = 1e-14,
    rtol = 1e-14,
    adaptive = true,
    dtmax = 0.01,
    kwargs...
)
    function flow_t(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = Imlambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
        du[2] = Relambdaflowfun(p, k, massfun(k), T, u[1], u[2], ϵ)
    end
    alg = Tsit5()
    kspan = (UVScale, IRScale)
    p = 1.0
    prob = ODEProblem(flow_t, u0, kspan, p)
    parameterList = range(1.0, stop = 500.0, length = numberOfParameters)
    # function parameterChange(prob, i, repeat)
    #     prob.p = parameterList[i]
    #     prob
    # end
    parameterChange = (prob, i, repeat) -> remake(prob, p = parameterList[i])
    ensembleProb = EnsembleProblem(prob, prob_func = parameterChange)
    solve(
        ensembleProb,
        alg,
        EnsembleThreads(),
        trajectories = numberOfParameters,
        atol = atol,
        rtol = rtol,
        adaptive = adaptive,
        dtmax = dtmax,
        kwargs...,
    )
end

function tchanelSolveTwoPointParallel(
    T,
    IRScale,
    UVScale,
    massfun,
    lamdafun,
    Imlambdaflowfun,
    Relambdaflowfun, numberOfParameters,
    ϵ = 0.1;
    u0,
    atol = 1e-14,
    rtol = 1e-14,
    adaptive = true,
    dtmax = 0.01,
    kwargs...
)
    function flow_t(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = Imlambdaflowfun(p, k, massfun(k),lamdafun(k), T, u[1], u[2], ϵ)
        du[2] = Relambdaflowfun(p, k, massfun(k),lamdafun(k), T, u[1], u[2], ϵ)
    end
    alg = Tsit5()
    kspan = (UVScale, IRScale)
    p = 0.0
    prob = ODEProblem(flow_t, u0, kspan, p)
    parameterList = range(0, stop = 400.0, length = numberOfParameters)
    # function parameterChange(prob, i, repeat)
    #     prob.p = parameterList[i]
    #     prob
    # end
    parameterChange = (prob, i, repeat) -> remake(prob, p = parameterList[i])
    function imrefun(sol,i)
        p0=parameterList[i];
        imgamm2=-hcubature(
            x -> twopointflow2(p0, x[1], massfun(x[1]), T, y->sol(y)[1]),
            [1.0],
            [800.0],
            atol = 1e-6,
            rtol = 1e-6,
        )[1];
        regamm2=hcubature(
            x -> twopointflow2(p0, x[1], massfun(x[1]), T,  y->sol(y)[2]),
            [1.0],
            [800.0],
            atol = 1e-6,
            rtol = 1e-6,
        )[1]  +p0^2 - massfun(800.0);
        return ([imgamm2,regamm2],false)
    end
    ensembleProb = EnsembleProblem(prob, prob_func = parameterChange,output_func=imrefun)
    solve(
        ensembleProb,
        alg,
        EnsembleThreads(),
        trajectories = numberOfParameters,
        atol = atol,
        rtol = rtol,
        adaptive = adaptive,
        dtmax = dtmax,
        kwargs...,
    )
end


text:
"""
    test1_setup()

Call `gcmgrid`, initialize a single point,
rely on `dxy_dt`, and just output `sol` at the end.

```
using IndividualDisplacements, MeshArrays, OrdinaryDiffEq
𝑃,sol=test1_setup()
```
"""
function test1_setup()

    mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write)
    XC=MeshArray(mygrid,Float32); XC[1]=vec(2500.:5000.:397500.0)*ones(1,42);
    XG=MeshArray(mygrid,Float32); XG[1]=vec(0.:5000.:395000.0)*ones(1,42);
    YC=MeshArray(mygrid,Float32); YC[1]=ones(80,1)*transpose(vec(2500.:5000.:207500.0));
    YG=MeshArray(mygrid,Float32); YG[1]=ones(80,1)*transpose(vec(0.:5000.:205000.0));

    dx=5000.0
    t0=0.0; t1=18001.0*3600.0
    u=-(YG.-YC[1][40,21])/2000000.
    v=(XG.-XC[1][40,21])/2000000.
    u0=u./dx; u1=u./dx
    v0=v./dx; v1=v./dx

    𝑃 = (u0=u0, u1=u1, v0=v0, v1=v1, 𝑇=[t0,t1], XC=XC, YC=YC, ioSize=(80,42))
    u0=[200000.0;0.0]./dx
    du=fill(0.0,2);
    prob = ODEProblem(dxy_dt,u0,[0.0,2998*3600.0],𝑃)
    sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

    return 𝑃,sol
end

"""
    test2_periodic_domain(np = 12, nq = 12)

Call `simple_periodic_domain`, initialize 6x6 point cloud,
rely on `dxy_dt!`, and call `postprocess_xy` at the end.

```
using IndividualDisplacements, MeshArrays, OrdinaryDiffEq
df,𝑃=test2_periodic_domain()

using Plots
@gif for t in 𝑃.t0:1.0:𝑃.t1
   scatter_subset(𝑃,df,t)
end
```
"""
function test2_periodic_domain(np = 12, nq = 12)
    #domain and time parameters
    Γ = simple_periodic_domain(np, nq)
    Γ = IndividualDisplacements.dict_to_nt(Γ)

    u = 0.1 ./ Γ.DXC
    v = 0.3 ./ Γ.DYC
    (u, v) = exchange(u, v, 1)
    𝑃 = (u0=u, u1=u, v0=v, v1=v, 𝑇=[0.0,400.0], dt=0.1,
         XC=Γ.XC, YC=Γ.YC, ioSize=(np,nq))

    #initial conditions
    x0 = np * (0.4:0.04:0.6)
    y0 = nq * (0.4:0.04:0.6)
    x0 = vec(x0) * ones(1, length(y0))
    y0 = ones(size(x0, 1), 1) * transpose(vec(y0))
    u0 = permutedims([[x0[i];y0[i];1.0] for i in eachindex(x0)])
    du=0*u0

    #solve for trajectories
    prob = ODEProblem(dxy_dt!, u0, 𝑃.𝑇, 𝑃)
    sol = solve(prob,Euler(),dt=𝑃.dt)

    return postprocess_xy(sol, 𝑃),𝑃
end

"""
    scatter_subset(𝑃,df,t=missing,dt=1.0)

```
@gif for t in 𝑃.t0:1.0:𝑃.t1
   scatter_subset(𝑃,df,t)
end
```
"""
function scatter_subset(𝑃,df,t=missing,dt=1.0)
    ismissing(t) ? t=maximum(df[!,:t]) : nothing
    df_t = df[ (df.t.>t-dt).&(df.t.<=t) , :]
    nx,ny=size(𝑃.XC[1])
    scatter(df_t.x,df_t.y,markersize=2.0,c=:red,
    xlims=(0,nx),ylims=(0,ny),leg=:none,marker = (:circle, stroke(0)))
end


text:
Base.show(io::IO, ::MIME"text/markdown", markdown::OrgDoc) = (markdown(io, markdown), nothing)

function markdown(io::IO, o::OrgDoc)
    for component in o.contents
        markdown(io, component)
    end
end

markdown(o::Union{OrgDoc, Component}) = markdown(stdout, o)

markdown(io::IO, c::Component) = html(io, c)

function markdown(io::IO, heading::Heading)
    print(io, '#'^(heading.level + 1), ' ')
    if !isnothing(heading.keyword)
        print(io, heading.keyword, ' ')
    end
    if !isnothing(heading.priority)
        print(io, "[#", heading.priority, "] ")
    end
    markdown.(Ref(io), heading.title)
    print(io, '\n')
    if !isnothing(heading.section)
        print(io, '\n')
        markdown(io, heading.section)
    end
end

function markdown(io::IO, section::Section)
    for component in section.contents
        markdown(io, component)
        component === last(section.contents) || print(io, "\n\n")
    end
    print(io, '\n')
end

function markdown(io::IO, afkw::AffiliatedKeywordsWrapper)
    @nospecialize afkw
    markdown(io, afkw.element)
end

# ---------------------
# Greater Elements
# ---------------------

function markdown(io::IO, specialblock::SpecialBlock)
    print(io, "!!! ", specialblock.name, '\n')
    content = join(sprint.(markdown, specialblock.content), '\n')
    for line in split(content, '\n')
        print(io, "    ", line, '\n')
    end
end

function markdown(io::IO, centerb::CenterBlock)
    print(io, "
\n")
    for el in centerb.contents
        markdown(io, el)
        print(io, '\n')
    end
    print(io, "
\n")
end

function markdown(io::IO, quoteb::QuoteBlock)
    content = join(sprint.(markdown, quoteb.contents), '\n')
    for line in split(content, '\n')
        print(io, "> ", line, '\n')
    end
end

function markdown(io::IO, drawer::Drawer)
    for el in drawer.contents
        markdown(io, el)
        el === last(drawer.contents) || print(io, '\n')
    end
end

function markdown(io::IO, dynblock::DynamicBlock)
    for el in dynblock.contents
        markdown(io, el)
        el === last(dynblock.contents) || print(io, '\n')
    end
end

function markdown(io::IO, fn::FootnoteDefinition)
    print(io, "[^", fn.label, "]:\n\n")
    content = join(sprint.(markdown, fn.definition), '\n')
    for line in split(content, '\n')
        print(io, "    ", line, '\n')
    end
end

# InlineTask

function markdown(io::IO, item::Item, indent::Integer=0, offset::Integer=0)
    print(io, ' '^indent,
          if item.bullet in ("+", "-", "*")
              '-' else item.bullet end,
          ' ')
    offset += length(item.bullet) + 1
    if !isnothing(item.counterset)
        print(io, " [@", item.counterset, "]")
        offset += length(item.counterset) + 4
    end
    if !isnothing(item.checkbox)
        print(io, " [", item.checkbox, ']')
        offset += 4
    end
    if !isnothing(item.tag)
        tagbuf = IOContext(IOBuffer(), :color => get(io, :color, false),
                           :displaysize => (displaysize(io)[1],
                                            displaysize(io)[2] - indent - 2))
        for obj in item.tag
            markdown(tagbuf, obj)
        end
        taglines = wraplines(String(take!(tagbuf.io)),
                             displaysize(io)[2] - indent - 2, offset)
        for line in taglines
            print(io, line)
            line === last(taglines) || print(io, '\n', ' '^(indent+2))
        end
        print(io, " ::")
        offset += length(taglines[end]) + 3
    end
    if length(item.contents) > 0
        print(io, ' ')
        offset += 1
        contentbuf = IOContext(IOBuffer(), :color => get(io, :color, false),
                               :displaysize => (displaysize(io)[1],
                                                displaysize(io)[2] - indent - 2))
        parlines = if item.contents[1] isa Paragraph
            for obj in item.contents[1]; markdown(contentbuf, obj) end
            contents = String(take!(contentbuf.io))
            wraplines(contents, displaysize(io)[2] - indent - 2, offset)
        else
            [""]
        end
        components = @view item.contents[if item.contents[1] isa Paragraph 2 else 1 end:end]
        for component in components
            markdown(contentbuf, component, indent)
            component === last(components) || print(contentbuf, '\n')
        end
        otherlines = split(String(take!(contentbuf.io)), '\n')
        lines = vcat(parlines, if otherlines == [""]; [] else otherlines end)
        for line in lines
            print(io, line)
            line === last(lines) || print(io, '\n', ' '^(indent+2))
        end
    end
end

function markdown(io::IO, list::List, indent::Integer=0)
    @nospecialize list
    for item in list.items
        markdown(io, item, indent)
        item === last(list.items) || print(io, '\n')
    end
end

markdown(::IO, ::PropertyDrawer) = nothing

const table_charset_markdown =
    Dict('|' => '|',
         '>' => '|',
         '<' => '|',
         '[' => ' ',
         ']' => ' ',
         '-' => '-',
         '+' => '|')

markdown(io::IO, table::Table) =
    layouttable(io, markdown, table, table_charset_markdown, 0)

# ---------------------
# Elements
# ---------------------

markdown(::IO, ::BabelCall) = nothing

markdown(::IO, ::CommentBlock) = nothing

function markdown(io::IO, block::ExampleBlock)
    print(io, "```\n",
          join(block.contents, '\n'),
          "\n```")
end

function markdown(io::IO, block::SourceBlock)
    print(io, "```", block.lang, "\n",
          join(block.contents, '\n'),
          "\n```")
end

# Block
# Clock
# Planning
# DiarySexp

markdown(::IO, ::Comment) = nothing

function markdown(io::IO, fw::FixedWidth)
    print(io, "```\n",
          join(fw.contents, '\n'),
          "\n```")
end

markdown(io::IO, ::HorizontalRule) = print(io, "---")

function markdown(io::IO, keyword::Keyword)
    if keyword.key == "title"
        print(io, "# ", keyword.value)
    elseif keyword.key == "subtitle"
        print(io, "### ", keyword.value)
    elseif keyword.key in ("author", "date")
        print(io, "By ", keyword.value)
    end
end

function markdown(io::IO, env::LaTeXEnvironment)
    print(io, "\$\$\n\\begin{", env.name, "}\n",
          join(env.contents, '\n'),
          "\n\\end{", env.name, "}\n\$\$")
end

# Node Property

function markdown(io::IO, par::Paragraph)
    for obj in par.contents
        markdown(io, obj)
    end
end

markdown(io::IO, entity::Entity) = print(io, Entities[entity.name].utf8)

markdown(io::IO, latex::LaTeXFragment) =
    print(io, "\$", latex.contents, "\$")

markdown(io::IO, snippet::ExportSnippet) =
    if snippet.backend == "markdown" print(io, snippet.snippet) end

function markdown(io::IO, fn::FootnoteReference)
    label = something(fn.label, string(rand(UInt32), base=36))
    print(io, "[^", label, "]")
    if !isnothing(fn.definition)
        print(io, '\n')
        markdown(io, FootnoteDefinition(
            SubString(label), Element[Paragraph(fn.definition)]))
        print(io, '\n')
    end
end

function markdown(io::IO, keycite::CitationReference)
    foreach(p -> markdown(io, p), keycite.prefix)
    print(io, keycite.key)
    foreach(s -> markdown(io, s), keycite.suffix)
end

function markdown(io::IO, cite::Citation)
    foreach(p -> markdown(io, p), cite.globalprefix)
    for keycite in cite.citerefs
        markdown(io, keycite)
        keycite === last(cite.citerefs) || print(io, ", ")
    end
    foreach(s -> markdown(io, s), cite.globalsuffix)
end

markdown(::IO, ::InlineBabelCall) = nothing

markdown(io::IO, src::InlineSourceBlock) =
    print(io, "`", src.body, "`")

markdown(io::IO, ::LineBreak) = print(io, "\n")

markdown(io::IO, link::RadioLink) = markdown(io, link.radio)

function markdown(io::IO, path::LinkPath)
    if path.protocol in keys(link_uri_schemes)
        pathuri = link_uri_schemes[path.protocol](path.path)
        print(io, pathuri)
    end
end

function markdown_image(io::IO, path::LinkPath)
    print(io, "![image](", sprint(markdown, path), ")")
end

function markdown(io::IO, link::Union{PlainLink, AngleLink})
    if splitext(link.path.path)[end] in org_html_image_extensions
        markdown_image(io, link.path)
    else
        print(io, '[', sprint(markdown, link.path), "](",
              string(link.path), ')')
    end
end

function markdown(io::IO, link::RegularLink)
    if isnothing(link.description) || length(link.description) == 0
        if splitext(link.path.path)[end] in org_html_image_extensions
            markdown_image(io, link.path)
        else
            markdown(io, PlainLink(link.path))
        end
    else
        print(io, '[')
        foreach(o -> markdown(io, o), link.description)
        print(io, "](", string(link.path), ")")
    end
end

function markdown(io::IO, o::OrgDoc, mac::Macro)
    expanded = macroexpand(o, mac)
    if isnothing(expanded)
        print(io, "{{{", mac.name,
              '(', join(mac.arguments, ","), ")}}}")
    else
        foreach(o -> markdown(io, o),
                parseorg(expanded, org_object_matchers, org_object_fallbacks))
    end
end

function markdown(io::IO, target::RadioTarget)
    foreach(o -> markdown(io, o), target.contents)
end

function markdown(io::IO, cell::TableCell)
    foreach(o -> markdown(io, o), cell.contents)
end

const markdown_markup_codes =
    Dict(:bold => "**",
         :italic => "*",
         :underline => "__",
         :verbatim => "`",
         :code => "`",
         :strikethrough => "~~")

function markdown(io::IO, markup::TextMarkup)
    print(io, markdown_markup_codes[markup.formatting])
    if markup.contents isa SubString{String}
        print(io, markup.contents)
    else
        markdown.(Ref(io), markup.contents)
    end
    print(io, markdown_markup_codes[markup.formatting])
end

function markdown(io::IO, text::TextPlain)
    tsub = replace(replace(replace(replace(text.text,
                                           "..." => "…"),
                                   r"---([^-])" => s"—\1"),
                           r"--([^-])" => s"—\1"),
                   r"\\-" => "")
    print(io, tsub)
end


text:
################################################################################

@testset verbose = true "Screening" begin

  ################################################################################

  # testing
  @testset "Sensitivity & Specificity" begin
    @test sensitivitySpecificity([10 40; 5 45]) == (sensitivity = 0.6666666666666666, specificity = 0.5294117647058824)
    @test sensitivitySpecificity([20 33; 10 37]) == (sensitivity = 0.6666666666666666, specificity = 0.5285714285714286)
  end

  ################################################################################

  # testing
  @testset "Predicitive Values" begin
    @test predictiveValue([10 40; 5 45]) == (positive = 0.2, negative = 0.9)
    @test predictiveValue([20 33; 10 37]) == (positive = 0.37735849056603776, negative = 0.7872340425531915)
  end

  ################################################################################

end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
module SimpleSolve

import SolverStubs: nlsolve, fixedpoint, seqsolve

export nlsolve, fixedpoint, seqsolve

include("nlsolve.jl")
include("fixedpoint.jl")
include("seqsolve.jl")

end # module


text:
using DryRun
using Test

@testset "DryRun.jl" begin
    include("dryer.jl")
end


text:
export Gridded

struct Gridded{D<:Degree} <: InterpolationType
    degree::D
end

function Base.show(io::IO, g::Gridded)
    print(io, "Gridded(")
    show(io, degree(g))
    print(io, ')')
end

const GridIndex{T} = Union{AbstractVector{T}, Tuple}

# Note: the term "knots" is probably misused, these are more properly "nodes"

struct GriddedInterpolation{T,N,TCoefs,IT<:DimSpec{Gridded},K<:Tuple{Vararg{AbstractVector}}} <: AbstractInterpolation{T,N,IT}
    knots::K
    coefs::Array{TCoefs,N}
    it::IT
end

function Base.:(==)(o1::GriddedInterpolation, o2::GriddedInterpolation)
    o1.it == o2.it &&
    o1.knots == o2.knots &&
    o1.coefs == o2.coefs
end

function GriddedInterpolation(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{TCoefs,N}, it::IT) where {N,TCoefs,TWeights<:Real,IT<:DimSpec{Gridded},pad}
    isconcretetype(IT) || error("The b-spline type must be a leaf type (was $IT)")
    isconcretetype(TCoefs) || @warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")

    check_gridded(it, knots, axes(A))
    c = zero(TWeights)
    if isempty(A)
        T = Base.promote_op(*, typeof(c), eltype(A))
    else
        T = typeof(c * first(A))
    end
    GriddedInterpolation{T,N,TCoefs,IT,typeof(knots)}(knots, A, it)
end

@inline function check_gridded(itpflag, knots, axs)
    flag, ax1, k1 = getfirst(itpflag), axs[1], knots[1]
    if flag isa NoInterp
        k1 == ax1 || error("for NoInterp knot vector should be $ax1, got $k1")
    else
        axes(k1, 1) == ax1 || throw(DimensionMismatch("knot vectors must have the same axes as the corresponding dimension of the array"))
    end
    degree(flag) isa Union{NoInterp,Constant,Linear} || error("only Linear, Constant, and NoInterp supported, got $flag")
    length(k1) == 1 && error("dimensions of length 1 not yet supported")  # FIXME
    issorted(k1) || error("knot-vectors must be sorted in increasing order")
    check_gridded(getrest(itpflag), Base.tail(knots), Base.tail(axs))
end
check_gridded(::Any, ::Tuple{}, ::Tuple{}) = nothing
degree(flag::Gridded) = flag.degree

Base.parent(A::GriddedInterpolation) = A.coefs
coefficients(A::GriddedInterpolation) = A.coefs

size(A::GriddedInterpolation) = size(A.coefs)
axes(A::GriddedInterpolation) = axes(A.coefs)

itpflag(A::GriddedInterpolation) = A.it

function interpolate(::Type{TWeights}, ::Type{TCoefs}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {TWeights,TCoefs,Tel,N,IT<:DimSpec{Gridded}}
    GriddedInterpolation(TWeights, knots, A, it)
end
"""
    itp = interpolate((nodes1, nodes2, ...), A, interpmode)

Interpolate an array `A` on a non-uniform but rectangular grid specified by the given `nodes`,
in the mode determined by `interpmode`.

`interpmode` may be one of

- `NoInterp()`
- `Gridded(Constant())`
- `Gridded(Linear())`

It may also be a tuple of such values, if you want to use different interpolation schemes along each axis.
"""
function interpolate(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {Tel,N,IT<:DimSpec{Gridded}}
    interpolate(tweight(A), tcoef(A), knots, A, it)
end

interpolate!(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {TWeights,Tel,N,IT<:DimSpec{Gridded}} =
    GriddedInterpolation(TWeights, knots, A, it)
function interpolate!(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {Tel,N,IT<:DimSpec{Gridded}}
    interpolate!(tweight(A), tcoef(A), knots, A, it)
end

lbounds(itp::GriddedInterpolation) = first.(itp.knots)
ubounds(itp::GriddedInterpolation) = last.(itp.knots)
lbound(ax::AbstractVector, gr::Gridded)   = lbound(ax, degree(gr))
ubound(ax::AbstractVector, gr::Gridded)   = ubound(ax, degree(gr))

include("indexing.jl")


text:
# Specialized routines for multiplying a subset of the columns of a sparse matrix by a vector or matrix

export colsmul!, tcolsmul!

"""
    colsmul!(C::AbstractMatrix, A::SparseArrays.AbstractSparseMatrixCSC, B::AbstractMatrix, cols, α=1, β=0)

Efficient implementation of the multiplication `α .* A[:, cols] * B[cols, :] .+ β .* C`, where `A`
is sparse. The result is stored in-place in `C`.
"""
function colsmul!(C::AbstractMatrix, A::SparseArrays.AbstractSparseMatrixCSC, B::AbstractMatrix, cols, α=1, β=0)
    length(size(C)) == length(size(B)) || throw(DimensionMismatch("C has dimensions $(size(C)), but B has dimensions $(size(B))"))
    size(C, 1) == size(A, 1) || throw(DimensionMismatch("C has dimensions $(size(C)), but A has dimensions $(size(A))"))
    size(A, 2) == size(B, 1) || throw(DimensionMismatch("A has dimensions $(size(A)), but B has dimensions $(size(B))"))
    0 < minimum(cols) <= maximum(cols) <= size(A, 2) || throw(ArgumentError("A has dimensions $(size(A)), but cols is $cols"))
    C .*= β
    rows = rowvals(A)
    vals = nonzeros(A)
    for col in cols
        for k in 1:size(C, 2)
            @inbounds v = α * B[col, k]
            for i in nzrange(A, col)
                @inbounds row = rows[i]
                @inbounds val = vals[i]
                v1 = val * v
                @inbounds C[row, k] += v1
            end
        end
    end
    C
end

function colsmul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, cols, args...)
    mul!(C, view(A, :, cols), view(B, cols, :), args...)
end

function colsmul!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector, args...)
    colsmul!(reshape(c, length(c), 1), A, reshape(b, length(b), 1), args...)
    c
end

"""

Efficient implementation of the multiplication `α .* Transpose(A[:, cols]) * B .+ β .* C`, where
`A` is sparse. The result is stored in-place in `C[cols, :]`.
"""
function tcolsmul!(C::AbstractMatrix, A::SparseArrays.AbstractSparseMatrixCSC, B::AbstractMatrix, cols, α=1, β=0)
    length(size(C)) == length(size(B)) || throw(DimensionMismatch("C has dimensions $(size(C)), but B has dimensions $(size(B))"))
    size(C, 1) == size(A, 2) || throw(DimensionMismatch("C has dimensions $(size(C)), but A has dimensions $(size(A))"))
    size(A, 1) == size(B, 1) || throw(DimensionMismatch("A has dimensions $(size(A)), but B has dimensions $(size(B))"))
    0 < minimum(cols) <= maximum(cols) <= size(A, 2) || throw(ArgumentError("A has dimensions $(size(A)), but cols is $cols"))
    rows = rowvals(A)
    vals = nonzeros(A)
    C .*= β
    for col in cols
        for k in 1:size(C, 2)
            v1 = zero(eltype(C))
            for i in nzrange(A, col)
                @inbounds row = rows[i]
                @inbounds val = vals[i]
                @inbounds v2 = B[row, k]
                v1 += val*v2
            end
            v1 *= α
            @inbounds C[col, k] += v1
        end
    end
    C
end

function tcolsmul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, cols, args...)
    @views mul!(C[cols, :], transpose(A[:, cols]), B, args...)
end

function tcolsmul!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector, args...)
    tcolsmul!(reshape(c, length(c), 1), A, reshape(b, length(b), 1), args...)
    c
end


text:
# ---
# title: 714. Best Time to Buy and Sell Stock with Transaction Fee
# id: problem714
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/best-time-to-buy-and-sell-stock-with-transaction-fee/description/>
# hidden: true
# ---
#
# Your are given an array of integers `prices`, for which the `i`-th element is
# the price of a given stock on day `i`; and a non-negative integer `fee`
# representing a transaction fee.
#
# You may complete as many transactions as you like, but you need to pay the
# transaction fee for each transaction. You may not buy more than 1 share of a
# stock at a time (ie. you must sell the stock share before you buy again.)
#
# Return the maximum profit you can make.
#
# **Example 1:**
#
#
#
#     Input: prices = [1, 3, 2, 8, 4, 9], fee = 2
#     Output: 8
#     Explanation: The maximum profit can be achieved by:
#
#     * Buying at prices[0] = 1
#     * Selling at prices[3] = 8
#     * Buying at prices[4] = 4
#     * Selling at prices[5] = 9
#     The total profit is ((8 - 1) - 2) + ((9 - 4) - 2) = 8.
#
#
# **Note:**
#
# * `0 < prices.length <= 50000`.
# * `0 < prices[i] < 50000`.
# * `0 <= fee < 50000`.
#
#
## @lc code=start
using LeetCode
function max_profit714(price::Vector{Int}, fee::Int)::Int
    dp1, dp2 = 0, -50000
    for p in price
        dp1, dp2 = max(dp1, dp2 + p - fee), max(dp2, dp1 - p)
    end
    return dp1
end
## @lc code=end


text:
using Pkg
using CSV
using DataFrames
using Distributions
using Random
using StatsBase
using NamedArrays
using Plots
using PyPlot

prices = [20, 25, 18, 20, 15]
score = [0.0, 0, 0, 0]
<<<<<<< HEAD

##Rozkłady
#Ilość samochodów na minutę tworzymy rozkład zeby uzyc sample w funkcji
cars = [1, 2, 3, 4, 5]
prob = [0.15, 0.25, 0.2, 0.2, 0.1]
w = Weights(prob)

##Główna funkcja kroku
#jeden krok = jedna minuta

# function OneStep()
#     #wylosować z rozkładu samochdów na minutę ile samochdów przyjedzie w tym kroku
#     cars_num = sample(cars, w)
#     #wylosować ilość osób w kazdym samochodzie, 1 samochod = jedna zmienna?
#     num_ppl = Dict(i => rand(1:cars_num) for i in 1:cars_num)
#     # wylosować zamówienia dla kazdego samochodu, zapisujemy je jako data DataFrame
#     #gdzie kolumny to produkty zamowione a rzedy to nr samochodu
#     orders = DataFrame(BigMac = zeros(cars_num), Drwal = zeros(cars_num),
#     McNuggets = zeros(cars_num), McChicken = zeros(cars_num), Wege = zeros(cars_num))

#     #losujemy zamowienia dla kazdego samochodu
#     for i in 1:nrow(orders)
#         for j in 1:length(num_ppl)
#             prod = rand(1:ncol(orders))
#             orders[i, prod] += 1
#         end
#     end

function OneStep(start_inv)
    cars_num = sample(cars, w)
=======

function OneStep(start_inv, start, num_cars, weights, start_val, val)
    cars_num = sample(num_cars, weights)
>>>>>>> 0c7659ff753bbe001737e16c6caf863e93042415
    #wylosować ilość osób w kazdym samochodzie, 1 samochod = jedna zmienna?
    num_ppl = Dict(i => rand(1:cars_num) for i in 1:cars_num)
    # wylosować zamówienia dla kazdego samochodu, zapisujemy je jako data DataFrame
    #gdzie kolumny to produkty zamowione a rzedy to nr samochodu
    products = ["BigMac", "Drwal", "McNuggets", "McChicken", "Wege"]
    orders = DataFrame(zeros(Int64, cars_num, length(products)), products)

    #losujemy zamowienia dla kazdego samochodu
    for i in 1:nrow(orders)
        for j in 1:length(num_ppl)
                prod = rand(1:ncol(orders)) # można zmienić na jakiś sensowniejszy rozkład
                orders[i, prod] += 1 # każda osoba wybiera po jednym produkcie?
        end
    end

    #sumujemy zamówienia po kolumnach z wszystkich samochodów
    sum_order = sum.(eachcol(orders))
    #odejmujemy od inventory zamówione składniki
    end_inv = start - sum_order
    #vektor ile zamówień nie udalo sie wypelnic na razie 0
    unfulfilled = zeros(5)

    #jesli inventory zeszlo ponizej zera to tyle ile jest na minusie tyle zamowien nie uda sie spelnic
    for i in 1:length(end_inv)
        if(end_inv[i] < 0)
            unfulfilled[i] -= end_inv[i]
            end_inv[i] = 0
        end
    end
    #spelnione zamowienia to wszystkie zamowienia - niespelnione
    fulfilled = sum_order - unfulfilled
    #przychod
    gain = fulfilled .* prices
    #kara za niespelnienie zamowien
    penalty = unfulfilled .* (0.25*prices)
    #zmniejszamy o jedna minute waznosc produktow
    val .-= 1
    #to co musielismy wyrzucic na razie zera
    thrown_out = zeros(5)
    #jesli waznosc spadnie do 0 to musimy wyrzucic to co zostalo i zrobic nowe oraz zaktualizowac date wwaznosci
    for i in 1:length(val)
        if(val[i] <= 0)
            thrown_out[i] = end_inv[i]
            end_inv[i] += start_inv[i]
            val[i] = start_val[i]
        end
    end
    #strata policzona na podstawie wyrzuconych
    #mozna zmienic te 0.25 i powiedziec ze co sie stanie jak produkty zdrozeja
    loss = thrown_out.* (0.25*prices)
    # zysk
    profit = gain - penalty - loss
    #wyniki
    score[1] = sum(profit) #zysk
    score[2] = sum(thrown_out)#ile wyrzucilismy w tej minucie
    score[3] = abs(sum(unfulfilled)) #ile zamowien zostalo niespelnionych w  tej minucie
    score[4] = sum(fulfilled) # ile zamowien spelnionych w tej minucie

    #nadpisujemy startowe inventory
    start = end_inv
    #zwracamy wynik koncowe inventory oraz obecna date waznosci do nastepnej iteracji
    return score, end_inv, val
end

#horizon to ilosc minut jaka chcemy miec w symulacji najlepiej 60
#val to daty waznosci
#start to inventory na poczatku minuty
#mean_val to ile razy odpalamy symulacje zeby policzyc srednie wartosci
function SimulateOneRun(horizon, val, start, mean_val, prob)
    #Ceny produktów
    prices = [20, 30, 23, 25, 18]
    #pusty wektor do wynikow
    summed_out = zeros(4)
    #tutaj jest petla zeby wyciagac srednia z wynikow, poniewaz losujemy wiele rzeczy
    #chcemy zeby wyniki byly srednia z wielu przebiegow
    #mean_val to ilosc iteracji do liczenia sredniej (im wiecej tym lepiej)
    for i in 1:mean_val
        #obecne daty waznosci
        _val = val
        #stałe poczatkowe daty
        start_val = deepcopy(_val)
        score = [0.0, 0.0, 0.0, 0.0]
        #obecny stan inventory
        _start = start
        #staly poczatkowy stan inventory
        start_inventory = deepcopy(_start)
        #Ilość samochodów na minutę tworzymy rozkład zeby uzyc sample w funkcji
        cars = [1, 2, 3, 4, 5]
        #prob = [0.15, 0.25, 0.2, 0.2, 0.1]
        w = Weights(prob)
        #Wynik symulacji jako dataframe
        df = DataFrame(Profit = 0.0, Thrown_out = 0.0, Unfulfilled = 0.0, Fulfilled = 0.0)
        #symulacja, horizon powinno byc 60 bo symulujemy godzine, mozna zwiekszac
        for i in 1:horizon
            #out to wyniki kroku, _start to stan inventory na koniec kroku/minuty,
            #_val to daty waznosci na koniec kroku/minuty
            out, _start, _val = OneStep(start_inventory, _start, cars, w, start_val, _val)
            #pushujemy do data frame
            push!(df, out[1:4])
        end
        #bierzemy wszystko od drugiego wyniku bo pierwszy to zera,
        #bo tworzylismy data frame z rzedem zer wczesniej
        df_out = df[2:end,:]
        #tutaj musimy skleic rzedami wyniki z kazdej z tych symulacji zeby potem policzyc srednia po rzedach
        summed_out = hcat(summed_out, sum.(eachcol(df_out)))
    end

    #println(start_val)
    #liczymy srednia po rzedach od drugiej do ostatniej kolumny, bo pierwsza to zera
    return mean(summed_out[:, 2:end], dims = 2)
end

#mean_val to ilosc razy jaka odpalimy simulateOneRun zeby policzyc srednia, tam wyzej bylo opisane
function SimulatemanyRuns(mean_val, prob)
    #daty waznosci
    validity = [10, 10, 10, 10, 10]
    #data frame z wynikami
    temp = DataFrame(Profit = 0.0, Thrown_out = 0.0, Unfulfilled = 0.0, Fulfilled = 0.0)
    #przykladowe stany inventory od 1 sztuki do 40 sztuk
    x = [fill(i, 5) for i in 1:40]
    #petla dla kazdego stanu inventory
    for i in x
       push!(temp, SimulateOneRun(60, validity, i, mean_val, prob))
    end
    return temp
end

<<<<<<< HEAD

final = SimulateOneRun(60, [20, 20, 20, 20, 20])

=======
#uruchomienie symulacji
@time xd = SimulatemanyRuns(100, [0.15, 0.25, 0.2, 0.2, 0.1])[2:end, :]
#to wyswietla caly data frame bez ucinania rzedow
show(stdout, "text/plain", xd)
#---------------------------------------------------------
###sprawdzmy dla jakich wartosci jest najwiekszy profit
#wylaczamy kolumne profit
profit = xd[!, "Profit"]
#wykres slupkowy
plt.clf()
plt.bar([1:1:40;], profit, label = "Profit")
#splotujmy jeszcze raz od 10 do 30
profit_2 = profit[10:30]
plt.bar([10:1:30;], profit_2)
profit_3 = profit[18:22]
plt.bar([18:1:22;], profit_3)
xlabel("Number of products")
ylabel("Profit")
gcf()
#---------------------------------------------
#sprawdzmy ile było wyrzucane jedzenia
thrown = xd[!, "Thrown_out"]
plt.clf()
plt.bar([1:1:40;], thrown, color = "red", alpha = 0.7)
xlabel("Number of products")
ylabel("Loss from wasted products")
gcf()

#nakladamy na profit
plt.clf()
plt.bar([1:1:40;], profit, label = "Profit")
#splotujmy jeszcze raz od 10 do 30
profit_2 = profit[10:30]
plt.bar([10:1:30;], profit_2)
profit_3 = profit[18:22]
plt.bar([18:1:22;], profit_3)
plt.bar([1:1:40;], thrown, alpha=0.5)
gcf()

#ilosc spelnionych i niespelnionych zamowien
fulfilled = xd[!, "Fulfilled"]
unfulfilled = xd[!, "Unfulfilled"]
plt.clf()
plt.bar([1:1:40;], fulfilled, color = "green")
plt.bar([1:1:40;], unfulfilled, color = "red")
legend(["Fulfilled", "Unfulfilled"])
xlabel("Number of products")
ylabel("Orders")
gcf()

#zmiana ruchu samochodów na bardziej intensywny i mniej intensywny
@time xd_2 = SimulatemanyRuns(100, [0.1, 0.2, 0.25, 0.25, 0.2])[2:end, :]
@time xd_3 = SimulatemanyRuns(100, [0.25, 0.3, 0.2, 0.15, 0.1])[2:end, :]
profit = xd[!, "Profit"]
profit_22 = xd_2[!, "Profit"]
profit_33 = xd_3[!, "Profit"]
#wykres slupkowy
plt.clf()

plt.bar([1:1:40;], profit_22, alpha = 0.5, color = "red")
plt.bar([1:1:40;], profit, color = "green")
plt.bar([1:1:40;], profit_33, alpha = 0.8, color = "blue")
xlabel("Number of products")
ylabel("Profit")
legend(["Low Traffic", "High Traffic", "Medium traffic"])
gcf()
>>>>>>> 0c7659ff753bbe001737e16c6caf863e93042415

max_stars_repo_name: myrchas/food-orders-simulation

max_issues_repo_name: myrchas/food-orders-simulation

max_forks_repo_name: myrchas/food-orders-simulation


text:
export initFarfieldWav
export farfieldWav

"""
    y = farfieldWav(wav::Array, fs::Real=16000.0; maxvalue::Real=0.6, minstep::Real=-0.6)
simulate far-field effects by auto-gain-ctrl method.
"""
function farfieldWav(wav::Array, fs::Real=16000.0; maxvalue::Real=0.6, minstep::Real=-0.6)
    MAX = maximum(abs.(wav)) + 3e-5
    wav .*= 1.0 / MAX;
    minstep  = (-1.0 < minstep < 0.0) ? minstep : -0.6
    maxvalue = (0.0 < maxvalue < 1.0) ? maxvalue : 0.6
    gain   = 1.0
    winlen = floor(Int, 0.016 * fs)
    frames = div(length(wav)-winlen, winlen) + 1
    firstIds = (0:(frames-1)) .* winlen .+ 1       # 帧起始下标
    lasstIds = firstIds .+ (winlen - 1)            # 帧结束下标
    maxvalue = 1 / maxvalue
    for t = 1:frames
        index = firstIds[t]:lasstIds[t];
        frame = wav[index];
        fmax = maximum(abs.(frame));
        step = 1.0 - fmax*maxvalue;
        step = step * abs(step);
        step = max(minstep, min(step, 0.0));
        gain = 0.8*gain + 0.2*(1 + step);
        wav[index] .*= gain * MAX;
    end
    return wav
end

"""
    initFarfieldWav(fs::Real, maxvalueSpan::NTuple{2,Number}) -> farfieldwav(wav::Array)
init farfield effect function
+ `fs` sampling rate e.g. 16000.0
+ `maxvalueSpan` e.g. (0.2, 1.0)
"""
function initFarfieldWav(fs::Real, maxvalueSpan::NTuple{2,Number})
    MinV, MaxV = maxvalueSpan
    @assert MinV <= MaxV
    function farfieldwav(wav::Array)
        v = rand()*(MaxV - MinV) + MinV
        return farfieldWav(wav, fs, maxvalue=v)
    end
    return farfieldwav
end

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
module TestDay14

import AOC2021.Day14.part1, AOC2021.Day14.part2
using Test

test_input = """
NNCB

CH -> B
HH -> N
CB -> H
NH -> C
HB -> C
HC -> B
HN -> C
NN -> C
BH -> H
NC -> B
NB -> B
BN -> B
BB -> N
BC -> B
CC -> N
CN -> C"""

@testset "Day 14" begin
    @testset "part 1" begin
        @test part1(test_input) == 1588
    end
    @testset "part 2" begin
        @test part2(test_input) == 2188189693529
    end
end

end # module


text:
using Luxor

include("randgen.jl")

const RES_W = 64
const RES_H = 64

@enum DistributionMode begin
    EasyMode = 0
    HardMode = 1
    ExtremeMode = 2
    MemoryMode = 10
end

mutable struct StepData
    reward::Float32
    done::Bool
    level_complete::Bool

    function StepData()

        reward = 0.0
        done = false
        level_complete = false

        new(
            reward,
            done,
            level_complete,
        )
    end
end

mutable struct GameOptions
    use_monochrome_assets::Bool
    restrict_themes::Bool
    use_backgrounds::Bool
    center_agent::Bool
    distribution_mode::DistributionMode
    use_sequential_levels::Bool

    function GameOptions()

        use_monochrome_assets = false
        restrict_themes = false
        use_backgrounds = true
        center_agent = false
        distribution_mode = HardMode
        use_sequential_levels = false

        new(
            use_monochrome_assets,
            restrict_themes,
            use_backgrounds,
            center_agent,
            distribution_mode,
            use_sequential_levels,
        )
    end
end

abstract type AbstractGame end

mutable struct Game <: AbstractGame
    game_name::String

    options::GameOptions

    grid_step::Bool
    level_seed_low::Int
    level_seed_high::Int

    level_seed_rand_gen::RandGen
    rand_gen::RandGen

    step_data::StepData
    action::Int

    timeout::Int

    current_level_seed::Int
    prev_level_seed::Int
    episodes_remaining::Int
    episode_done::Bool

    last_reward_timer::Int
    last_reward::Float32
    default_action::Int

    fixed_asset_seed::Int

    cur_time::Int

    reset_count::Int
    total_reward::Float32

    function Game(name::String)

        game_name = name

        options = GameOptions()

        grid_step = false
        level_seed_low = 0
        level_seed_high = 1

        level_seed_rand_gen = RandGen()
        rand_gen = RandGen()

        step_data = StepData()
        step_data.reward = 0.0
        step_data.done = true
        step_data.level_complete = false
        action = 0

        timeout = 1000

        current_level_seed = 0
        prev_level_seed = 0
        episodes_remaining = 0
        episode_done = false

        last_reward_timer = 0
        last_reward = -1.0
        default_action = 0

        fixed_asset_seed = 0

        cur_time = 0

        reset_count = 0
        total_reward = 0.0

        new(
            game_name,
            options,
            grid_step,
            level_seed_low,
            level_seed_high,
            level_seed_rand_gen,
            rand_gen,
            step_data,
            action,
            timeout,
            current_level_seed,
            prev_level_seed,
            episodes_remaining,
            episode_done,
            last_reward_timer,
            last_reward,
            default_action,
            fixed_asset_seed,
            cur_time,
            reset_count,
            total_reward,
        )
    end
end

function get_game(abstract_game::AbstractGame)

    return abstract_game
end

function render_to_buf(abstract_game::AbstractGame, w::Int, h::Int)

    rectangle = rect(Point(0, 0), w, h)
    game_draw(abstract_game, rectangle)
    return nothing
end

function reset(abstract_game::AbstractGame)

    game = get_game(abstract_game)

    game.reset_count += 1

    if game.episodes_remaining == 0
        if game.options.use_sequential_levels && game.step_data.level_complete
            game.current_level_seed = game.current_level_seed + 997
        else
            game.current_level_seed = randint(game.level_seed_rand_gen, game.level_seed_low, game.level_seed_high)
        end

        game.episodes_remaining = 1
    else
        game.step_data.reward = 0
        game.step_data.done = false
        game.step_data.level_complete = false
    end

    seed(game.rand_gen, game.current_level_seed)
    game_reset(abstract_game)

    game.cur_time = 0
    game.total_reward = 0
    game.episodes_remaining -= 1
    game.action = game.default_action
    return nothing
end

function step(abstract_game::AbstractGame)

    game = get_game(abstract_game)

    game.cur_time += 1
    will_force_reset = false

    if game.action == -1
        game.action = game.default_action
        will_force_reset = true
    end

    game.step_data.reward = 0
    game.step_data.done = false
    game.step_data.level_complete = false
    game_step(abstract_game)

    game.step_data.done = game.step_data.done || will_force_reset || (game.cur_time >= game.timeout)
    game.total_reward += game.step_data.reward

    if game.step_data.reward != 0
        game.last_reward_timer = 10
        game.last_reward = game.step_data.reward
    end

    game.prev_level_seed = game.current_level_seed

    if game.step_data.done
        reset(abstract_game)
    end

    if game.options.use_sequential_levels && game.step_data.level_complete
        game.step_data.done = false
    end

    game.episode_done = game.step_data.done

    observe(abstract_game)
    return nothing
end

function observe(abstract_game::AbstractGame)

    render_to_buf(abstract_game, RES_W, RES_H)
    return nothing
end

function game_init(abstract_game::AbstractGame)

    return nothing
end

function game_reset(abstract_game::AbstractGame)

    return nothing
end

function game_step(abstract_game::AbstractGame)

    return nothing
end

function game_draw(abstract_game::AbstractGame, rectangle)

    return nothing
end


text:
#
# Spline
#

struct Spline
  x :: Vector{Float64}
  y :: Vector{Float64}
end

"""

`fitspline(x,y)`

Computes a spline.

Use `fit = fitspline(x,y)`. `fit.x` and `fit.y` will contain
the data points of the spline.

# Example

```jldoctest
julia> x = sort(rand(10)); y = rand(10);

julia> fit = fitspline(x,y)

 -------- Spline fit ---------------------------

 x spline: x = [0.05986115284955605, 0.0755869317846746...
 y spline: y = [0.6749928175217872, 0.5711266179726132...

 -----------------------------------------------
```
"""
function fitspline(X :: AbstractArray{<:Real}, Y :: AbstractArray{<:Real}, options :: Options)
  X, Y = checkdata(X,Y,options)
  t = 1:length(X)
  A = hcat(X,Y)
  itp = Interpolations.scale(interpolate(A,
            (BSpline(Interpolations.Cubic(Natural(OnGrid()))), NoInterp())), t, 1:2)
  tfine = 1:length(X)/options.fine:length(X)
  x, y = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
  return Spline(x,y)
end
fitspline(X :: AbstractArray{<:Real}, Y :: AbstractArray{<:Real}) = fitspline(X,Y,Options())

function Base.show( io :: IO, fit :: Spline )
  println("")
  println(" -------- Spline fit --------------------------- ")
  println("")
  println(" x spline: x = [",fit.x[1],", ",fit.x[2],"...")
  println(" y spline: y = [",fit.y[1],", ",fit.y[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

export fitspline


text:
using Plots
using NBodySimulator
using StaticArrays

struct CustomPotentialParameters <: PotentialParameters
    a::AbstractFloat
end

import NBodySimulator.get_accelerating_function
function get_accelerating_function(p::CustomPotentialParameters, simulation::NBodySimulation)
    ms = get_masses(simulation.system)
    (dv, u, v, t, i) -> begin custom_accel = SVector(0.0, 0.0, p.a); dv .= custom_accel*ms[i] end
end

parameters = CustomPotentialParameters(-9.8)
system = PotentialNBodySystem(bodies, Dict(:custom_potential_params => parameters))

function MakeParticles(N, m, q, size, offset, velocity)
    coords = size * rand(3, N) .+ offset
    speeds = velocity * rand(3, N)

    particles = Array{ChargedParticle{cType, Float64, Float64} where cType<:Real, 1}()

    for i in 1:N
        push!(particles, ChargedParticle(SVector(coords[1, i], coords[2, i], coords[3, i]), SVector(speeds[1, i], speeds[2, i], speeds[3, i]), m, q))
    end

    return particles
end

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function denman_beavers(A)
    return denman_beavers!(copy(A))
end

@inline function denman_beavers!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    max_memslots=4
    memslots=Vector{Matrix{T}}(undef,max_memslots)
    n=size(A,1)
    for j=1:max_memslots
        memslots[j]=Matrix{T}(undef,n,n)
    end
    # The first slots are precomputed nodes [:A]
    memslots[1]=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A_shift Xinv0 Y1 X1 Yinv1 Xinv1 Y2 Yinv2 X3
    # Computing A_shift = x*A+x*I
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle A
    memslots[1] .= coeff1.*memslots[1]
    mul!(memslots[1],true,I*coeff2,true,true)
    # Computing Xinv0 with operation: ldiv
    memslots[2]=inv(memslots[1])
    # Computing Y1 = x*I+x*Xinv0
    coeff1=0.5
    coeff2=0.5
    # Smart lincomb recycle Xinv0
    memslots[2] .= coeff2.*memslots[2]
    mul!(memslots[2],true,I*coeff1,true,true)
    # Computing X1 = x*A_shift+x*I
    coeff1=0.5
    coeff2=0.5
    # Smart lincomb recycle A_shift
    memslots[1] .= coeff1.*memslots[1]
    mul!(memslots[1],true,I*coeff2,true,true)
    # Computing Yinv1 with operation: ldiv
    memslots[3]=inv(memslots[2])
    # Computing Xinv1 with operation: ldiv
    memslots[4]=inv(memslots[1])
    # Computing Y2 = x*Y1+x*Xinv1
    coeff1=0.5
    coeff2=0.5
    # Smart lincomb recycle Y1
    memslots[2] .= coeff1.*memslots[2] .+ coeff2.*memslots[4]
    # Deallocating Xinv1 in slot 4
    # Computing Yinv2 with operation: ldiv
    memslots[4]=inv(memslots[2])
    # Deallocating Y2 in slot 2
    # Computing X3 = x*X1+x*Yinv1+x*Yinv2
    coeff1=0.25
    coeff2=0.25
    coeff3=0.5
    # Smart lincomb recycle X1
    memslots[1] .= coeff1.*memslots[1] .+ coeff2.*memslots[3] .+ coeff3.*memslots[4]
    # Deallocating Yinv1 in slot 3
    # Deallocating Yinv2 in slot 4
    return memslots[1] # Returning X3
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


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

mutable struct BootDiagnostics <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    storageUri::Any # spec type: Union{ Nothing, String } # spec name: storageUri

    function BootDiagnostics(;enabled=nothing, storageUri=nothing)
        o = new()
        validate_property(BootDiagnostics, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(BootDiagnostics, Symbol("storageUri"), storageUri)
        setfield!(o, Symbol("storageUri"), storageUri)
        o
    end
end # type BootDiagnostics

const _property_map_BootDiagnostics = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("storageUri")=>Symbol("storageUri"))
const _property_types_BootDiagnostics = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("storageUri")=>"String")
Base.propertynames(::Type{ BootDiagnostics }) = collect(keys(_property_map_BootDiagnostics))
Swagger.property_type(::Type{ BootDiagnostics }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BootDiagnostics[name]))}
Swagger.field_name(::Type{ BootDiagnostics }, property_name::Symbol) =  _property_map_BootDiagnostics[property_name]

function check_required(o::BootDiagnostics)
    true
end

function validate_property(::Type{ BootDiagnostics }, name::Symbol, val)
end


text:
__precompile__()

module Integrators

    using DelimitedFiles
    using Documenter
    using LinearAlgebra
    using OffsetArrays

    using Base: TwicePrecision

    using ..CommonFunctions
    using ..Config
    using ..Interpolation
    using ..Utils

    using ..BasisFunctions
    using ..Quadratures
    using ..Discontinuities
    using ..Equations
    using ..Solutions
    using ..Solvers

    export InitialGuessODE, InitialGuessPODE, initialize!, update!

    include("integrators/initial_guess/extrapolation.jl")
    include("integrators/initial_guess/initial_guess_ode.jl")
    include("integrators/initial_guess/initial_guess_pode.jl")

    export AbstractCoefficients, AbstractTableau

    include("integrators/abstract_coefficients.jl")
    include("integrators/abstract_tableau.jl")

    export Integrator, DeterministicIntegrator, StochasticIntegrator, IntegratorCache
    export integrate, integrate!, equation, timestep
    export NonlinearFunctionParameters, function_stages!

    include("integrators/abstract_integrator.jl")
    include("integrators/integrator_cache.jl")

    include("integrators/integrators_common.jl")

    export CoefficientsRK, AbstractTableauRK, AbstractTableauIRK, AbstractTableauPRK,
           IntegratorRK, writeTableauToFile

    export get_symplectic_conjugate_coefficients, symplecticize,
           check_symplecticity, check_symmetry, compute_symplecticity_error,
           check_order_conditions_B, check_order_conditions_C, check_order_conditions_D

    include("integrators/rk/abstract_integrator_rk.jl")
    include("integrators/rk/coefficients.jl")
    include("integrators/rk/tableaus.jl")

    export IntegratorERK, TableauERK, readTableauERKFromFile

    include("integrators/rk/integrators_erk.jl")

    export IntegratorDIRK, TableauDIRK
    export IntegratorFIRK, TableauFIRK

    include("integrators/rk/integrators_dirk.jl")
    include("integrators/rk/integrators_firk.jl")

    export IntegratorEPRK, TableauEPRK
    export IntegratorIPRK, TableauIPRK
    export IntegratorFLRK

    include("integrators/rk/integrators_eprk.jl")
    include("integrators/rk/integrators_iprk.jl")
    include("integrators/rk/integrators_flrk.jl")

    export IntegratorPGLRK, CoefficientsPGLRK

    include("integrators/rk/integrators_pglrk.jl")

    export IntegratorSERK, TableauSERK
    export IntegratorSIRK, TableauSIRK
    export IntegratorSIPRK, TableauSIPRK
    export IntegratorSISPRK, TableauSISPRK
    export IntegratorWERK, TableauWERK
    export IntegratorWIRK, TableauWIRK

    include("integrators/stochastic_integrators_common.jl")
    include("integrators/stochastic/integrators_serk.jl")
    include("integrators/stochastic/integrators_sirk.jl")
    include("integrators/stochastic/integrators_siprk.jl")
    include("integrators/stochastic/integrators_sisprk.jl")
    include("integrators/stochastic/integrators_werk.jl")
    include("integrators/stochastic/integrators_wirk.jl")

    export IntegratorVPRK, IntegratorVPRKpNone, TableauVPRK

    export IntegratorVPRKdegenerate

    export IntegratorVPRKpStandard, IntegratorVPRKpSymplectic,
           IntegratorVPRKpInternal, IntegratorVPRKpMidpoint,
           IntegratorVPRKpSymmetric,
           IntegratorVPRKpSecondary, IntegratorVPRKpVariational,
           IntegratorVPRKpVariationalQ, IntegratorVPRKpVariationalP

    export IntegratorVPRKpLegendre, TableauVPRKpLegendre

    include("integrators/vprk/integrators_vprk_abstract.jl")
    include("integrators/vprk/integrators_vprk_cache.jl")
    include("integrators/vprk/integrators_vprk_common.jl")
    include("integrators/vprk/integrators_vprk.jl")
    include("integrators/vprk/integrators_vprk_degenerate.jl")
    include("integrators/vprk/integrators_vprk_pinternal.jl")
    include("integrators/vprk/integrators_vprk_pmidpoint.jl")
    include("integrators/vprk/integrators_vprk_pstandard.jl")
    include("integrators/vprk/integrators_vprk_psecondary.jl")
    include("integrators/vprk/integrators_vprk_psymmetric.jl")
    include("integrators/vprk/integrators_vprk_pvariational.jl")
    include("integrators/vprk/integrators_vprk_plegendre.jl")

    export CoefficientsARK, CoefficientsPRK, CoefficientsMRK, CoefficientsIRK,
           CoefficientsSPARK

    export AbstractIntegratorSPARK
    export AbstractTableauSPARK, TableauSPARK

    export IntegratorHPARK, TableauHPARK
    export IntegratorVPARK, TableauVPARK

    export IntegratorHSPARK, TableauHSPARK
    export IntegratorHSPARKprimary, TableauHSPARKprimary
    export IntegratorHSPARKsecondary, TableauHSPARKsecondary

    export IntegratorVSPARK, TableauVSPARK
    export IntegratorVSPARKprimary, TableauVSPARKprimary
    export IntegratorVSPARKsecondary, TableauVSPARKsecondary

    include("integrators/spark/abstract_integrator_spark.jl")
    include("integrators/spark/coefficients.jl")

    include("integrators/spark/integrators_spark_cache.jl")
    include("integrators/spark/integrators_spark_common.jl")
    include("integrators/spark/integrators_spark_tableau.jl")

    include("integrators/spark/integrators_vpark.jl")
    include("integrators/spark/integrators_hpark.jl")

    include("integrators/spark/integrators_vspark_common.jl")
    include("integrators/spark/integrators_vspark.jl")
    include("integrators/spark/integrators_vspark_primary.jl")
    include("integrators/spark/integrators_vspark_secondary.jl")

    include("integrators/spark/integrators_hspark_common.jl")
    include("integrators/spark/integrators_hspark.jl")
    include("integrators/spark/integrators_hspark_primary.jl")
    include("integrators/spark/integrators_hspark_secondary.jl")

    export TableauGLM

    include("integrators/glm/integrators_glm.jl")

    export IntegratorSplitting, AbstractTableauSplitting,
           TableauSplittingGS, TableauSplittingNS, TableauSplittingSS

    include("integrators/splitting/integrators_splitting.jl")

    export IntegratorCGVI, IntegratorDGVI, IntegratorDGVIEXP,
           IntegratorDGVIPI, IntegratorDGVIP0, IntegratorDGVIP1

    include("integrators/cgvi/integrators_cgvi.jl")
    include("integrators/dgvi/integrators_dgvi.jl")
    include("integrators/dgvi/integrators_dgvi_experimental.jl")
    include("integrators/dgvi/integrators_dgvi_path_integral.jl")
    include("integrators/dgvi/integrators_dgvi_projection_initial.jl")
    include("integrators/dgvi/integrators_dgvi_projection_final.jl")

    include("integrators/integrators.jl")

    function __init__()
        default_params = (
            (:ig_interpolation, HermiteInterpolation),
            (:ig_extrapolation_stages, 5),
            (:int_show_progress_nmin,  1000),
        )

        for param in default_params
            add_config(param...)
        end
    end

end


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

name = "SLATEC"
version = v"4.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sabjohnso/slatec.git",
              "417db9e31c49eba4aee5ab9bb719093f6886bcee"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/slatec/src

# Manually set LAPACK and BLAS libraries
echo "list( REMOVE_ITEM SLATEC_SOURCE_FILES \${BLAS_SOURCE_FILES} )" >> CMakeLists.txt
echo "list( APPEND SLATEC_SOURCE_FILES \${RETAINED_BLAS_SOURCE_FILES} )" >> CMakeLists.txt
echo "list( REMOVE_ITEM SLATEC_SOURCE_FILES \${EXCLUDE_FROM_SHARED} )" >> CMakeLists.txt
echo "include(\$ENV{prefix}/lib/cmake/lapack-3.9.0/lapack-config.cmake)" >> CMakeLists.txt
echo "include(\$ENV{prefix}/lib/cmake/openblas/OpenBLASConfig.cmake)" >> CMakeLists.txt
echo "add_library( slatec_shared SHARED \${SLATEC_SOURCE_FILES} )" >> CMakeLists.txt
echo "target_link_libraries( slatec_shared \${LAPACK_LIBRARIES} \${OpenBLAS_LIBRARIES})" >> CMakeLists.txt
echo "set_target_properties( slatec_shared PROPERTIES OUTPUT_NAME slatec )" >> CMakeLists.txt
echo "install( TARGETS slatec_shared LIBRARY DESTINATION lib )" >> CMakeLists.txt

if [[ "${target}" == aarch64-apple-* ]]; then
    # Fix issue due to GCC 10+.
    #     [  2%] Building Fortran object src/CMakeFiles/slatec_shared.dir/dwnlsm.f.o
    #     cd /workspace/srcdir/slatec/build/src && /opt/bin/aarch64-apple-darwin20-libgfortran5-cxx11/aarch64-apple-darwin20-gfortran --sysroot=/opt/aarch64-apple-darwin20/aarch64-apple-darwin20/sys-root -Dslatec_shared_EXPORTS  -O3 -DNDEBUG -O3 -fPIC   -c /workspace/srcdir/slatec/src/dwnlsm.f -o CMakeFiles/slatec_shared.dir/dwnlsm.f.o
    #     /workspace/srcdir/slatec/src/dwnlsm.f:440:28:
    #
    #       119 |       CALL DCOPY (N, 1.D0, 0, D, 1)
    #           |                     2
    #     ......
    #       440 |          CALL DCOPY (NSOLN, Z, 1, X, 1)
    #           |                            1
    #     Error: Rank mismatch between actual argument at (1) and actual argument at (2) (scalar and rank-1)
    #     /workspace/srcdir/slatec/src/dwnlsm.f:608:25:
    #
    #       119 |       CALL DCOPY (N, 1.D0, 0, D, 1)
    #           |                     2
    #     ......
    #       608 |       CALL DCOPY (NSOLN, Z, 1, X, 1)
    #           |                         1
    #     Error: Rank mismatch between actual argument at (1) and actual argument at (2) (scalar and rank-1)
    export FFLAGS="-fallow-argument-mismatch"
fi

mkdir ../build
cd ../build/
# Above on the fly added CMake code buils shared library with OpenBLAS
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DUSE_OPTIMIZED_BLAS=0 \
    -DBUILD_SHARED_LIBRARY=0 \
    ..
make -j${nproc}
make install

install_license ../../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libslatec", :libslatec)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="LAPACK_jll", uuid="51474c39-65e3-53ba-86ba-03b1b862ec14")),
    Dependency(PackageSpec(name="OpenBLAS_jll", uuid="4536629a-c528-5b80-bd46-f80d51c5b363")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
# For the time being need LLVM 11 because of <https://github.com/JuliaPackaging/BinaryBuilderBase.jl/issues/158>.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_llvm_version=v"11")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ------------------------------------------------------------------------------------------
# # Building a single neural network layer using `Flux.jl`
#
# In this notebook, we'll move beyond binary classification. We'll try to distinguish
# between three fruits now, instead of two. We'll do this using **multiple** neurons
# arranged in a **single layer**.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Read in and process data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can start by loading the necessary packages and getting our data into working order
# with similar code we used at the beginning of the previous notebooks, except that now we
# will combine three different apple data sets, and will add in some grapes to the fruit
# salad!
# ------------------------------------------------------------------------------------------

# Load apple data in with `readdlm` for each file
apples1, applecolnames1 = readdlm("data/Apple_Golden_1.dat", '\t', header = true)
apples2, applecolnames2 = readdlm("data/Apple_Golden_2.dat", '\t', header = true)
apples3, applecolnames3 = readdlm("data/Apple_Golden_3.dat", '\t', header = true)

# Check that the column names are the same for each apple file
println( applecolnames1 == applecolnames2 == applecolnames3)

# ------------------------------------------------------------------------------------------
# Since each apple file has columns with the same headers, we know we can concatenate these
# columns from the different files together:
# ------------------------------------------------------------------------------------------

apples = vcat(apples1, apples2, apples3)

# ------------------------------------------------------------------------------------------
# And now let's build an array called `x_apples` that stores data from the `red` and `blue`
# columns of `apples`. From `applecolnames1`, we can see that these are the 3rd and 5th
# columns of `apples`:
# ------------------------------------------------------------------------------------------

applecolnames1[3], applecolnames1[5]

length(apples[:, 1])

x_apples  = [ [apples[i, 3], apples[i, 5]] for i in 1:length(apples[:, 3]) ]

# ------------------------------------------------------------------------------------------
# Similarly, let's create arrays called `x_bananas` and `x_grapes`:
# ------------------------------------------------------------------------------------------

# Load data from *.dat files
bananas, bananacolnames = readdlm("data/Banana.dat", '\t', header = true)
grapes1, grapecolnames1 = readdlm("data/Grape_White.dat", '\t', header = true)
grapes2, grapecolnames2 = readdlm("data/Grape_White_2.dat", '\t', header = true)

# Concatenate data from two grape files together
grapes = vcat(grapes1, grapes2)

# Check that column 3 and column 5 refer to the "red" and "blue" columns from each file
println("All column headers are the same: ", bananacolnames == grapecolnames1 == grapecolnames2 == applecolnames1)

# Build x_bananas and x_grapes from bananas and grapes
x_bananas  = [ [bananas[i, 3], bananas[i, 5]] for i in 1:length(bananas[:, 3]) ]
x_grapes = [ [grapes[i, 3], grapes[i, 5]] for i in 1:length(grapes[:, 3]) ]

# ------------------------------------------------------------------------------------------
# ## One-hot vectors
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now we wish to classify *three* different types of fruit. It is not clear how to encode
# these three types using a single output variable; indeed, in general this is not possible.
#
# Instead, we have the idea of encoding $n$ output types from the classification into
# *vectors of length $n$*, called "one-hot vectors":
#
# $$
# \textrm{apple} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix};
# \quad
# \textrm{banana} = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix};
# \quad
# \textrm{grape} = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}.
# $$
#
# The term "one-hot" refers to the fact that each vector has a single $1$, and is $0$
# otherwise.
#
# Effectively, the first neuron will learn whether or not (1 or 0) the data corresponds to
# an apple, the second whether or not (1 or 0) it corresponds to a banana, etc.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `Flux.jl` provides an efficient representation for one-hot vectors, using advanced
# features of Julia so that it does not actually store these vectors, which would be a waste
# of memory; instead `Flux` just records in which position the non-zero element is. To us,
# however, it looks like all the information is being stored:
# ------------------------------------------------------------------------------------------

using Flux: onehot

onehot(1, 1:3)

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Make an array `labels` that gives the labels (1, 2 or 3) of each data point. Then use
# `onehot` to encode the information about the labels as a vector of `OneHotVector`s.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Single layer in Flux
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's suppose that there are two pieces of input data, as in the previous single neuron
# notebook. Then the network has 2 inputs and 3 outputs:
# ------------------------------------------------------------------------------------------

include("draw_neural_net.jl")
draw_network([2, 3])

# ------------------------------------------------------------------------------------------
# `Flux` allows us to express this again in a simple way:
# ------------------------------------------------------------------------------------------

using Flux

model = Dense(2, 3, σ)

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Now what do the weights inside `model` look like? How does this compare to the diagram of
# the network layer above?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Training the model
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Despite the fact that the model is now more complicated than the single neuron from the
# previous notebook, the beauty of `Flux.jl` is that the rest of the training process
# **looks exactly the same**!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Implement training for this model.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
#
# Visualize the result of the learning for each neuron. Since each neuron is sigmoidal, we
# can get a good idea of the function by just plotting a single contour level where the
# function takes the value 0.5, using the `contour` function with keyword argument
# `levels=[0.5, 0.501]`.
# ------------------------------------------------------------------------------------------

plot()

contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[1], levels=[0.5, 0.501], color = cgrad([:blue, :blue]))
contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[2], levels=[0.5,0.501], color = cgrad([:green, :green]))
contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[3], levels=[0.5,0.501], color = cgrad([:red, :red]))

scatter!(first.(x_apples), last.(x_apples), m=:cross, label="apples")
scatter!(first.(x_bananas), last.(x_bananas), m=:circle, label="bananas")
scatter!(first.(x_grapes), last.(x_grapes), m=:square, label="grapes")

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# Interpret the results by checking which fruit each neuron was supposed to learn and what
# it managed to achieve.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/16. Tools - Using Flux to build a single layer neural net.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
addprocs(3)

trials = 4000000
hits = @parallel (+) for i = 1:trials
  Int(sum(rand(2).^2) < 1)
end

approximation = hits/trials * 4
@assert(3.1 < approximation < 3.2)


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

const elementary_charge = Float64(1.602176487e-19)
const ϵ0  = Float64(8.8541878176e-12)

const material_properties = Dict{Symbol, NamedTuple}()

material_properties[:Vacuum] = (
    E_ionisation = 0.0u"eV",
    f_fano = 0.0,
    ϵ_r = 1.0,
    ρ = 0Unitful.g / (Unitful.cm^3),
    name = "Vacuum"
)

material_properties[:HPGe] = (
    E_ionisation = 2.95u"eV",
    f_fano = 0.13,
    ϵ_r = 16.0,
    ρ = 5.323Unitful.g / (Unitful.cm^3), # u"g/cm^3" throws warnings in precompilation
    name = "High Purity Germanium"
)

# These values might just be approximationsl
material_properties[:Si] = (
    E_ionisation = 3.62u"eV",
    f_fano = 0.11,
    ϵ_r = 11.7,
    ρ = 2.3290Unitful.g / (Unitful.cm^3), # u"g/cm^3" throws warnings in precompilation
    name = "Silicon"
)
material_properties[:Al] = (
    name = "Aluminium",
    ϵ_r = 10.8, # Aluminium Foil
    ρ = 2.6989Unitful.g / (Unitful.cm^3) # u"g/cm^3" throws warnings in precompilation
)

material_properties[:LAr] = (
    name = "liquid Argon",
    E_ionisation = 0u"eV",
    f_fano = 0.107,
    ϵ_r = 1.505,
    ρ = 1.396Unitful.g / (Unitful.ml) # u"g/cm^3" throws warnings in precompilation
)

material_properties[:Co] = (
    name = "Copper",
    ϵ_r = 20,
    ρ = 8.96Unitful.g / (Unitful.ml) # u"g/cm^3" throws warnings in precompilation
)

material_properties[:CdZnTe] = (
    name = "Cadmium zinc telluride",
    E_ionisation = 4.64u"eV",
    ϵ_r = 10.9,
    ρ = 5.78Unitful.g / (Unitful.ml) # u"g/cm^3" throws warnings in precompilation
)
# Add new materials above this line
# and just put different spellings into the dict `materials` below

materials = Dict{String, Symbol}(
    "HPGe" => :HPGe,
    "vacuum" => :Vacuum,
    "Vacuum" => :Vacuum,
    "Copper" => :Co,
    "copper" => :Co,
    "Al"  => :Al,
    "CZT" => :CdZnTe,
    "Si" => :Si
)

for key in keys(material_properties)
    if !haskey(materials, string(key))
        push!(materials, string(key) => key)
    end
end


text:
"""
    HistoryFunction(h, integrator)

Wrap history function `h` and integrator `integrator` to create a common
interface for retrieving values at any time point with varying accuracy.

Before the initial time point of the solution of the `integrator` values are calculated by
history function `h`, for time points up to the final time point of the solution
interpolated values of the solution are returned, and after the final time point an inter-
or extrapolation of the current state of the `integrator` is retrieved.
"""
mutable struct HistoryFunction{H,I<:DEIntegrator} <: DiffEqBase.AbstractHistoryFunction
  h::H
  integrator::I
  isout::Bool
end

HistoryFunction(h, integrator) = HistoryFunction(h, integrator, false)

function (f::HistoryFunction)(p, t, ::Type{Val{deriv}}=Val{0}; idxs=nothing) where deriv
  @unpack integrator = f
  @unpack tdir, sol = integrator

  tdir_t = tdir * t

  @inbounds if tdir_t < tdir * sol.t[1]
    if deriv == 0 && idxs === nothing
      return f.h(p, t)
    elseif idxs === nothing
      return f.h(p, t, Val{deriv})
    elseif deriv == 0
      return f.h(p, t; idxs = idxs)
    else
      return f.h(p, t, Val{deriv}; idxs = idxs)
    end
  end

  tdir_solt = tdir * sol.t[end]
  if tdir_t <= tdir_solt
    return sol.interp(t, idxs, Val{deriv}, p)
  end

  # history function is evaluated at time point past the final time point of
  # the current solution
  if tdir_t - tdir_solt > 10 * eps(tdir_t)
    f.isout = true
  end

  if integrator.t == sol.prob.tspan[1]
    # handle extrapolations at initial time point
    return constant_extrapolant(t, integrator, idxs, Val{deriv})
    #TODO only difference form DelayDiffEq here
  else
    return StochasticDiffEq.current_interpolant(t, integrator, idxs, Val{deriv})
    #TODO only difference form DelayDiffEq here
  end
end

function (f::HistoryFunction)(val, p, t, ::Type{Val{deriv}}=Val{0}; idxs=nothing) where deriv
  @unpack integrator = f
  @unpack tdir, sol = integrator

  tdir_t = tdir * t

  @inbounds if tdir_t < tdir * sol.t[1]
    if deriv == 0 && idxs === nothing
      return f.h(val, p, t)
    elseif idxs === nothing
      return f.h(val, p, t, Val{deriv})
    elseif deriv == 0
      return f.h(val, p, t; idxs = idxs)
    else
      return f.h(val, p, t, Val{deriv}; idxs = idxs)
    end
  end

  tdir_solt = tdir * sol.t[end]
  if tdir_t <= tdir_solt
    return sol.interp(val, t, idxs, Val{deriv}, p)
  end

  # history function is evaluated at time point past the final time point of
  # the current solution
  if tdir_t - tdir_solt > 10 * eps(tdir_t)
    f.isout = true
  end

  if integrator.t == sol.prob.tspan[1]
    # handle extrapolations at initial time point
    return constant_extrapolant!(val, t, integrator, idxs, Val{deriv})
  else
    return StochasticDiffEq.current_interpolant!(val, t, integrator, idxs, Val{deriv})
    #TODO only difference form DelayDiffEq here
  end
end


text:
#=
  Changepoint detection of text messages counts in webppl.

  From http://nbviewer.jupyter.org/github/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter1_Introduction/Ch1_Introduction_PyMC3.ipynb
  """
  You are given a series of daily text-message counts from a user of your system.
  The data, plotted over time, appears in the chart below. You are curious to
  know if the user's text-messaging habits have changed over time, either
  gradually or suddenly. How can you model this? (This is in fact my own
  text-message data. Judge my popularity as you wish.)
  """

  Identify the changepoint in text message counts.

  Cf changepoint_coal_miners.jl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function change_point_text_messages(y)
    len = length(y)
    alpha = mean(y) # 1/0.581151832460733  # 1/mean(data)
    lambda_1 ~ Exponential(alpha)
    lambda_2 ~ Exponential(alpha)
    tau ~ DiscreteUniform(1,len)

    # Data is according to Poisson with two different lambdas,
    # one before the change point (tau) and one after the change point
    for i in 1:len
        if tau > i
            y[i] ~ Poisson(lambda_1)
        else
            y[i] ~ Poisson(lambda_2)
        end
    end

end

# http://nbviewer.jupyter.org/github/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter1_Introduction/Ch1_Introduction_PyMC3.ipynb
data =[13, 24, 8, 24, 7, 35, 14, 11, 15, 11, 22, 22, 11, 57, 11, 19, 29, 6,
        19, 12, 22, 12, 18, 72, 32, 9, 7, 13, 19, 23, 27, 20, 6, 17, 13, 10,
        14, 6, 16, 15, 7, 2, 15, 15, 19, 70, 49, 7, 53, 22, 21, 31, 19, 11,
        18, 20, 12, 35, 17, 23, 17, 4, 2, 31, 30, 13, 27, 0, 39, 37, 5, 14,
        13, 22];

# display(plot(data))
model = change_point_text_messages(data)

#=
println("Prior distribution")
chns_prior = sample(model, Prior(), MCMCThreads(), 1000, num_chns)
display(chns_prior)
display(plot(chns_prior))
=#

println("\nPosterior distribution")
num_chns = 4

# 2.06s tau 42+/- 8.37 rhat ~1.5
# chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)
chns = sample(model, MH(), MCMCThreads(), 30_000, num_chns)

# 1.94s tau 37.5+/- 21.2 (rhat ok)
# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)

# chns = sample(model, PG(20), MCMCThreads(), 10_000, num_chns,drop_warmup=true)

# chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)

## LoadError: MethodError: Cannot `convert` an object of type Float64 to an object of type Symbol
# chns = sample(model, Gibbs(NUTS(0.65,5,1,1,0,:lambda_1,:lambda_2), PG(20,:tau) ), 1000)
# chns = sample(model, Gibbs(HMC(0.1,5,:lambda_1,:lambda_2), PG(20,:tau) ), MCMCThreads(), 1000, num_chns)

# Error
# chns = sample(model, Gibbs(HMCDA(0.15,0.65,0.3,:lambda_1,:lambda_2), PG(20,:tau) ), MCMCThreads(), 1000, num_chns) # Error...

display(chns)
# plot(chns)
show_var_dist_pct(chns, :tau)
show_var_dist_pct(chns, :lambda_1,20)
show_var_dist_pct(chns, :lambda_2,20)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Base.Test
using Lora

println("    Testing MHState constructors...")

v = Float32(1.5)
pstate = BasicContUnvParameterState(v, [:accept], [true])
sstate = MHState(deepcopy(pstate))

@test eltype(sstate) == BasicContUnvParameterState{eltype(v)}
@test isequal(sstate.pstate, pstate)
@test sstate.tune.accepted == 0
@test sstate.tune.proposed == 0
@test isnan(sstate.tune.rate)
@test isnan(sstate.ratio)

v = Float64[-6.55, 2.8]
pstate = BasicContMuvParameterState(v)
sstate = MHState(deepcopy(pstate), VanillaMCTune(10, 100, 0.1))

@test eltype(sstate) == BasicContMuvParameterState{eltype(v)}
@test isequal(sstate.pstate, pstate)
@test sstate.tune.accepted == 10
@test sstate.tune.proposed == 100
@test sstate.tune.rate == 0.1
@test isnan(sstate.ratio)

v = Float16[3.16, -2.97, -8.53]
pstate = BasicContMuvParameterState(v)
sstate = MHState(deepcopy(pstate), VanillaMCTune(), 0.27)

@test eltype(sstate) == BasicContMuvParameterState{eltype(v)}
@test isequal(sstate.pstate, pstate)
@test sstate.tune.accepted == 0
@test sstate.tune.proposed == 0
@test isnan(sstate.tune.rate)
@test sstate.ratio == 0.27