text:
# #############################################################################
# #############################################################################
#Cavity refractive index distribution

#This file contains functions for specifying the refractive index distribution of the cavity, by subtypes of the abstract type RefractiveIndex. The following are the available constructors, which can take both positional or keyword arguments:

#   UniformIdx(base::Float64=2.0)
#   constructs a uniform distribution of the specified base value (default 2.0).

#   GaussianIdx(base::Float64=2.0,center::Array{Float64,1}=[0.0,0.0], amplitude::Float64=-0.2,spread::Array{Float64,1}=[0.5,0.5],rotation::Float64=0.0)
#   constructs a distribution with a Gaussian perturbation of specified center position (default at (x,y)=(0.0,0.0)), amplitude (default -0.2), standard deviation of spread in orthogonal directions (default dx=dy=0.5), and angle of rotation (default=0.0), over a specified base value (default 2.0).

#   TaylorIdx(coeff::Array{Float64,2}=[2.0 0.0;0.0 0.4])
#   constructs a general distribution specified by the Taylor expansion coefficients with respect to the x and y coordinates (default is n = 2.0 + 0.4*x*y).

#   ChebFourIdx(coeff_1::Array{Float64,1}=[2.0],coeff_c::Array{Float64,2}=[0.1 -0.1; -0.1 0.1],coeff_s::Array{Float64,2}=Array(Float64,(0,0)))
#   constructs a general differentiable distribution specified by the coefficients for the cross terms of the odd Chebyshev expansion of the radial dependence and the Fourier expansion of the angular dependence. The even terms of the Chebyshev expansion are omitted because they are nonzero at the origin, which would make the distribution discontinuous at the origin if there is also a theta dependence.
#   (Default is n = 2.0*T_0 + 0.1*T_1(r)*cos(theta) - 0.1*T_1(r)*cos(2*theta) - 0.1*T_3(r)*cos(theta) + 0.1*T_3(r)*cos(2*theta))

#   SumIdx(idxs...)
#   constructs a distribution that is the sum of the distributions specified by the arguments. Useful for multiple independent Gaussian perturbations.

#Information about the index distribution in polar coordinates can be extracted by the following functions, which dispatch according to the distribution function subtype to the appropriate calculations:

#   n = nfunc(idx::RefractiveIndex,r,theta)
#   returns the value of the refractive index distribution <idx> at the polar position (<r>,<theta>). <r> can be a Float64 or Array{Float64,N}, but <theta> must follow its type and dimensions, and the result <n> will also follow its type and dimensions.

#   n,dr_n,dtheta_n = nderiv(idx::RefractiveIndex,r,theta)
#   returns a tuple (<n>,<dr_n>,<dtheta_n>) of the value <n> and partial derivatives <dr_n> and <dtheta_n> of the refractive index distribution <idx> at the polar position (<r>,<theta>). <r> can be a Float64 or Array{Float64,N}, but <theta> must follow its type and dimensions, and each entry in the tuple result will also follow its type and dimensions. There are no functions to compute the derivatives alone since the calculations involved almost always requires the refractive index value to also be computed.

#   nderiv!(idx::RefractiveIndex,r::Float64,theta::Float64,result::Array{Float64,1})
#   This is the in-place version of nderiv storing the results as the first three entries of the array <result>.

# #############################################################################
# #############################################################################
#Initialize

#Dependencies
#require("util.jl")
#require("boundary.jl")
import Base: ==, hash

#General index distribution type
abstract RefractiveIndex

# #############################################################################
# #############################################################################
#Hashing and equality

#Extend base hash function with method to hash equivalent RefractiveIndex objects identically.
hash(idx::RefractiveIndex) = fieldhashing(idx)
hash(idx::RefractiveIndex,h::UInt64) = hash(hash(idx),h)

#Equality is equality of type and field values
==(idx1::RefractiveIndex,idx2::RefractiveIndex) = fieldequality(idx1,idx2)

# #############################################################################
# #############################################################################
#Functions for generic index distribution

#These functions are vectorized/in-place versions of nfunc and nderiv for a general RefractiveIndex. The implementations for specific index distribution subtypes, if present, are slightly more efficient.

#Index value ------------------------------------------------------------------
function nfunc{N}(idx::RefractiveIndex,r::Array{Float64,N},theta::Array{Float64,N})
    #assert(size(r)==size(theta))
    #Initialize array
    n::Array{Float64,N} = similar(r)
    #Set value to nfunc output
    for i=1:length(r)
        n[i] = nfunc(idx,r[i],theta[i])
    end
    return n::Array{Float64,N}
end

#Index value and derivatives --------------------------------------------------
function nderiv{N}(idx::RefractiveIndex,r::Array{Float64,N},theta::Array{Float64,N})
    #assert(size(r)==size(theta))
    #Initialize arrays
    n::Array{Float64,N} = similar(r)
    dr_n::Array{Float64,N} = similar(r)
    dtheta_n::Array{Float64,N} = similar(r)
    #Set values to nderiv output
    for i=1:length(theta)
        n[i],dr_n[i],dtheta_n[i] = nderiv(idx,r[i],theta[i])
    end
    return (n,dr_n,dtheta_n)
end

#In-place storing version of nderiv
function nderiv!(idx::RefractiveIndex,r::Float64,theta::Float64,result::Array{Float64,1})
    result[1],result[2],result[3] = nderiv(idx,r,theta)
    return nothing
end

# #############################################################################
# #############################################################################
#List of primitive refractive index distributions

# #############################################################################
#Uniform

#Uniform distribution (base value still needs to be specified to determine losses at boundary)
type UniformIdx <: RefractiveIndex
    #Base value
    base::Float64

    #Construct by fieldname
    UniformIdx(;base::Float64=2.0) = new(base)
    #Construct by position
    UniformIdx(base::Float64=2.0) = new(base)
end

#Index value ------------------------------------------------------------------
nfunc(idx::UniformIdx,r::Float64,theta::Float64) = idx.base::Float64
function nfunc{N}(idx::UniformIdx,r::Array{Float64,N},theta::Array{Float64,N})
    #assert(size(r)==size(theta))
    n::Array{Float64,N} = fill(idx.base,size(r))
    return n::Array{Float64,N}
end

#Index value and derivatives --------------------------------------------------
nderiv(idx::UniformIdx,r::Float64,theta::Float64) = (idx.base,0.0,0.0)
function nderiv{N}(idx::UniformIdx,r::Array{Float64,N},theta::Array{Float64,N})
    #assert(size(r)==size(theta))
    return (fill(idx.base,size(r)), zeros(Float64,size(r)), zeros(Float64,size(r)))
end

# #############################################################################
#Gaussian

#A Gaussian perturbation to a uniform base value.
type GaussianIdx <: RefractiveIndex
    #Base value
    base::Float64
    #Center position of the Gaussian pertubation (Array of x and y coordinate)
    center::Array{Float64,1}
    #Amplitude of the perturbation
    amplitude::Float64
    #The spreading of the Gaussian perturbation (Array of x and y stddev values)
    spread::Array{Float64,1}
    #Angle of rotation
    rotation::Float64

    #Construct by fieldname
    function GaussianIdx(;base::Float64=2.0,center::Array{Float64,1}=[0.0,0.0],amplitude::Float64=-0.2,spread::Array{Float64,1}=[0.5,0.5],rotation::Float64=0.0)
        return new(base,center,amplitude,spread,rotation)
    end
    #Construct by position
    function GaussianIdx(base::Float64=2.0,center::Array{Float64,1}=[0.0,0.0],amplitude::Float64=-0.2,spread::Array{Float64,1}=[0.5,0.5],rotation::Float64=0.0)
        return new(base,center,amplitude,spread,rotation)
    end
end

#Index value ------------------------------------------------------------------
function nfunc(idx::GaussianIdx,r::Float64,theta::Float64)
    x::Float64 = r*cos(theta)
    y::Float64 = r*sin(theta)
    dxr::Float64 = (x-idx.center[1])*cos(idx.rotation) - (y-idx.center[2])*sin(idx.rotation)
    dyr::Float64 = (y-idx.center[2])*cos(idx.rotation) + (x-idx.center[1])*sin(idx.rotation)
    n::Float64 = idx.base + idx.amplitude*exp(-0.5*(sq(dxr/idx.spread[1]) + sq(dyr/idx.spread[2])))
    return n
end

#Index value and derivatives --------------------------------------------------
function nderiv(idx::GaussianIdx,r::Float64,theta::Float64)
    x::Float64 = r*cos(theta)
    y::Float64 = r*sin(theta)
    dxr::Float64 = (x-idx.center[1])*cos(idx.rotation) - (y-idx.center[2])*sin(idx.rotation)
    dyr::Float64 = (y-idx.center[2])*cos(idx.rotation) + (x-idx.center[1])*sin(idx.rotation)
    gaussian::Float64 = idx.amplitude*exp(-0.5*(sq(dxr/idx.spread[1]) + sq(dyr/idx.spread[2])))

    n::Float64 = idx.base + gaussian
    dr_n::Float64 = -(cos(theta+idx.rotation)*dxr/sq(idx.spread[1]) + sin(theta+idx.rotation)*dyr/sq(idx.spread[2]))*gaussian
    dtheta_n::Float64 = ((y*cos(idx.rotation)+x*sin(idx.rotation))*dxr/sq(idx.spread[1]) + (-x*cos(idx.rotation)+y*sin(idx.rotation))*dyr/sq(idx.spread[2]))*gaussian
    return (n,dr_n,dtheta_n)
end

# #############################################################################
#Taylor

#Generic distribution generated by Tayor expansion along the x and y directions about the origin
type TaylorIdx <: RefractiveIndex
    #coeff is a 2D array where the (i,j)th element is the coefficient of the x^(i-1)*y^(j-1) term in the Taylor expansion about the origin of the index distribution. coeff needs only to be large enough to include the significant coefficients
    coeff::Array{Float64,2}

    #Construct by fieldname
    TaylorIdx(;coeff::Array{Float64,2}=[2.0 0.0;0.0 0.4]) = new(coeffmatrix(coeff))
    #Construct by position
    TaylorIdx(coeff::Array{Float64,2}=[2.0 0.0;0.0 0.4]) = new(coeffmatrix(coeff))
end

#Index value ------------------------------------------------------------------
function nfunc(idx::TaylorIdx,r::Float64,theta::Float64)
    x::Float64 = r*cos(theta)
    y::Float64 = r*sin(theta)
    n::Float64 = 0.0
    for i = 1:size(idx.coeff,1)
        for j = 1:size(idx.coeff,2)
            n += idx.coeff[i,j]*gsl_pow_int(y,i-1)*gsl_pow_int(x,j-1)
        end
    end
    return n::Float64
end

function nfunc{N}(idx::TaylorIdx,r::Array{Float64,N},theta::Array{Float64,N})
    assert(size(r)==size(theta))
    n::Array{Float64,N} = similar(r)
    for k = 1:length(r)
        x::Float64 = r[k]*cos(theta[k])
        y::Float64 = r[k]*sin(theta[k])
        n[k] = 0.0
        for i = 1:size(idx.coeff,1)
            for j = 1:size(idx.coeff,2)
                n[k] += idx.coeff[i,j]*gsl_pow_int(y,i-1)*gsl_pow_int(x,j-1)
            end
        end
    end
    return n::Array{Float64,N}
end

#Index value and derivatives --------------------------------------------------
function nderiv(idx::TaylorIdx,r::Float64,theta::Float64)
    x::Float64 = r*cos(theta)
    y::Float64 = r*sin(theta)
    n::Float64 = 0.0
    dr_n::Float64 = 0.0
    dtheta_n::Float64 = 0.0
    for i = 1:size(idx.coeff,1)
        for j = 1:size(idx.coeff,2)
            n += idx.coeff[i,j]*gsl_pow_int(y,i-1)*gsl_pow_int(x,j-1)
            if i >= 2
                dr_n += idx.coeff[i,j]*(i-1)*gsl_pow_int(y,i-2)*gsl_pow_int(x,j-1)*sin(theta)
                dtheta_n += idx.coeff[i,j]*(i-1)*gsl_pow_int(y,i-2)*gsl_pow_int(x,j)
            end
            if j >= 2
                dr_n += idx.coeff[i,j]*gsl_pow_int(y,i-1)*(j-1)*gsl_pow_int(x,j-2)*cos(theta)
                dtheta_n -= idx.coeff[i,j]*gsl_pow_int(y,i)*(j-1)*gsl_pow_int(x,j-2)
            end
        end
    end
    return (n,dr_n,dtheta_n)
end

function nderiv{N}(idx::TaylorIdx,r::Array{Float64,N},theta::Array{Float64,N})
    assert(size(r)==size(theta))
    n::Array{Float64,N} = similar(r)
    dr_n::Array{Float64,N} = similar(r)
    dtheta_n::Array{Float64,N} = similar(r)
    for k = 1:length(r)
        x::Float64 = r[k]*cos(theta[k])
        y::Float64 = r[k]*sin(theta[k])
        n[k] = 0.0
        dr_n[k] = 0.0
        dtheta_n[k] = 0.0
        for i = 1:size(idx.coeff,1)
            for j = 1:size(idx.coeff,2)
                n[k] += idx.coeff[i,j]*gsl_pow_int(y,i-1)*gsl_pow_int(x,j-1)
                if i >= 2
                    dr_n[k] += idx.coeff[i,j]*(i-1)*gsl_pow_int(y,i-2)*gsl_pow_int(x,j-1)*sin(theta[k])
                    dtheta_n[k] += idx.coeff[i,j]*(i-1)*gsl_pow_int(y,i-2)*gsl_pow_int(x,j)
                end
                if j >= 2
                    dr_n[k] += idx.coeff[i,j]*gsl_pow_int(y,i-1)*(j-1)*gsl_pow_int(x,j-2)*cos(theta[k])
                    dtheta_n[k] -= idx.coeff[i,j]*gsl_pow_int(y,i)*(j-1)*gsl_pow_int(x,j-2)
                end
            end
        end
    end
    return (n,dr_n,dtheta_n)
end

# #############################################################################
#Chebyshev-Fourier

#General distribution generated by the product of Chebyshev and Fourier expansion terms
type ChebFourIdx <: RefractiveIndex
    #coeff_1 is a 1D array with the ith entry the coefficient for the T_{i-1}(r) term.
    coeff_1::Array{Float64,1}
    #Necessary condition for continuity at origin: coefficients for even orders of i in T_i(r)*cos(j*theta) or T_i(r)*sin(j*theta) must be 0, since the value of T_i(0) is nonzero, which would produce a discontinuity at the origin. So we define: coeff_c is a 2D array with the (i,j)th entry the coefficient for the T_{2*i-1}(r)*cos(j*theta) term, while coeff_s is a 2D array with the (i,j)th entry the coefficient for the T_{2*i-1}(r)*sin(j*theta) term
    coeff_c::Array{Float64,2}
    coeff_s::Array{Float64,2}

    #Construct by fieldname
    ChebFourIdx(;coeff_1::Array{Float64,1}=[2.0],coeff_c::Array{Float64,2}=[0.1 -0.1; -0.1 0.1],coeff_s::Array{Float64,2}=Array(Float64,(0,0))) = new(coeffmatrix(coeff_1),coeffmatrix(coeff_c),coeffmatrix(coeff_s))
    #Construct by position
    ChebFourIdx(coeff_1::Array{Float64,1}=[2.0],coeff_c::Array{Float64,2}=[0.1 -0.1; -0.1 0.1],coeff_s::Array{Float64,2}=Array(Float64,(0,0))) = new(coeffmatrix(coeff_1),coeffmatrix(coeff_c),coeffmatrix(coeff_s))
end

#Index value ------------------------------------------------------------------
function nfunc(idx::ChebFourIdx,r::Float64,theta::Float64)
    #Computing Chebyshev terms using recurrence relation
    n::Float64 = 0.0

    #Begin Tval at T_0
    Tval_2::Float64 = 2.0*sq(r)-1; Tval_1::Float64 = r; Tval::Float64=1.0

    for i = 1:max(length(idx.coeff_1),2*size(idx.coeff_c,1),2*size(idx.coeff_s,1))
        #Compute Chebyshev value T_{i-1}(r) using recurrence relation
        Tval = 2.0*r*Tval_1 - Tval_2
        Tval_2 = Tval_1; Tval_1 = Tval

        if i <= length(idx.coeff_1)
            n += idx.coeff_1[i]*Tval
        end
        if mod(i,2) == 0
            for j = 1:max(size(idx.coeff_c,2),size(idx.coeff_s,2))
                if i <= 2*size(idx.coeff_c,1) && j <= size(idx.coeff_c,2)
                    n += idx.coeff_c[i/2,j]*Tval*cos(j*theta)
                end
                if i <= 2*size(idx.coeff_s,1) && j <= size(idx.coeff_s,2)
                    n += idx.coeff_s[i/2,j]*Tval*sin(j*theta)
                end
            end
        end
    end
    return n::Float64
end

function nfunc{N}(idx::ChebFourIdx,r::Array{Float64,N},theta::Array{Float64,N})
    assert(size(r)==size(theta))
    #Computing Chebyshev terms using recurrence relation
    n::Array{Float64,N} = zeros(Float64,size(r))

    #Begin Tval at T_0
    Tval_2::Array{Float64,N} = 2.0*sq(r)-1.0; Tval_1::Array{Float64,N} = r; Tval::Array{Float64,N}=ones(Float64,size(r))

    for i = 1:max(length(idx.coeff_1),2*size(idx.coeff_c,1),2*size(idx.coeff_s,1))
        #Compute Chebyshev value T_{i-1}(r) using recurrence relation
        Tval = 2.0*r.*Tval_1 - Tval_2
        Tval_2 = Tval_1; Tval_1 = Tval

        if i <= length(idx.coeff_1)
            n += idx.coeff_1[i]*Tval
        end
        if mod(i,2) == 0
            for j = 1:max(size(idx.coeff_c,2),size(idx.coeff_s,2))
                if i <= 2*size(idx.coeff_c,1) && j <= size(idx.coeff_c,2)
                    n += idx.coeff_c[i/2,j]*Tval.*cos(j*theta)
                end
                if i <= 2*size(idx.coeff_s,1) && j <= size(idx.coeff_s,2)
                    n += idx.coeff_s[i/2,j]*Tval.*sin(j*theta)
                end
            end
        end
    end
    return n::Array{Float64,N}
end

#Index value and derivatives --------------------------------------------------
function nderiv(idx::ChebFourIdx,r::Float64,theta::Float64)
    n::Float64 = 0.0
    dr_n::Float64 = 0.0
    dtheta_n::Float64 = 0.0

    #Begin Tval and Uval at T_0 and U_{-1} respectively
    Tval_2::Float64 = 2.0*sq(r)-1; Tval_1::Float64 = r; Tval::Float64 = 1.0
    Uval_2::Float64 = -2.0*r; Uval_1::Float64 = -1.0; Uval::Float64 = 0.0

    for i = 1:max(length(idx.coeff_1),2*size(idx.coeff_c,1),2*size(idx.coeff_s,1))
        #Compute Chebyshev value T_{i-1}(r) using recurrence relation
        Tval = 2.0*r*Tval_1 - Tval_2
        Tval_2 = Tval_1; Tval_1 = Tval
        #Compute (second kind) Chebyshev value U_{i-1}(r) using recurrence relation, to compute d(T_i)/dx = i*U_{i-1}
        Uval = 2.0*r*Uval_1 - Uval_2
        Uval_2 = Uval_1; Uval_1 = Uval

        if i <= length(idx.coeff_1)
            n += idx.coeff_1[i]*Tval
            dr_n += idx.coeff_1[i]*(i-1)*Uval; #Uval is one term behind
        end
        if mod(i,2) == 0
            for j = 1:max(size(idx.coeff_c,2),size(idx.coeff_s,2))
                if i <= 2*size(idx.coeff_c,1) && j <= size(idx.coeff_c,2)
                    n += idx.coeff_c[i/2,j]*Tval*cos(j*theta)
                    dr_n += idx.coeff_c[i/2,j]*(i-1)*Uval*cos(j*theta)
                    dtheta_n -= idx.coeff_c[i/2,j]*Tval*j*sin(j*theta)
                end
                if i <= 2*size(idx.coeff_s,1) && j <= size(idx.coeff_s,2)
                    n += idx.coeff_s[i/2,j]*Tval*sin(j*theta)
                    dr_n += idx.coeff_s[i/2,j]*(i-1)*Uval*sin(j*theta)
                    dtheta_n += idx.coeff_s[i/2,j]*Tval*j*cos(j*theta)
                end
            end
        end
    end
    return (n,dr_n,dtheta_n)
end

function nderiv{N}(idx::ChebFourIdx,r::Array{Float64,N},theta::Array{Float64,N})
    assert(size(r)==size(theta))
    n::Array{Float64,N} = similar(r)
    dr_n::Array{Float64,N} = similar(r)
    dtheta_n::Array{Float64,N} = similar(r)

    for k=1:length(r)
        n[k] = 0.0
        dr_n[k] = 0.0
        dtheta_n[k] = 0.0
        #Begin Tval and Uval at T_0 and U_{-1} respectively
        Tval_2::Float64 = 2.0*sq(r[k])-1.0; Tval_1::Float64 = r[k]; Tval::Float64 = 1.0
        Uval_2::Float64 = -2.0*r[k]; Uval_1::Float64 = -1.0; Uval::Float64 = 0.0

        for i = 1:max(length(idx.coeff_1),2*size(idx.coeff_c,1),2*size(idx.coeff_s,1))
            #Compute Chebyshev value T_{i-1}(r) using recurrence relation
            Tval = 2.0*r[k]*Tval_1 - Tval_2
            Tval_2 = Tval_1; Tval_1 = Tval
            #Compute (second kind) Chebyshev value U_{i-1}(r) using recurrence relation, to compute d(T_i)/dx = i*U_{i-1}
            Uval = 2.0*r[k]*Uval_1 - Uval_2
            Uval_2 = Uval_1; Uval_1 = Uval

            if i <= length(idx.coeff_1)
                n[k] += idx.coeff_1[i]*Tval
                dr_n[k] += idx.coeff_1[i]*(i-1)*Uval; #Uval is one term behind
            end
            if mod(i,2) == 0
                for j = 1:max(size(idx.coeff_c,2),size(idx.coeff_s,2))
                    if i <= 2*size(idx.coeff_c,1) && j <= size(idx.coeff_c,2)
                        n[k] += idx.coeff_c[i/2,j]*Tval*cos(j*theta[k])
                        dr_n[k] += idx.coeff_c[i/2,j]*(i-1)*Uval*cos(j*theta[k])
                        dtheta_n[k] -= idx.coeff_c[i/2,j]*j*Tval*sin(j*theta[k])
                    end
                    if i <= 2*size(idx.coeff_s,1) && j <= size(idx.coeff_s,2)
                        n[k] += idx.coeff_s[i/2,j]*Tval*sin(j*theta[k])
                        dr_n[k] += idx.coeff_s[i/2,j]*(i-1)*Uval*sin(j*theta[k])
                        dtheta_n[k] += idx.coeff_s[i/2,j]*j*Tval*cos(j*theta[k])
                    end
                end
            end
        end
    end
    return (n,dr_n,dtheta_n)
end

# #############################################################################
# #############################################################################
#Combinator refractive index distribution

#Sum of refractive index distributions
type SumIdx <: RefractiveIndex
    #List of refractive index distributions
    idx_list::Array{RefractiveIndex,1}

    SumIdx(idxs...) = new([idxs...])
end

function nfunc(idx::SumIdx,r::Float64,theta::Float64)
    n::Float64 = 0.0
    for i=1:length(idx.idx_list)
        n += nfunc(idx.idx_list[i],r,theta)
    end
    return n
end

function nderiv(idx::SumIdx,r::Float64,theta::Float64)
    N::Float64 = 0.0
    dr_N::Float64 = 0.0
    dtheta_N::Float64 = 0.0
    for i=1:length(idx.idx_list)
        n::Float64,dr_n::Float64,dtheta_n::Float64 = nderiv(idx.idx_list[i],r,theta)
        N += n
        dr_N += dr_n
        dtheta_N += dtheta_n
    end
    return (N,dr_N,dtheta_N)
end


text:
# Parallel instance of fwi_objective function # Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
export TWRIOptions

"""
    twri_objective(model, source, dobs; options=Options(), optionswri=TWRIOptions())

Evaluate the time domain Wavefield reconstruction inversion objective function. Returns a tuple with function value and \\
gradient(s) w.r.t to m and/or y. `model` is a `Model` structure with the current velocity model and `source` and `dobs` are the wavelets and \\
observed data of type `judiVector`.

Example
=======

    function_value, gradient = fwi_objective(model, source, dobs)

"""
function twri_objective(model::Model, source::judiVector, dObs::judiVector, y::Union{judiVector, Nothing};
                        options=Options(), optionswri=TWRIOptions())
# fwi_objective function for multiple sources. The function distributes the sources and the input data amongst the available workers.

    results = judipmap(j -> twri_objective(model, source[j], dObs[j], subsample(y, j), subsample(options,j), subsample(optionswri,j)), 1:dObs.nsrc)

    # Collect and reduce gradients
    objective = results[1][1]
    gradientm = results[1][2]
    gradienty = results[1][3]

    for j=2:dObs.nsrc
        ~isnothing(gradientm) && (gradientm .+= results[j][2])
        ~isnothing(gradienty) && (gradienty = [gradienty; results[j][3]])
        objective += results[j][1]
    end

    # first value corresponds to function value, the rest to the gradient
    optionswri.params == :m && return objective, gradientm
    optionswri.params == :y && return objective, gradienty
    optionswri.params == :all && return objective, gradientm, gradienty
    return objective
end

# TWRI options
mutable struct TWRIOptions
    grad_corr::Bool
    comp_alpha::Bool
    weight_fun
    eps
    params
    Invq::String
end

"""
    TWRIOptions
        grad_corr::Bool
        comp_alpha::Bool
        weight_fun
        eps
        params::Symbol
        Invq::String

Options structure for TWRI.

`grad_corr`: Whether to add the gradient correction J'(m0, q)*∇_y

`comp_alpha`: WHether to compute optimal alpha (alpha=1 if not)

`weight_fun`: Whether to apply focusing/weighting function to F(m0)'*y and its norm

`eps`: Epsilon (moise level) value (default=0)

`Invq`: How to compute F'Y, either as full field or as a rank 1 approximation `w(t)*q(x)` using the source wavelet for w

`param`: WHich gradient to compute. Choices are `nothing` (objective only), `:m`, `:y` or `:all`

Constructor
==========

All arguments are optional keyword arguments with the following default values:

TWRIOptions(;grad_corr=false, comp_alpha=true, weight_fun=nothing, eps=0, params=:m)
"""

TWRIOptions(;grad_corr=false, comp_alpha=true,
            weight_fun=nothing, eps=0, params=:m, Invq="standard")=
            TWRIOptions(grad_corr, comp_alpha, weight_fun, eps, params, Invq)

function subsample(opt::TWRIOptions, srcnum::Int)
    eloc = length(opt.eps) == 1 ? opt.eps : opt.eps[srcnum]
    return TWRIOptions(opt.grad_corr, opt.comp_alpha, opt.weight_fun, eloc, opt.params, opt.Invq)
end

ext: jl


text:
#=Chaos Randomization
nickhigham juliacon 2018

here's another idea, Randomization, useful & simple

#Chaoscarries on explaining chaos in his own words:
sth nasty happens ()
in order to avoid that by adding extra noise to the Problem
(in other words we explain deterministic chaos i.e. noise )
@23:05
try to avoid tricky cases (i thought we're modelling chaos!)

1. (it) breaks down when A is not Diagonalizable [tested, spot-on]
Numericall Unstable when A has Nearly Defective

Point:(you have to ) Understand how it words (in addition to having)
theoretical Understanding (of what you gonna get )
(if you do add noise to the problem)

Example: matrix functions
given n x m matrix A

A = Xdiag(λi)*X^-1 => f(A) = X*diag(f(λi))*X^-1
if A is Diagonalizable [positive, squared]
I could compute Eigen Decomposition

lambda is

if that was possible, then f()

    -eigen values stay the same multiply out

    put things better again
    given Diagonalization
    you can gen any function (in the matrix)

        Example: the exponential
        or the square root

=#

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using UnicodePlots
import Pkg
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
using SpikeSynchrony
using Statistics
using JLD
using Distributed
using SharedArrays
using Plots
using UnicodePlots
using Evolutionary
using Distributions
using LightGraphs
using Metaheuristics

using SpikingNN
#using StatsBase
using Random
using SparseArrays
using Revise
##
# Override to function to include a state.
##
SNN.@load_units
unicodeplots()

###
# Network 1.
###

const Ne = 200;
const Ni = 50

# https://github.com/RainerEngelken/JuliaCon2017/blob/master/code/generating_large_random_network_topology.ipynb
# needs less memory (n*k instead of n^2)
function gensparsetopo(n,k,seed)
    # seed random
    #srand(seed)
    p = k/(n-1)
    #A = sprand(Bool,n,n,p)
    #weights = rand(Uniform(n,k),n,n)

end
#gensparsetopo (generic function with 1 method)
# https://gist.github.com/flcong/2eba0189d7d3686ea9633a6d14398931
const seed = 10
const k = 0.5
const this_size = 50
#const Ground_sparse = gensparsetopo(size,k,seed)

using JLD
using Logging
#const ground_weights
function get_constant_gw()
    try
        filename = string("ground_weights.jld")
        ground_weights = load(filename,"ground_weights")

        return ground_weights

    catch e
        #this_size = 50
        ground_weights = rand(Uniform(-2,1),this_size,this_size)
        filename = string("ground_weights.jld")
        save(filename, "ground_weights", ground_weights)
        return ground_weights
    end
end

    #=
    try

        filename = string("ground_synapsese.jld")
        synapses = load(filename,"nsynapse")
    catch e
        low = ConstantRate(0.0)
        high = ConstantRate(0.99)
        switch(t; dt = 1) = (t < Int(T/2)) ? low(t; dt = dt) : high(t; dt = dt)
        nsynapse = QueuedSynapse(Synapse.Alpha())
        T = 1000
        syn_current = [switch(t) for t = 1:T]
        excite!(nsynapse, filter(x -> x != 0, syn_current))
        #ai = [ nsynapse for i in 1:size]
        @show(nsynapse)
        filename = string("ground_synapsese.jld")

        save(filename, "nsynapse", nsynapse)
        println("got here")
    end
    return nsynapse, ground_weights
    =#
    #ground_weights

    #finally
    #    ground_weights = rand(Uniform(-2,1),size,size)
    #    println("gets here no!")
    #
    #    filename = string("ground_weights.jld")
    #    save(filename, "ground_weights", ground_weights)
        # What to do unconditionally when try/catch block exits.
    #    return ground_weights

function syn(T)
    low = ConstantRate(0.0)
    high = ConstantRate(0.99)
    switch(t; dt = 1) = (t < Int(T/2)) ? low(t; dt = dt) : high(t; dt = dt)
    the_synapses = QueuedSynapse(Synapse.Alpha())
    syn_current = [switch(t) for t = 1:T]
    excite!(the_synapses, filter(x -> x != 0, syn_current))
    the_synapses
end

function sim_net_darsnack(weight_gain_factor)
    ##
    # Plan network structure stays constant, only synaptic gain varies.
    #
    ##
    ground_weights = get_constant_gw()

    # neuron parameters
    vᵣ = 0
    τᵣ = 1.0
    vth = 1.0
    weights = ground_weights .* weight_gain_factor
    #@show(weights)
    pop_lif = Population(weights; cell = () -> LIF(τᵣ, vᵣ),
                              synapse = Synapse.Alpha,
                              threshold = () -> Threshold.Ideal(vth))

    T = 1000

    the_synapses = syn(T)
    ai = [ the_synapses for i in 1:this_size]
    spikes = simulate!(pop_lif, T; inputs = ai)
    println("\n simulated: \n")
    rasterplot(spikes)|>display#, label = ["Input 1"])#, "Input 2"])

    return spikes,weights,the_synapses
end

function sim_net_darsnack_learn(weight_gain_factor)
    ##
    # Plan network structure stays constant, only synaptic gain varies.
    #
    ##
    ground_weights = get_constant_gw()

    # neuron parameters
    #vᵣ = 0
    #τᵣ = 1.0
    #vth = 1.0
    weights = ground_weights .* weight_gain_factor
    #@show(weights)
    η₀ = 5.0
    τᵣ = 1.0
    vth = 1.0

    pop = Population(weights; cell = () -> SRM0(η₀, τᵣ),
                              synapse = Synapse.Alpha,
                              threshold = () -> Threshold.Ideal(vth),
                              learner = STDP(0.5, 0.5, size(weights, 1)))

    # create step input currents
    ai = InputPopulation([ConstantRate(0.8) for i in 1:this_size])
    #ai = [ the_synapses for i in 1:this_size]

    # create network
    net = Network(Dict([:input => ai, :pop => pop]))
    connect!(net, :input, :pop; weights = weights, synapse = Synapse.Alpha)

    # simulate
    #w = Float64[]
    T = 1000

    @time output = simulate!(net, T; dense = true)

    #pop_lif = Population(weights; cell = () -> LIF(τᵣ, vᵣ),
    #                          synapse = Synapse.Alpha,
    #                          threshold = () -> Threshold.Ideal(vth))

    the_synapses = syn(T)
    spikes = output[:pop]
    #spikes = simulate!(pop_lif, T; inputs = ai)
    println("\n simulated: \n")
    rasterplot(spikes)|>display#, label = ["Input 1"])#, "Input 2"])

    return spikes,weights,ai
end

function sim_net_darsnack_used(weight_gain_factor)
    ##
    # Plan network structure stays constant, only synaptic gain varies.
    #
    ##
    ground_weights = get_constant_gw()

    # neuron parameters
    vᵣ = 0
    τᵣ = 1.0
    vth = 1.0
    weights = ground_weights .* weight_gain_factor
    #@show(weights)
    #η₀ = 5.0
    #τᵣ = 1.0
    #vth = 1.0

    pop = Population(weights; cell = () ->LIF(τᵣ, vᵣ),
                              synapse = Synapse.Alpha,
                              threshold = () -> Threshold.Ideal(vth))
                              # threshold = Threshold.Ideal
                              #learner = STDP(0.5, 0.5, size(weights, 1)))

    # create step input currents
    ai = InputPopulation([ConstantRate(0.8) for i in 1:this_size])
    #ai = [ the_synapses for i in 1:this_size]

    # create network
    net = Network(Dict([:input => ai, :pop => pop]))
    connect!(net, :input, :pop; weights = weights, synapse = Synapse.Alpha)

    # simulate
    #w = Float64[]
    T = 1000

    @time output = simulate!(net, T; dense = true)

    #pop_lif = Population(weights; cell = () -> LIF(τᵣ, vᵣ),
    #                          synapse = Synapse.Alpha,
    #                          threshold = () -> Threshold.Ideal(vth))

    the_synapses = syn(T)
    spikes = output[:pop]
    #spikes = simulate!(pop_lif, T; inputs = ai)
    println("\n simulated: \n")
    rasterplot(spikes)|>display#, label = ["Input 1"])#, "Input 2"])

    #clear(pop)
    #clear(net)
    #clear(ai)
    #clear(weights)

    return spikes,weights,ai
end

function get_trains_dars(train_dic::Dict)
    valued = [v for v in values(train_dic)]
    keyed = [k for k in keys(train_dic)]
    cellsa = Array{Union{Missing,Any}}(undef, length(keyed))#, Int(last(findmax(valued)[1])))
    #nac = Int(last(findmax(valued)[1]))

    for (inx, cell_id) in enumerate(1:length(keyed))
        cellsa[inx] = []
    end
    @inbounds for cell_id in keys(train_dic)
        @inbounds for time in train_dic[cell_id]
            append!(cellsa[Int(cell_id)], time*ms)
        end
    end
    #@show(cellsa)
    #cellsa
    cellsb = cellsa[:,1]
    cellsb
end

#=
function make_net_SNeuralN()
    weights = rand(Uniform(-2,1),25,25)
    pop = Population(weights; cell = () -> LIF(τᵣ, vᵣ),
                              synapse = Synapse.Alpha,
                              threshold = () -> Threshold.Ideal(vth))
    # create input currents
    low = ConstantRate(0.0)
    high = ConstantRate(0.1499)
    switch(t; dt = 1) = (t < Int(T/2)) ? low(t; dt = dt) : high(t; dt = dt)
    n1synapse = QueuedSynapse(Synapse.Alpha())
    n2synapse = QueuedSynapse(Synapse.Alpha())
    excite!(n1synapse, filter(x -> x != 0, [low(t) for t = 1:T]))
    excite!(n2synapse, filter(x -> x != 0, [switch(t) for t = 1:T]))
    voltage_array_size = size(weights)[1]
    voltages = Dict([(i, Float64[]) for i in 1:voltage_array_size])
    cb = () -> begin
        for id in 1:size(pop)
            push!(voltages[id], getvoltage(pop[id]))
        end
    end
    input1 = [ (t; dt) -> 0 for i in 1:voltage_array_size/3]
    input2 = [ n2synapse for i in voltage_array_size/3+1:2*voltage_array_size/3]
    input3 = [ (t; dt) -> 0 for i in 2*voltage_array_size/3:voltage_array_size]
    input = vcat(input2, input1, input3)
    return input,cb,voltages
end
=#
    #outputs = simulate!(pop, T; cb = cb, inputs=input)

const Ne = 200;
const Ni = 50
const σee = 1.0
const pee = 0.5
const σei = 1.0
const pei = 0.5

global E
global spkd_ground

function make_net_from_graph_structure(xx)#;

    xx = Int(round(xx))
    @show(xx)
    #h = turan_graph(xx, xx)#, seed=1,cutoff=0.3)

    h = circular_ladder_graph(xx)#, xx)#, seed=1,cutoff=0.3)
    hi = circular_ladder_graph(xx)#, seed=1,cutoff=0.3)
    E = SNN.IZ(; N = Ne, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
    I = SNN.IZ(; N = Ni, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
    #EE = SNN.SpikingSynapse(E, E, :v; σ = σee, p = 1.0)
    EI = SNN.SpikingSynapse(E, I, :v; σ = σei, p = 1.0)
    IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1.0)
    II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 1.0)
    # PINningSynapse
    P = [E, I]#, EEA]
    C = [EI, IE, II]#, EEA]
    #EE = SNN.PINningSynapse(E, E, :v; σ=0.5, p=0.8)
    #for n in 1:(N - 1)
    #    SNN.connect!(EE, n, n + 1, 50)
    #end
    #for (i,j) in enumerate(h.fadjlist) println(i,j) end
    EE = SNN.SpikingSynapse(E, E, :v; σ=0.5, p=0.8)

    @inbounds for (i,j) in enumerate(h.fadjlist)
        @inbounds for k in j
            SNN.connect!(EE,i, k, 10)
        end
    end

    @inbounds for (i,j) in enumerate(hi.fadjlist)
        @inbounds for k in j
            if i<Ni && k<Ni

                SNN.connect!(EI,i, k, 10)
                SNN.connect!(IE,i, k, 10)
                SNN.connect!(II,i, k, 10)
            end
        end
    end

    #for (i,j) in enumerate(h.fadjlist)
    #    for k in j
    #        SNN.connect!(EI, i, k, 50)
    #    end
    #end
    P = [E, I]#, EEA]
    C = [EE, EI, IE, II]#, EEA]
    return P, C

end

function make_net_SNN(Ne, Ni; σee = 1.0, pee = 0.5, σei = 1.0, pei = 0.5)
     Ne = 200;
     Ni = 50

     E = SNN.IZ(; N = Ne, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
     I = SNN.IZ(; N = Ni, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
     EE = SNN.SpikingSynapse(E, E, :v; σ = σee, p = pee)
     EI = SNN.SpikingSynapse(E, I, :v; σ = σei, p = pei)
     IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 0.5)
     II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.5)
     P = [E, I]#, EEA]
     C = [EE, EI, IE, II]#, EEA]
     return P, C
end
function get_trains(p)
    fire = p.records[:fire]
    x, y = Float32[], Float32[]
    for time in eachindex(fire)
        for neuron_id in findall(fire[time])
            push!(x, time)
            push!(y, neuron_id)
        end
    end
    cellsa = Array{Union{Missing,Any}}(undef, 1, Int(findmax(y)[1]))
    nac = Int(findmax(y)[1])
    for (inx, cell_id) in enumerate(1:nac)
        cellsa[inx] = []
    end
    @inbounds for cell_id in unique(y)
        @inbounds for (time, cell) in collect(zip(x, y))
            if Int(cell_id) == cell
                append!(cellsa[Int(cell_id)], time)

            end

        end
    end

    cellsa

end

#P, C = make_net(Ne, Ni, σee = 0.5, pee = 0.8, σei = 0.5, pei = 0.8, a = 0.02)
#sggcu =[ CuArray(convert(Array{Float32,1},sg)) for sg in spkd_ground ]

#Flux.SGD
#Flux.gpu

function rmse(spkd)
    error = Losses(mean(spkd),spkd;agg=mean)
end

function rmse_depr(spkd)
    total = 0.0
    @inbounds for i = 1:size(spkd, 1)
        total += (spkd[i] - mean(spkd[i]))^2.0
    end
    return sqrt(total / size(spkd, 1))
end

function raster_difference(spkd0, spkd_found)
    maxi0 = size(spkd0)[2]
    maxi1 = size(spkd_found)[2]
    mini = findmin([maxi0, maxi1])[1]
    spkd = ones(mini)
    maxi = findmax([maxi0, maxi1])[1]

    if maxi > 0
        if maxi0 != maxi1
            return sum(ones(maxi))

        end
        if isempty(spkd_found[1, :])
            return sum(ones(maxi))
        end
    end
    spkd = ones(mini)
    @inbounds for i in eachindex(spkd)
        if !isempty(spkd0[i]) && !isempty(spkd_found[i])
            maxt1 = findmax(spkd0[i])[1]
            maxt2 = findmax(spkd_found[i])[1]
            maxt = findmax([maxt1, maxt2])[1]
            if maxt1 > 0.0 && maxt2 > 0.0
                t, S = SpikeSynchrony.SPIKE_distance_profile(
                    unique(sort(spkd0[i])),
                    unique(sort(spkd_found[i]));
                    t0 = 0.0,
                    tf = maxt,
                )
                spkd[i] = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1])
            end
        end
    end
    spkd
end

#=
using Base.Threads

function raster_difference_threads(spkd0, spkd_found)
    maxi0 = size(spkd0)[2]
    maxi1 = size(spkd_found)[2]
    mini = findmin([maxi0, maxi1])[1]
    spkd = ones(mini)
    maxi = findmax([maxi0, maxi1])[1]
    if maxi > 0
        if maxi0 != maxi1
            return sum(ones(maxi))

        end
        if isempty(spkd_found[1, :])
            return sum(ones(maxi))
        end
    end
    spkd = ones(mini)
    N = length(spkd)
    Threads.foreach(
        Channel(nthreads()) do chnl
            for i in 1:N
                val = 1.0
                if !isempty(spkd0[i]) && !isempty(spkd_found[i])
                    maxt1 = findmax(spkd0[i])[1]
                    maxt2 = findmax(spkd_found[i])[1]
                    maxt = findmax([maxt1, maxt2])[1]
                    if maxt1 > 0.0 && maxt2 > 0.0
                        t, S = SpikeSynchrony.SPIKE_distance_profile(
                            unique(sort(spkd0[i])),
                            unique(sort(spkd_found[i]));
                            t0 = 0.0,
                            tf = maxt,
                        )
                        val = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1])
                    end
                end
                put!(chnl, val)
                println("thread ", threadid(), ": ", "put ", val)
            end
        end
    end
    spkd
end
=#

#=
function loss(model)
    @show(Ne,Ni)
    @show(model)
    P1, C1 = make_net_SNN(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])
        SNN.sim!(P1, C1, 1ms)
    end
    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")
    SNN.raster([E1]) |> display
    error = raster_difference(spkd_ground, spkd_found)
    error = sum(error)
    @show(error)

    error

end
=#

function loss(model)
    @show(Ne,Ni)
    @show(model)

    σee = model[1]
    pee = model[2]
    σei = model[3]
    pei = model[4]
    P1, C1 = make_net_SNN(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")
    SNN.raster([E1]) |> display

    error = raster_difference(spkd_ground, spkd_found)
    error = sum(error)
    @show(error)

    error

end

function eval_best(params)
    #xx = Int(round(params[1]))
    #@show(xx)
    #P1, C1 = make_net_SNN(xx)#,a=a)
    #println("found ",xx)

    σee = params[1]
    pee = params[2]
    σei = params[3]
    pei = params[4]
    P1, C1 = make_net_SNN(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])#vec(E_stim[t,:])#[i]#3randn(Ne)
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth: \n")
    SNN.raster([E]) |> display
    println("candidate: \n")

    SNN.raster([E1]) |> display
    #error = raster_difference(spkd_ground,spkd_found)
    E1, spkd_found

end

#using Evolutionary, MultivariateStats
#range = Any[lower,upper]

# mutation = domainrange(fill(0.5,4))
#
#function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::GA, options)
#    idx = sortperm(state.fitpop)
#record["population"] = population
#    state.fitpop[idx[1:5]]
#    record["fitpop"] = state.fitpop[idx[1:5]]
#end

#using Evolutionary

function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::GA, options)
    idx = sortperm(state.fitpop)
    record["fitpop"] = state.fitpop[:]#idx[1:last(idx)]]
    record["pop"] = population[:]
    #record["σ"] = state.
end
#=
function Evolutionary.value!(::Val{:serial}, fitness, objfun, population::AbstractVector{IT}) where {IT}
    @show(typeof(fitness))
    fitness = SharedArrays.SharedArray{Float32}(fitness)
    @time @sync @distributed for i in 1:length(population)

    #pmap(f, [::AbstractWorkerPool], c...; distributed=true, batch_size=) -> collection
        fitness[i] = value(objfun, population[i])
        #println("I'm worker $(myid()), working on i=$i")
    end
    @show(typeof(fitness))

    fitness
    #fitness = Array(fitness)
    #@show(fitness)
end
=#


text:
using Random

import Base.rand

include("xos256p.jl")

"""
return a random Float64 in [0, 1)
"""
function rand(xos::Xoshiro256p)
    r = next(xos)
    return r/(UInt64(0)-1) # return a random UInt64 divided by the largest UInt64 possible (UInt64(0) - 1)
end

"""
Seeds a given random number generator with seed x.
"""
function srand(rng::Xoshiro256p,seed::Vector{UInt64})
    for i = 1:6
        rng.Cg[i] = rng.Bg[i] = rng.Ig[i] = seed[i]
    end
end

"""
Resets a given random number generator to the beginning of the current stream.
"""
function reset_stream!(rng::Xoshiro256p)
    for i = 1:4
        rng.Cg[i] = rng.Bg[i] = rng.Ig[i]
    end
end

function reset_substream!(rng::Xoshiro256p)
    for i = 1:4
        rng.Cg[i] = rng.Bg[i]
    end
end

function next_substream!(rng::Xoshiro256p)
    short_jump(rng)
end

mutable struct Xoshiro256plusGen <: AbstractRNGStream
    nextSeed::Vector{UInt64}
    function Xoshiro256plusGen(x::Vector{UInt64})
        new(copy(x))
    end
end

function show(io::IO,rng_gen::Xoshiro256plusGen)
    print(io, "Seed for next Xoshiro256plus generator:\n$(rng_gen.nextSeed)")
end

function show(io::IO,rng::Xoshiro256p)
    print(io,"Full state of Xoshiro256plus generator:\nCg = $(rng.Cg)\nBg = $(rng.Bg)\nIg = $(rng.Ig)")
end

get_state(rng_gen::MRG32k3aGen) = copy(rng_gen.nextSeed)
get_state(rng::MRG32k3a) = copy(rng.Cg)

"""
Reseeds the RNG generator object.
"""
function srand(rng_gen::Xoshiro256plusGen,seed::Vector{UInt64})
    for i = 1:4
        rng_gen.nextSeed[i] = seed[i]
    end
end

"""
Given an RNG generator object, returns the next RNG stream.
"""
function next_stream(rng_gen::Xoshiro256plusGen)
    rng = Xoshiro256p(copy(rng_gen.nextSeed))
    rng2 = Xoshiro256p(copy(rng_gen.nextSeed))
    long_jump(rng)
    rng_gen.nextSeed[:] = rng.Cg
    return rng2
end


text:
negate(term::Const) = Const(Symbol(:¬, term.name))
negate(term::Compound) = Compound(Symbol(:¬, term.name), term.args)

function reify_negations(term::Compound, domain::Domain)
    if term.name == :not && is_pred(term.args[1], domain)
        negate(term.args[1])
    else
        Compound(term.name, [reify_negations(a, domain) for a in term.args])
    end
end
reify_negations(term::Var, domain) =
    term
reify_negations(term::Const, domain) =
    term

function get_negation_clauses(domain::Domain)
    return map(values(get_predicates(domain))) do sig
        pred = convert(Term, sig)
        Clause(negate(pred), [Compound(:not, [pred])])
    end
end


text:
"""
    SAA

Default `SampledSolver`. Generates a `StochasticSolution` using the sample average approximation (SAA) method, to the desired confidence level.
"""
struct SAA{S <: SPSolverType} <: AbstractSampledSolver
    internal_solver::S

    function SAA(solver::SPSolverType)
        if isa(solver, JuMP.UnsetSolver)
            error("Cannot solve emerging SAA problems without functional solver.")
        end
        S = typeof(solver)
        return new{S}(solver)
    end
end
"""
    SAA(; solver::SPSolverType = JuMP.UnsetSolver())

Return an SAA where the emerging SAA problems are solved using `solver`.
"""
function SAA(; solver::SPSolverType = JuMP.UnsetSolver())
    return SAA(solver)
end

mutable struct SAAModel{M <: StochasticModel, S <: SPSolverType} <: AbstractSampledModel
    stochasticmodel::M
    solver::S
    solution::StochasticSolution
    saa::StochasticProgram

    function SAAModel(stochasticmodel::StochasticModel, solver::SPSolverType)
        M = typeof(stochasticmodel)
        S = typeof(solver)
        return new{M, S}(stochasticmodel, solver, EmptySolution())
    end
end

function SampledModel(stochasticmodel::StochasticModel, solver::SAA)
    return SAAModel(stochasticmodel, solver.internal_solver)
end

function optimize_sampled!(saamodel::SAAModel, sampler::AbstractSampler, confidence::AbstractFloat; M::Integer = 10, T::Integer = 10, Ñ::Integer = 1000, tol::AbstractFloat = 1e-2, Ninit::Int = 16, Nmax::Integer = 5000, solver_config::Function = (solver,N)->nothing, log = true)
    sm = saamodel.stochasticmodel
    solver = saamodel.solver
    N = Ninit
    α = 1-confidence
    progress = ProgressThresh(tol, 0.0, "SAA gap")
    log && ProgressMeter.update!(progress, Inf,
                                 showvalues = [
                                     ("Confidence interval", NaN),
                                     ("Relative error", Inf),
                                     ("Sample size", NaN),
                                     ("Current sample size", N)
                                 ])
    while true
        CI = confidence_interval(sm, sampler; solver = solver, confidence = 1-α, N = N, M = M, Ñ = max(N, Ñ), T = T, log = log, keep = false, offset = 6, indent = 4)
        Q = (upper(CI) + lower(CI))/2
        gap = length(CI)/abs(Q+1e-10)
        log && ProgressMeter.update!(progress, gap,
                                     showvalues = [
                                         ("Confidence interval", CI),
                                         ("Relative error", gap),
                                         ("Sample size", N),
                                         ("Current sample size", 2*N)
                                     ])
        if gap <= tol
            sp = sample(sm, sampler, N)
            optimize!(sp, solver = solver)
            Q = optimal_value(sp)
            while !(Q ∈ CI)
                sp = sample(sm, sampler, N)
                optimize!(sp, solver = solver)
                Q = optimal_value(sp)
            end
            saamodel.solution = StochasticSolution(optimal_decision(sp), Q, N, CI)
            saamodel.saa = sp
            return :Optimal
        end
        N = N * 2
        if N > Nmax
            return :LimitReached
        end
        solver_config(solver, N)
    end
end

function internal_solver(solver::SAA)
    return internal_solver(solver.internal_solver)
end

function stochastic_solution(saamodel::SAAModel)
    return saamodel.solution
end


text:
module CUDATests

using Test
using HDF5
using FFTW

using PicoQuant
include("test_utils.jl")

using CUDA
try
    device = CUDA.device()

    @testset "CUDA tests" begin
        @testset "Test state preparation algorithm" begin
            @test begin
                qubits = 3
                circ = create_simple_preparation_circuit(qubits, 4, 42)

                ψ′ = CuArray(get_statevector_using_qiskit(circ, big_endian=true))

                tng = convert_qiskit_circ_to_network(circ, InteractiveBackend{CuArray{ComplexF64}}())
                add_input!(tng, "000")
                out_node = full_wavefunction_contraction!(tng, "vector")

                ψ = load_tensor_data(tng, out_node)

                overlap = (a, b) -> abs(transpose(conj.(a)) * b)

                overlap(ψ, ψ′) ≈ 1.
            end

            @test begin
                qubits = 3
                circ = create_simple_preparation_circuit(qubits, 4, 42)

                ψ′ = CuArray(get_statevector_using_qiskit(circ, big_endian=true))

                tng = convert_qiskit_circ_to_network(circ, InteractiveBackend{CuArray{ComplexF64}}())
                add_input!(tng, "000")
                mps_nodes = contract_mps_tensor_network_circuit!(tng)
                calculate_mps_amplitudes!(tng, mps_nodes)

                ψ = load_tensor_data(tng, :result)

                overlap = (a, b) -> abs(transpose(conj.(a)) * b)

                overlap(ψ, ψ′) ≈ 1.
            end
        end
    end
catch
    println("Skipping CUDA tests as CUDA.jl failed to initialize")
end

end


text:
module SfFluidDynamics

using Unitful

import PhysicalConstants.CODATA2018: g_n

import ..SfPhysics: Angle, to_angle
import ..SfPhysics: projectile_displacement, projectile_velocity, projectile_flight_time,
	projectile_peak_displacement, projectile_range, projectile_angle

export drag_force, reynolds_number

const Viscosity{T} = Unitful.AbstractQuantity{T, Unitful.𝐌 * Unitful.𝐋^-1 *  Unitful.𝐓^-1, typeof(u"Pa*s")}

"""
    drag_force(drag_coefficient::Real, ρ_fluid::Unitful.Density, reference_area::Unitful.Area, v::Unitful.Velocity)

Drag acting on a body with `drag_coefficient` and `reference_area` travelling at `v` through a fluid with density `ρ_fluid`.
"""
drag_force(drag_coefficient::Real, ρ_fluid::Unitful.Density, reference_area::Unitful.Area, v::Unitful.Velocity) =
	upreferred(0.5ρ_fluid * v^2 * drag_coefficient * reference_area)

"""
    reynolds_number(v::Unitful.Velocity, ρ_fluid::Unitful.Density, l::Unitful.Length, μ::Viscosity)

Reynolds number of a body with characteristic dimension `l` travelling at velocity `v` through a fluid of density `ρ_fluid` and viscosity `μ`.
"""
reynolds_number(v::Unitful.Velocity, ρ_fluid::Unitful.Density, l::Unitful.Length, μ::Viscosity) = (v * ρ_fluid * l) / μ |>u"s/s"

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
using Microgrids
using CSV, DataFrames

# Importing load and solar data
data = DataFrame(CSV.File("$(@__DIR__)/data/Ouessant_data_2016.csv"))

# Simulation steps
ntimestep = length(data.Load)

# Components parameters
# Project
lifetime = 25
discount_rate = 0.05
timestep = 1
# Load
Pload = data."Load"[1:ntimestep]
# Photovoltaic
power_rated_PV = 4106.82251423571
fPV = 1.
IT = data."Ppv1k"[1:ntimestep] ./ 1000
IS = ones(ntimestep)
investiment_cost_PV = 1200.
om_cost_PV = 20.
replacement_cost_PV = 1200.
salvage_cost_PV = 1200.
lifetime_PV = 25
# Battery
energy_initial = 0.
energy_max = 6839.87944197573
energy_min = 0
power_min = -1.114*energy_max
power_max = 1.002*energy_max
loss = 0.05
investiment_cost_BT = 350.
om_cost_BT = 10.
replacement_cost_BT = 350.
salvage_cost_BT = 350.
lifetime_BT = 15
lifetime_thrpt = 3000
# Diesel generator
power_rated_DG = 1800.
min_load_ratio = 0
F0 = 0.0
F1 = 0.240
fuel_cost = 1.
investiment_cost_DG = 400.
om_cost_DG = 0.02
replacement_cost_DG = 400.
salvage_cost_DG = 400.
lifetime_DG = 15000

# Create microgrid components
project = Project(lifetime, discount_rate, timestep)
dieselgenerator = DieselGenerator(power_rated_DG, min_load_ratio, F0, F1, fuel_cost, investiment_cost_DG, om_cost_DG, replacement_cost_DG, salvage_cost_DG, lifetime_DG)
photovoltaic = Photovoltaic(power_rated_PV, fPV, IT, IS, investiment_cost_PV, om_cost_PV, replacement_cost_PV, salvage_cost_PV, lifetime_PV)
battery = Battery(energy_initial, energy_max, energy_min, power_min, power_max, loss, investiment_cost_BT, om_cost_BT, replacement_cost_BT, salvage_cost_BT, lifetime_BT, lifetime_thrpt)

# Create microgrid
microgrid = Microgrid(project, Pload, dieselgenerator, battery, [photovoltaic])

# Run simulation
results = simulate(microgrid)


text:
# Insert graph data into an adjacency table.
using JLD

# Iterate through files
(@isdefined Nfile) || (Nfile = 8)
for i = 1:Nfile
    inserttime = @elapsed begin

    # Create filename
    fname = joinpath(Base.source_dir(),"data",string(i))
    println("On file: "*string(i))

    # Load associative array.
    #A = loadassoc(fname*".A.jld")
    A = load(fname*".A.jld")["A"]

    # Insert associative array.
    put(Tadj,A)

    #Compute in and out degree.
    Aout_i = putCol(sum(A,2),"OutDeg,")
    Ain_i = putCol(sum(A,1)',"InDeg,")

    # Accumulate in and out degrees.
    put(TadjDeg,Aout_i)
    put(TadjDeg,Ain_i)

    end
    println("Time: "*string(inserttime)*", Edges/sec: "*string((2*nnz(A)+nnz(Aout_i)+nnz(Ain_i))./inserttime))
    #println(['Time: ' num2str(insertTime) ', Edges/sec: ' num2str((2*nnz(A)+nnz(Aout_i)+nnz(Ain_i))./insertTime)]);
end

println("Table entries: "*string(nnz(Tadj)));


text:
using MathProgBase

"""
Computes the estimates of the logit regression of y on x.
Args:
    y (array): (n by 1) array.
    x (array): (n by k) array.
"""
function logit(y::AbstractArray,x::AbstractArray)
    println("Start logit")
    cnv = 0
    b0 = zeros(size(x)[2],1)
    b1 = zeros(size(x)[2],1)
    while cnv==0
        ind = x*b0
        P = (exp.(ind))./(1+exp.(ind))
        grd = sum(repmat(y-P,1,size(x)[2]).*x,1)
        hes = -x'*(repmat(P.*(1-P),1,size(x)[2]).*x)
        b1 = b0 - inv(hes)*grd'
        dev = maximum(abs.(b1.-b0))
        if dev < 1e-8
            println("true")
            cnv = 1
        end
        b0 = b1
    end
    return b1
end

"""
Refines the bounds of the unknown coefficients for the parameter space in warm start approach.
Args:
    x (array): (n by k) matrix of covariate data.
    beta0 (int): Coefficient for the first covariate in x.
    bnd (array): ((k-1) by 2) matrix where the first and second columns respectively
                 store the lower and upper bounds of the unknown coefficients
Yields:
    bound (array): ((k-1) by 2) matrix for the refined lower and upper bounds
                of the unknown coefficients for the parameter space used
                in the warm-start approach.
"""
function getBnd(y::AbstractArray,x::AbstractArray,beta0::Integer,bnd::AbstractArray)
    println("Start get bound")
    k=size(x)[2]-1

    p_hat = 1./(1+exp.(-x*logit(y,x)))

    constr=repmat(p_hat-0.5,1,size(x)[2]).*x
    bound= copy(bnd)

    A = constr[:,2:size(x)[2]]
    rhs = -constr[:,1]*beta0

    #Set Params
    tol=0.000001
    solver = GurobiSolver(OutputFlag=0,OptimalityTol=tol,FeasibilityTol=tol,IntFeasTol=tol)

    lb = bound[:,1]
    ub = bound[:,2]

    for i=1:k
        objcoef=zeros(1,k)
        objcoef[i]=1
        obj = objcoef

        try
            result= linprog(vec(obj), A, '>', vec(rhs), vec(lb), vec(ub), solver)
            bound[i,1]=result.objval
            catch
                println("Error reported")
        end
        lb = bound[:,1]
        try
            result= linprog(-vec(obj), A, '>', vec(rhs), vec(lb), vec(ub), solver)
            bound[i,2]=-result.objval
            catch
                println("Error reported")
        end
        ub = bound[:,2]
    end
    return bound
end

"""
Implements the warm start approach.
Args:
    y (array):      vector of binary outcomes.
    x_foc (array):  (n by k) matrix of focus covariates where the first column should
                    contain data of the continous regressor with respect to the scale normalization.
    x_aux (array):  (n by d) matrix of auxiliary covariates which will be selected based on best
                    subset covariate selection procedure.
    beta0 (int):    Coefficient taking value either 1 or -1 to normalize the scale for the first
                    first covariate in x_foc.
    q (int):        Cardinality constraint for the covariate selection.
    T (int):        Time limit specified for the MIO solver.
    bnd (array):    (((k-1)+d) by 2) matrix where the first and second columns respectively
                    store the lower and upper bounds of the unknown coefficients. First (k-1)
                    rows correspond to the bounds of the focused covariates excluding the first one.
                    Remaining d rows correspond to the bounds of the auxiliary covariates.
    mio (int) :     MIO formulation of the best subset maximum score problem
                    set mio = 1 for Method 1 of the paper
                    set mio = 2 for Method 2 of the paper
    tau (int):      Tuning parameter for enlarging the estimated bounds.
Yields:
    bhat (array):   Maximum score estimates for the unknown coefficients.
    score (array):  Value of maximum score objective function.
    gap   (array):  MIO optimization gap value in case of early termination
                    gap = 0 ==> optimal solution is found within the time limit.
    rtime (array):  Time used by the MIO solver in the estimation procedure.
"""
function WarmStartMaxScore(y::AbstractArray,x_foc::AbstractArray,x_aux::AbstractArray,beta0::Integer,q::Integer,T::Real,tol::Real,bnd::AbstractArray,mio::Integer,tau::Real)
    bnd_h = getBnd(y,[x_foc x_aux],beta0,bnd)
    bnd_abs = tau*maximum(abs.(bnd_h),2)
    bnd0 = [maximum([-bnd_abs bnd[:,1]],2) minimum([bnd_abs bnd[:,2]],2)] # this is the refind bound used for warm start MIO
    bhat,score,gap,rtime,ncount = MaxScore(y,x_foc,x_aux,beta0,q,T,tol,bnd0,mio)
    return bhat,score,gap,rtime,ncount
end

export WarmStartMaxScore


text:
using Distributions

function hyperplanes_poisson(R, λ)::Array{Line}
    """
    Inputs:
        max radius R
        intensity λ

    Outputs:
        N lines with parameters (θ, r) uniformly random,
        in the set [0, 2π)×[0,R), where N=Poi(2πRλ) .
    """
    n = rand(Poisson(2*pi*R*λ), 1)[1]
    α = (2*pi) .* rand(n) ; radii = R * rand(n)
    return Array([Line(α[k], radii[k]) for k in 1:n])
end

function hyperplanes_triangle(R, n::Int)::Array{Line}
    """
    Inputs:
        max radius R
        number of lines n

    Outputs:

    """
    radii = R * rand(n)
    angles = [pi/2, pi/2 + 2*pi/3, pi/2 + 4*pi/3]
    α = rand(angles ∪ -angles, n)
    return Array([Line(α[k], radii[k]) for k in 1:n])
end

function hyperplanes_anisotropic_poisson(R, K, Λ)::Array{Line}
    """
    Inputs:
        - radius R
        - 1d array of directions K
        - 1d array of intensities Lambda with same dimension as k
    """
    @assert length(K)==length(Λ) throw(DimensionMismatch)

    n = rand(Poisson(2*pi*R*Λ[1]), 1)[1]
    output = [Line(K[1], x) for x in 2*R*rand(n).-R]
    if length(K)==1
        return output
    else
        for i in 2:length(K)
            n = rand(Poisson(2*pi*R*Λ[i]), 1)[1]
            newdir = [Line(K[i], x) for x in 2*R*rand(n).-R]
            output = output ∪ newdir
        end
        return output
    end
end

max_stars_repo_name: SimonCoste/LineIntersections

max_issues_repo_name: SimonCoste/LineIntersections

max_forks_repo_name: SimonCoste/LineIntersections


text:
# ------------------------------------------------------------------------------------------
# # Learning to recognize handwritten digits using a neural network
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We have now reached the point where we can tackle a very interesting task: applying the
# knowledge we have gained with machine learning in general, and `Flux.jl` in particular, to
# create a neural network that can recognize handwritten digits! The data are from a data
# set called MNIST, which has become a classic in the machine learning world.
#
# [We could also try to apply the techniques to the original images of fruit instead.
# However, the fruit images are much larger than the MNIST images, which makes the learning
# a suitable neural network too slow.]
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Data munging
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# As we know, the first difficulty with any new data set is locating it, understanding what
# format it is stored in, reading it in and decoding it into a useful data structure in
# Julia.
#
# The original MNIST data is available [here](http://yann.lecun.com/exdb/mnist); see also
# the [Wikipedia page](https://en.wikipedia.org/wiki/MNIST_database). However, the format
# that the data is stored in is rather obscure.
#
# Fortunately, various packages in Julia provide nicer interfaces to access it. We will use
# the one provided by `Flux.jl`.
#
# The data are images of handwritten digits, and the corresponding labels that were
# determined by hand (i.e. by humans). Our job will be to get the computer to **learn** to
# recognize digits by learning, as usual, the function that relates the input and output
# data.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Loading and examining the data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# First we load the required packages:
# ------------------------------------------------------------------------------------------

using Flux, Flux.Data.MNIST

# ------------------------------------------------------------------------------------------
# Now we read in the data:
# ------------------------------------------------------------------------------------------

labels = MNIST.labels();
images = MNIST.images();  # the semi-colon (`;`) here is important: it prevents Julia from displaying the object

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Examine the `labels` data. Then examine the first few images. *Do not try to view the
# whole of the `images` object!* Try to drill down to discover how the data is laid out.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Convert the first image to a matrix of `Float64`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Munging the data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In the previous notebooks, we arranged the input data for Flux as a `Vector` of `Vector`s.
# Now we will use an alternative arrangement, as a matrix, since that allows `Flux` to use
# matrix operations, which are more efficient.
#
# The column $i$ of the matrix is a vector consisting of the $i$th data point
# $\mathbf{x}^{(i)}$.  Similarly, the desired outputs are given as a matrix, with the $i$th
# column being the desired output $\mathbf{y}^{(i)}$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# An image is a matrix of colours, but now we need a vector instead. To do so, we just
# arrange all of the elements of the matrix in a certain way into a single list;
# fortunately, Julia already provides the function `vec` to do so!
#
# 1. Which order does `vec` use? [This reflects the underlying way in which the matrix is
# stored in memory.]
#
# 2. How can you convert an image into a `Vector` of `Float64`?
#
# 3. Define a variable $n$ that is the length of these vectors.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
# Make a function `rewrite` that accepts a range and converts that range of images to
# floating-point vectors and stacks them horizontally using `hcat` and the "splat" operator
# `...`.
#
# We also want a matrix of one-hot vectors. `Flux` provides a function `onehotbatch` to do
# this (you will need to import it). It works like `onehot`, but takes in a vector of labels
# and outputs a matrix `Y`.
#
# Return the pair `(X, Y)`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Setting up the neural network
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now we must set up a neural network. Since the data is complicated, we may expect to need
# several layers.
# But we can start with a single layer.
#
# - The network will take as inputs the vectors $\mathbf{x}^{(i)}$, so the input layer has
# $n$ nodes.
#
# - The output will be a one-hot vector encoding the digit from 1 to 9 or 0 that is desired.
# There are 10 possible categories, so we need an output layer of size 10.
#
# It is then our task as neural network designers to insert layers between these input and
# output layers, whose weights will be tuned during the learning process. *This is an art,
# not a science*! But major advances have come from finding a good structure for the
# network.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Softmax
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We will make a network with a single layer; let's choose each neuron in the layer to use
# the `relu` activation function.
# The output `relu` can be arbitrarily large, but in the end we will wish to compare the
# network's output with one-hot vectors, i.e. values between $0$ and $1$.
#
# In order to make this work, we will thus use an extra function at the end that takes in a
# vector of arbitrary real numbers and maps it ("squashes it down") to a vector of numbers
# between $0$ and $1$.
#
# The most used function with this property is $\mathrm{softmax}$. Firstly we take the
# exponential of each input variable to make them positive. Then we divide by the sum to
# make sure they lie between $0$ and $1$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# $$\mathrm{softmax}(\mathbf{x})_i := \frac{\exp (x_i)}{\sum_j \exp(x_j)}$$
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Note that here we have written the result for the $i$th component of the function
# $\mathbf{R}^n \to \mathbf{R}^n$. Note also that the function returns a vector of numbers
# that are positive, and whose components sum to $1$. Thus, in fact, they can be thought of
# as probabilities.
#
# In the neural network context, using a `softmax` after the final layer thus allows us to
# interpret the outputs as probabilities, in our case the probability that the network
# assigns that a given image represents each possible output value ($0$-$9$)!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# Make a neural network with one single layer, using the function $\sigma$, and a final
# `softmax`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Training
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# As we know, **training** consists of iteratively adjusting the model's parameters to
# decrease the `loss` function. Which parameters need to be adjusted? All of them!
#
# Since the `loss` function contains a call to the `model` function, calling `back!` on the
# result of the loss function updates the information about the gradient of the loss
# function with respect to *every node in the network!*:
# ------------------------------------------------------------------------------------------

l = loss(X, Y)

Flux.Tracker.back!(l)

# ------------------------------------------------------------------------------------------
# This is what is going on inside the `train!` function.
# In fact, `train!(loss, data, opt)` iterates over each object in `data` and runs this
# function.
# For this reason, `data` must consist of an iterable object that returns pairs `(X, Y)` at
# each step.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The simplest possibility is
# ------------------------------------------------------------------------------------------

data = ((X, Y), )  # one-element tuple

# ------------------------------------------------------------------------------------------
# Alternatively, we can make one call to the `train!` function iterate over several copies
# of `data`, using `repeated`. This is an **iterator**; it does not copy the data 100 times,
# which would be very wasteful; it just gives an object that repeatedly loops over the same
# data:
# ------------------------------------------------------------------------------------------

dataset = Base.Iterators.repeated((X, Y), 100)

# ------------------------------------------------------------------------------------------
# #### Exercise 6
#
# Train the model on a subset of $N$ images with $N = 5000$.
# ------------------------------------------------------------------------------------------

N = 5_000
X, Y = rewrite(1:N)

# ------------------------------------------------------------------------------------------
# The function `loss` evaluated on the matrices gives the overall loss:
# ------------------------------------------------------------------------------------------

loss(X, Y)

@time Flux.train!(loss, data, opt)

@time Flux.train!(loss, dataset, opt)

# ------------------------------------------------------------------------------------------
# This is (approximately) equivalent to just doing a `for` loop to run the previous `train!`
# command 100 times.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Using callbacks
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The `train!` function can take an optional keyword argument, `cb` (short for
# "*c*all*b*ack"). A callback function is a function that you provide as an argument to a
# function `f`, which "calls back" your function every so often.
#
# This provides the possibility to provide a function that is called at each step or every
# so often during the training process.
# A common use case is to provide a visual trace of the training process by printing out the
# current value of the `loss` function:
# ------------------------------------------------------------------------------------------

callback() = @show(loss(X, Y))

Flux.train!(loss, data, opt; cb = callback)

Flux.train!(loss, dataset, opt; cb = callback)

# ------------------------------------------------------------------------------------------
# However, it is expensive to calculate the complete `loss` function and it is not necessary
# to output it every step. So `Flux` also provides a function `throttle`, that provides a
# mechanism to call a given function at most once every certain number of seconds:
# ------------------------------------------------------------------------------------------

Flux.train!(loss, dataset, opt; cb = Flux.throttle(callback, 1))

for i in 1:100
    Flux.train!(loss, dataset, opt; cb = Flux.throttle(callback, 1))
end

# ------------------------------------------------------------------------------------------
# ## Testing phase
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We now have trained a model, i.e. we have found the parameters `W` and `b` for the network
# layer(s). In order to **test** if the learning procedure was really successful, we check
# how well the resulting trained network performs when we test it with images that the
# network has not yet seen!
#
# Often, a dataset is split up into "training data" and "test (or validation) data" for this
# purpose, and indeed the MNIST data set has a separate pool of training data. We can
# instead use the images that we have not included in our reduced training process.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 7
#
# Take the next 100 images after those that were used for training. How well does it do?
# ------------------------------------------------------------------------------------------

X_test, Y_test = rewrite(N+1:N+100)

loss(X_test, Y_test)

display(images[N+1])
labels[N+1]

[model(X_test[:,1]) Y_test[:,1]]

loss(X_test[:,1], Y_test[:,1])

loss(X_test, Y_test)

# ------------------------------------------------------------------------------------------
# #### Exercise 8
#
# Use the `indmax` function to write a function `prediction` that reports which digit
# `model` predicts, as the index with the maximum weight.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 9
#
# Count the number of correct predictions over the whole data set, and hence the percentage
# of images that are correctly predicted. [This percentage is what is used to compare
# different machine learning techniques.]
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Improving the prediction
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# So far we have used a single layer. In order to improve the prediction, we probably need
# to use more layers.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 10
#
# Introduce an intermediate, hidden layer. Does it give a better prediction?
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/19. Recognizing handwriting using a neural network.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module TestGP
using GaussianProcesses, ScikitLearnBase
using GaussianProcesses: set_params!, get_params
using Test, Random
using LinearAlgebra: diag

Random.seed!(1)

@testset "GP" begin
    d, n = 3, 10

    X = 2π * rand(d, n)
    y = [sum(sin, view(X, :, i)) / d for i in 1:n]
    mZero = MeanZero()
    kern = SE(0.0, 0.0)

    ntest = 5
    Xtest = randn(d, ntest)

    @testset "GPE constructors" begin
        gp = GP(X, y, mZero, kern)
        gp = GPE(X, y, mZero, kern)
        gp = GPE(X, y, mZero, kern, 1.2)
        gp = GPE(X, y, mZero, kern, GaussianProcesses.Scalar(1.2))
    end

    gp = GP(X, y, mZero, kern)

    # Verify that predictive mean at input observations
    # are the same as the output observations
    @testset "Predictive mean at obs locations" begin
        y_pred, σ2 = predict_y(gp, X)
        @test maximum(abs, gp.y - y_pred) ≈ 0.0 atol=0.1
        y_pred, pred_cov = predict_y(gp, X; full_cov=true)
        @test maximum(abs, gp.y - y_pred) ≈ 0.0 atol=0.1
        @test σ2 ≈ diag(pred_cov)
    end
    @testset "Predictive mean at test locations" begin
        y_pred, sig = predict_y(gp, Xtest)
        y_pred, sig = predict_y(gp, Xtest; full_cov=true)
    end

    # ScikitLearn interface test
    @testset "ScikitLearn interface" begin
        gp_sk = ScikitLearnBase.fit!(GPE(), X', y)
        y_pred = ScikitLearnBase.predict(gp_sk, X')
        @test maximum(abs, gp_sk.y - y_pred) ≈ 0.0 atol=0.1
    end

    # Modify kernel and update
    @testset "Update" begin
        gp.kernel.ℓ2 = 4.0
        X_pred = 2π * rand(d, n)

        GaussianProcesses.update_target!(gp)
        y_pred, sig = predict_y(gp, X_pred)
    end

    #Check that the rand function works
    @testset "Random GP sampling" begin
        X_test = 2π * rand(d, n)
        samples = rand(gp, X_test)
    end

    @testset "params round trip" begin
        params_1 = deepcopy(get_params(gp))
        set_params!(gp, params_1)
        params_2 = get_params(gp)
        @test params_1 ≈ params_2
    end

end
end


text:
ENV["JULIA_DEBUG"] = "DynaPlanning"

include("../src/DynaPlanning.jl")
include("../src/DynaUtils.jl")
include("TestDynaPlanning_Instances.jl")

using LightGraphs, LightGraphsFlows, GraphPlot, Compose # to reduce scope
using SparseArrays                      # idem
using DataStructures                    # idem
using Random
using Test
using Logging
using .Threads
using Actors
import Actors: spawn

import .DynaPlanning
import .DynaUtils
import .TestDynaPlanning_Instances

# Logging
logger = LogLevel(Logging.Debug)

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

@testset "Constructors" begin
    @test typeof(DynaPlanning.EnvironmentModel()) == Main.DynaPlanning.DictSrv{Actors.Link{Channel{Any}}}

    @test typeof(DynaPlanning.Incentiver()) == Main.DynaPlanning.Incentiver
    @test typeof(
        DynaPlanning.Incentiver(Main.DynaPlanning.StackelbergStrategy()).strategy_followed,
    ) == Main.DynaPlanning.StackelbergStrategy
    @test typeof(
        DynaPlanning.Incentiver(Main.DynaPlanning.SelfishStrategy()).strategy_followed,
    ) == Main.DynaPlanning.SelfishStrategy

    @test typeof(DynaPlanning.Q()) == Main.DynaPlanning.Q

    @test typeof(DynaPlanning.DynaAgentPlus()) == Main.DynaPlanning.DynaAgentPlus

    @test typeof(DynaPlanning.StackelbergInstance()) ==
          Main.DynaPlanning.StackelbergInstance
end

@testset "DynaUtils" begin
    ρ = DynaUtils.price_of_anarchy([0.0, 0.0, 1.0, 0.0], TestDynaPlanning_Instances.c_e)
    @test ρ == 4 / 3
end

@testset "Stackelberg" begin
    @test nthreads() == TestDynaPlanning_Instances.N_DynaAgents

    # find_flow in first epoch
    f_e_test = DynaUtils.find_flow(TestDynaPlanning_Instances.testStackelbergInstance)

    @test maximum(f_e_test) == 0.05
end

max_stars_repo_name: mgoar/dynaplanning

max_issues_repo_name: mgoar/dynaplanning

max_forks_repo_name: mgoar/dynaplanning


text:
using Revise
using IsopycnalSurfaces, Test

@testset "IsopycnalSurfaces.jl" begin
    #############################
    # Test dedup! function
    @testset "column" begin
        ################################
        # Idealized mapping onto sigma1
        pz = collect(0.:500.:4000.) # pressure levels
        nz = length(pz)

        θz = collect(range(20,stop=10,length=nz))
        Sz = collect(range(36,stop=35,length=nz))
        ztest = sort(rand(2:8,2))
        σ₁true = sigma1column(θz[ztest],Sz[ztest],pz[ztest])

        σ₁grid = range(minimum(σ₁true),stop=maximum(σ₁true),length=20)
        σ₁=sigma1column(θz,Sz,pz)

        @testset "column_spline" begin
            splorder = 3
            sgood = findall(minimum(σ₁) .<= σ₁grid .<= maximum(σ₁))
            pσ = var2sigmacolumn(σ₁,pz,σ₁grid[sgood],splorder)

            @test isapprox(pσ[begin],pz[ztest[begin]])
            @test isapprox(pσ[end],pz[ztest[end]])
        end

        @testset "column_linear" begin
            splorder = 100
            sgood = findall(minimum(σ₁) .<= σ₁grid .<= maximum(σ₁))
            pσ = var2sigmacolumn(σ₁,pz,σ₁grid[sgood],splorder)

            @test isapprox(pσ[begin],pz[ztest[begin]])
            @test isapprox(pσ[end],pz[ztest[end]])
        end

    end
end


text:
module visualization

importall MicroLogging
import ImageCore: colorview, normedview
import Colors: RGB, distinguishable_colors
import PlotUtils: cgrad

"""
Apply a colormap to an image.

If vmin or vmax is provided, normalize the values in the given range. `cmap` can
either be a matplotlib colormap name or a number. In the latter case, a colormap
is generated with `cmap` discrete, as well as possible distinguishable colors.
`cmap_seed` can be used to set a different seed color for creating these values.
"""
function apply_colormap{T<:Real}(image::Array{T, 2};
                                 vmin=nothing, vmax=nothing, cmap=:viridis,
                                 cmap_seed=nothing)
    @debug "Input image shape: $(size(image))." *
        "Applying colormap with vmin: $(vmin), vmax: $(vmax) and " *
        "cmap $(cmap)."
    # Prepare for normalization.
    vis_image = convert(Array{Float32, 2}, image);
    @debug "Current min: $(minimum(vis_image)), " *
        "current max: $(maximum(vis_image))."
    # Normalization.
    if vmin != nothing
        imin = Float32(vmin);
        clamp!(vis_image, vmin, typemax(Float32));
    else
        imin = minimum(vis_image);
    end
    if vmax != nothing
        imax = Float32(vmax);
        clamp!(vis_image, typemin(Float32), vmax);
    else
        imax = maximum(vis_image);
    end
    @debug "After normalization: min: $(minimum(vis_image)), " *
           "max: $(maximum(vis_image))."
    vis_image -= imin;
    vis_image /= (imax - imin);
    # Visualization.
    result = colorview(RGB, normedview(zeros(UInt8,
                                             3,
                                             size(vis_image)[1],
                                             size(vis_image)[2])));
    @debug "Result size: $(size(result))."
    if isa(cmap, Integer)
        vis_image = round.(Int32, vis_image .* Float32(cmap - 1));
        if cmap_seed == nothing
            map = distinguishable_colors(cmap);
        else
            map = distinguishable_colors(cmap, cmap_seed);
        end
        for row_idx = 1:size(vis_image)[2]
            for col_idx = 1:size(vis_image)[1]
                @inbounds result[col_idx, row_idx] = map[vis_image[col_idx,
                                                                   row_idx] + 1];
            end
        end
    else
        gradient = cgrad(cmap);
        for row_idx = 1:size(vis_image)[2]
            for col_idx = 1:size(vis_image)[1]
                @inbounds result[col_idx, row_idx] = gradient[
                    vis_image[col_idx, row_idx]];
            end
        end
    end
    return result;
end
end  # module


text:
module Controlz

using Polynomials
using DifferentialEquations
using PyPlot
using Interpolations
using Printf
using Roots
using DataFrames

# hipster plot theme
PyPlot.matplotlib.style.use(normpath(joinpath(pathof(Controlz), "..", "hipster.mplstyle")))

include("tf.jl")
include("special_tfs.jl")
include("systems.jl")
include("sim.jl")
include("viz.jl")
include("controls.jl")
include("margins.jl")
include("closed_loops.jl")
include("show.jl")

export
    # tf.jl
    TransferFunction, zeros_poles_gain, zeros_poles_k, zero_frequency_gain, evaluate, s, proper, strictly_proper, pole_zero_cancellation, zpk_form, system_order,
    # special_tfs.jl
    first_order_system, second_order_system, time_constant, damping_coefficient,
    # systems.jl
    characteristic_polynomial,
    # show.jl
    # sim.jl
    simulate, interpolate,
    # viz.jl
    viz_response, viz_poles_and_zeros, nyquist_diagram, root_locus, bode_plot, mk_gif,
    # controlz.jl
    PController, PIController, PIDController,
    # margins.jl
    gain_phase_margins,
    # closedloops.jl
    ClosedLoopTransferFunction
end


text:
using Base: BitSignedSmall
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    MinDistanceSampling(α, ρ=0.65, δ=100, metric=Euclidean())

Generate sample from geometric object such that all pairs of
points are at least `α` units of distance away from each other.
Optionally specify the relative radius `ρ` for the packing
pattern, the oversampling factor `δ` and the `metric`.

This method is sometimes referred to as Poisson disk sampling
or blue noise sampling in the computer graphics community.

## References

* Lagae, A. & Dutré, P. 2007. [A Comparison of Methods for
  Generating Poisson Disk Distributions]
  (https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-8659.2007.01100.x)
* Bowers et al. 2010. [Parallel Poisson disk sampling with
  spectrum analysis on surfaces](https://dl.acm.org/doi/10.1145/1882261.1866188)
* Medeiros et al. 2014. [Fast adaptive blue noise on polygonal surfaces]
  (https://www.sciencedirect.com/science/article/abs/pii/S1524070313000313)
"""
struct MinDistanceSampling{T,M} <: ContinuousSamplingMethod
  α::T
  ρ::T
  δ::Int
  metric::M
end

MinDistanceSampling(α::T; ρ=T(0.65), δ=100, metric=Euclidean()) where {T} =
  MinDistanceSampling(α, ρ, δ, metric)

function sample(rng::AbstractRNG, Ω::DomainOrData, method::MinDistanceSampling)
  # retrive parameters
  α = method.α
  ρ = method.ρ
  δ = method.δ
  m = method.metric

  # total volume/area of the object
  V = sum(measure, Ω)

  # expected number of Poisson samples
  # for relative radius (Lagae & Dutré 2007)
  N = 2V/√3 * (ρ/α)^2

  # number of oversamples (Medeiros et al. 2014)
  O = ceil(Int, δ * N)

  # oversample the object
  points = sample(rng, Ω, HomogeneousSampling(O))

  # collect points into point set
  𝒫 = PointSet(collect(points))

  # discard points that do not satisfy distance criterion
  sample(rng, 𝒫, BallSampling(α, metric=m))
end


text:
export COFFDynamicLinks, COFFDynamicLink, COFFRPath

struct COFFDynamicLink{H <: COFFHandle} <: DynamicLink{H}
    path::String
end
path(dl::COFFDynamicLink) = dl.path

@io struct COFFImageImportDescriptor
    Characteristics::UInt32
    TimeDateStamp::UInt32
    ForwarderChain::UInt32
    Name::UInt32
    FirstThunk::UInt32
end

struct COFFDynamicLinks{H <: COFFHandle} <: DynamicLinks{H}
    handle::H
    links::Vector
end

function import_name(oh::COFFHandle, idata::COFFSectionRef, iid)
    seek(oh, iid.Name - section_address(idata) + section_offset(idata))
    return strip(readuntil(oh, '\0'), '\0')
end

function find_section_for_rva(oh::H, rva) where {H <: COFFHandle}
    for s in Sections(oh)
        if section_address(s) <= rva && section_address(s) + section_size(s) > rva
            return s
        end
    end
    error("Unable to find section for RVA $(repr(rva))")
end

function COFFDynamicLinks(oh::H) where {H <: COFFHandle}
    # Start by finding the virtual address of the import table
    import_table_rva = oh.opt_header.directories.ImportTable.VirtualAddress

    # Next figure out which section that belongs to:
    s = find_section_for_rva(oh, import_table_rva)

    # We'll load in all the ImageImportDescriptors we can
    iids = COFFImageImportDescriptor[]

    # Read in ImageImportDescriptors until it's all NULL
    seek(oh, (import_table_rva - section_address(s)) +  section_offset(s))
    while true
        iid = unpack(oh, COFFImageImportDescriptor)
        if iid.Name == 0 && iid.FirstThunk == 0 && iid.Characteristics == 0 && iid.ForwarderChain == 0 && iid.TimeDateStamp == 0
            break
        else
            push!(iids, iid)
        end
    end

    # Now, jump around and get all the strings
    links = [COFFDynamicLink{H}(import_name(oh, s, iid)) for iid in iids]
    return COFFDynamicLinks(oh, links)
end
DynamicLinks(oh::COFFHandle) = COFFDynamicLinks(oh)

handle(dls::COFFDynamicLinks) = dls.handle
lastindex(dls::COFFDynamicLinks) = lastindex(dls.links)
getindex(dls::COFFDynamicLinks, idx) = getindex(dls.links, idx)

"""
    COFFRPath

COFF RPath object; note that while COFF files do not have an RPath within them,
they _do_ seach the same directory as the loading binary (e.g. the `\$ORIGIN`).
We use `COFFRPath` to effect this, although strictly speaking there is no such
thing as a "COFF RPath".
"""
struct COFFRPath{H <: COFFHandle} <: RPath{H}
    handle::H
end

RPath(oh::COFFHandle) = COFFRPath(oh)
handle(crp::COFFRPath) = crp.handle
rpaths(crp::COFFRPath) = String[]

# COFF files don't have an RPath, but they _do_ always search the $ORIGIN
function canonical_rpaths(crp::COFFRPath)
    return [abspath(dirname(path(handle(crp))) * "/")]
end


text:
"""
    Measurements

"""
abstract type Measurement end

"""
    Mechanisms

"""

abstract type AbstractMechanism end

"""
    Species and Species Pools

"""
abstract type AbstractSpecies end
abstract type DiscreteSpeciesPool end
abstract type ContinuousSpeciesPool end

"""
    Landscapes

"""

"""
    Tensors

"""

"""
    Summary stats

"""


text:
module EmpiricalModeDecomposition

using AbstractFFTs
using Dierckx
using Distributions
using FFTW
using ForwardDiff
using LinearAlgebra
using Random
using Statistics

include("mode-decompositions/emd.jl")
include("mode-decompositions/eemd.jl")
include("extras/hilbert.jl")

export EMDSetting, emd, EEMDSetting, eemd, hilbert_transform, hht

end # module


text:
using IterativeSolvers
using Base.Test

import Base: size, A_mul_B!, Ac_mul_B!, eltype, similar, scale!, copy!, fill!, length, broadcast!
import Base.LinAlg: norm

# Type used in Dampenedtest
# solve (A'A + diag(v).^2 ) x = b
# using LSMR in the augmented space A' = [A ; diag(v)] b' = [b; zeros(size(A, 2)]
mutable struct DampenedVector{Ty, Tx}
    y::Ty
    x::Tx
end

eltype(a::DampenedVector) = promote_type(eltype(a.y), eltype(a.x))
norm(a::DampenedVector) = sqrt(norm(a.y)^2 + norm(a.x)^2)

function Base.Broadcast.broadcast!(f::Tf, to::DampenedVector, from::DampenedVector, args...) where {Tf}
    to.x .= f.(from.x, args...)
    to.y .= f.(from.y, args...)
    to
end

function copy!(a::DampenedVector{Ty, Tx}, b::DampenedVector{Ty, Tx}) where {Ty, Tx}
    copy!(a.y, b.y)
    copy!(a.x, b.x)
    a
end

function fill!(a::DampenedVector, α::Number)
    fill!(a.y, α)
    fill!(a.x, α)
    a
end

scale!(α::Number, a::DampenedVector) = scale!(a, α)

function scale!(a::DampenedVector, α::Number)
    scale!(a.y, α)
    scale!(a.x, α)
    a
end

similar(a::DampenedVector, T) = DampenedVector(similar(a.y, T), similar(a.x, T))
length(a::DampenedVector) = length(a.y) + length(a.x)

mutable struct DampenedMatrix{TA, Tx}
    A::TA
    diagonal::Tx
end

eltype(A::DampenedMatrix) = promote_type(eltype(A.A), eltype(A.diagonal))

function size(A::DampenedMatrix)
    m, n = size(A.A)
    l = length(A.diagonal)
    (m + l, n)
end

function size(A::DampenedMatrix, dim::Integer)
    m, n = size(A.A)
    l = length(A.diagonal)
    dim == 1 ? (m + l) :
    dim == 2 ? n : 1
end

function A_mul_B!(α::Number, mw::DampenedMatrix{TA, Tx}, a::Tx,
    β::Number, b::DampenedVector{Ty, Tx}) where {TA, Tx, Ty}
    if β != 1.
        if β == 0.
            fill!(b, 0.)
        else
            scale!(b, β)
        end
    end
    A_mul_B!(α, mw.A, a, 1.0, b.y)
    map!((z, x, y)-> z + α * x * y, b.x, b.x, a, mw.diagonal)
    return b
end

function Ac_mul_B!(α::Number, mw::DampenedMatrix{TA, Tx}, a::DampenedVector{Ty, Tx},
    β::Number, b::Tx) where {TA, Tx, Ty}
    if β != 1.
        if β == 0.
            fill!(b, 0.)
        else
            scale!(b, β)
        end
    end
    Ac_mul_B!(α, mw.A, a.y, 1.0, b)
    map!((z, x, y)-> z + α * x * y, b, b, a.x, mw.diagonal)
    return b
end

"""
Produces the m × n submatrix from
A = [ 1
      1 2
        2 3
          3 4
            ...
              n ]
suitably padded by zeros.
"""
function sol_matrix(m, n)
    mn = min(m, n)
    spdiagm((1.0 : mn - 1, 1.0 : mn), (-1, 0), m, n)
end

@testset "LSMR" begin
    srand(1234321)

    @testset "Small dense matrix" for T = (Float32, Float64)
        A = rand(T, 10, 5)
        b = rand(T, 10)
        x, history = lsmr(A, b, log = true)
        @test isa(history, ConvergenceHistory)
        @test norm(x - A\b) ≤ √eps(T)
    end

    @testset "SOL test" for (m, n, damp) = ((10, 10, 0), (20, 10, 0), (20, 10, 0.1))
        # Test adapted from the BSD-licensed Matlab implementation at
        #    http://www.stanford.edu/group/SOL/software/lsqr.html
        #              Michael Saunders, Systems Optimization Laboratory,
        #              Dept of MS&E, Stanford University.
        #-----------------------------------------------------------------------
        # 11 Apr 1996: First version for distribution with lsqr.m.
        #              Michael Saunders, Dept of EESOR, Stanford University.

        A = sol_matrix(m, n)
        x = float(n : -1 : 1)
        b = A * x
        x_lsmr = lsmr(A, b, atol = 1e-7, btol = 1e-7, conlim = 1e10, maxiter = 10n)
        @test norm(b - A * x) ≤ 1e-4
    end

    @testset "Dampened test" for (m, n) = ((10, 10), (20, 10))
        # Test used to make sure A, b can be generic matrix / vector
        b = rand(m)
        A = rand(m, n)
        v = rand(n)
        Adampened = DampenedMatrix(A, v)
        bdampened = DampenedVector(b, zeros(n))
        x, ch = lsmr(Adampened, bdampened, log=true)
        @test norm((A'A + diagm(v) .^ 2)x - A'b) ≤ 1e-3
    end
end


text:
## Matrix whole space
function sample_network!(res::MCMCGradientEvaluationCache,
  prob::HermitianMatrixProblem, net, σ, wholespace=false)

  lnψ, ∇lnψ = logψ_and_∇logψ!(res.∇lnψ, net, σ)
  E         = compute_Cloc(prob, net, σ, lnψ, res.σ)

  prob = wholespace ? exp(2*real(lnψ)) : 1.0
  res.Zave   += prob
  res.Eave   += prob  * E
  push!(res.Evalues, prob*E)

  for (i, _∇lnψ)= enumerate(∇lnψ.tuple_all_weights)
    res.Oave[i]  .+= prob .* _∇lnψ
    res.EOave[i] .+= prob .* E .* conj.(_∇lnψ)
  end
  return res
end

## Matrix whole space
function sample_network!(res::MCMCGradientLEvaluationCache, prob::LRhoSquaredProblem,
                         net, σ, wholespace=false)
  CLO_i = res.LLO_i

  lnψ, ∇lnψ = logψ_and_∇logψ!(res.∇lnψ, net, σ)
  C_loc = compute_Cloc!(CLO_i, res.∇lnψ2, prob, net, σ, lnψ, res.σ)

  prob = wholespace ? exp(2*real(lnψ)) : 1.0
  E = abs(C_loc)^2

  res.Zave += prob
  res.Eave += prob * E
  push!(res.Evalues, prob*E)

  for (i, _∇lnψ)= enumerate(∇lnψ.tuple_all_weights)
    res.Oave[i]  .+= prob .* _∇lnψ
    res.EOave[i] .+= prob .* E .* conj.(_∇lnψ)
    res.LLOave[i].+= prob .* conj(C_loc) .* CLO_i[i]
  end
  return res
end


text:
using Pumas, Plots, CSV, Random
Random.seed!(0)

single_dose_regimen = DosageRegimen(100, time=0)
first(single_dose_regimen.data)

s1 = Subject(id=1, evs=single_dose_regimen,cvs=(Wt=70,))

choose_covariates() = (Wt = rand(55:80),)

pop = Population(map(i -> Subject(id = i,evs = single_dose_regimen, cvs =  choose_covariates()),1:24))

pop[1].covariates

mymodel = @model begin
  @param   begin
    tvcl ∈ RealDomain(lower=0, init = 1.0)
    tvv ∈ RealDomain(lower=0, init = 20)
    tvka ∈ RealDomain(lower = 0, init= 1)
    Ω ∈ PDiagDomain(init=[0.09,0.09, 0.09])
    σ_prop ∈ RealDomain(lower=0,init=0.04)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    CL = tvcl * (Wt/70)^0.75 * exp(η[1])
    V  = tvv * (Wt/70) * exp(η[2])
    Ka = tvka * exp(η[3])
  end
  @covariates Wt

  @dynamics Depots1Central1
    #@dynamics begin
    #    Depot' =  -Ka*Depot
    #    Central' =  Ka*Depot - (CL/V)*Central
    #end

  @derived begin
      cp = @. 1000*(Central / V)
      dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
    end
end

param = init_param(mymodel)

obs = simobs(mymodel, pop, param, obstimes=0:1:72)
plot(obs)

simdf = DataFrame(obs)
first(simdf, 6)

simdf[!, :route] .= "ev"

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg"

ncadf = read_nca(simdf, id=:id, time=:time, conc=:cp, amt=:amt,
    route=:route,timeu=timeu, concu=concu, amtu=amtu, lloq=0.4concu)

plot(ncadf)

auc = NCA.auc(ncadf)

cmax = NCA.cmax(ncadf)

report = NCAReport(ncadf)
report = NCA.to_dataframe(report)
first(report,6)

simdf.cmt = ifelse.(ismissing.(simdf.cmt), 2, simdf.cmt)
est_df = simdf[.!((simdf.dv .== 0.0) .& (simdf.cmt .==2)),:]
first(est_df,6)

data = read_pumas(est_df ,cvs = [:Wt], dvs=[:dv])

res = fit(mymodel,data,param,Pumas.FOCEI())

infer(res)

preds = DataFrame(predict(res))
first(preds, 6)

resids = DataFrame(wresiduals(res))
first(resids, 6)

ebes = DataFrame(empirical_bayes(res))
first(ebes, 6)

resout = DataFrame(inspect(res))
first(resout, 6)