text:
# Test code on a relatively large molecule; my old friend C60

# Geom: D. Dragoe
# http://www.ccl.net/chemistry/resources/messages/2000/12/07.005-dir/
C60 = Quante.Molecule([
Atom( 6 , 0.726656 , -1.000157 , 3.300459),
Atom( 6 , 1.175755 , 0.382026 , 3.300459),
Atom( 6 , 1.410183 , -1.940951 , 2.581756),
Atom( 6 , 2.281725 , 0.741377 , 2.581756),
Atom( 6 , 2.281725 , 1.977639 , 1.817704),
Atom( 6 , 0.000000 , 1.236262 , 3.300459),
Atom( 6 , 0.000000 , 2.399147 , 2.581756),
Atom( 6 , 1.175755 , 2.781173 , 1.817704),
Atom( 6 , -0.683527 , -2.941108 , 1.817704),
Atom( 6 , -1.410183 , -1.940951 , 2.581756),
Atom( 6 , 0.683527 , -2.941108 , 1.817704),
Atom( 6 , -0.726656 , -1.000157 , 3.300459),
Atom( 6 , -1.175755 , 0.382026 , 3.300459),
Atom( 6 , -2.585938 , -1.558925 , 1.817704),
Atom( 6 , -3.008381 , -0.258779 , 1.817704),
Atom( 6 , -2.281725 , 0.741377 , 2.581756),
Atom( 6 , 0.726656 , -3.399304 , -0.581443),
Atom( 6 , -0.726656 , -3.399304 , -0.581443),
Atom( 6 , 1.410183 , -3.177213 , 0.581443),
Atom( 6 , -1.410183 , -3.177213 , 0.581443),
Atom( 6 , -2.585938 , -2.322977 , 0.581443),
Atom( 6 , -1.175755 , -2.781173 , -1.817704),
Atom( 6 , -2.281725 , -1.977639 , -1.817704),
Atom( 6 , -3.008381 , -1.741534 , -0.581443),
Atom( 6 , 3.008381 , -1.741534 , -0.581443),
Atom( 6 , 2.281725 , -1.977639 , -1.817704),
Atom( 6 , 2.585938 , -2.322977 , 0.581443),
Atom( 6 , 1.175755 , -2.781173 , -1.817704),
Atom( 6 , 0.000000 , -2.399147 , -2.581756),
Atom( 6 , 2.281725 , -0.741377 , -2.581756),
Atom( 6 , 1.175755 , -0.382026 , -3.300459),
Atom( 6 , 0.000000 , -1.236262 , -3.300459),
Atom( 6 , 3.008381 , -0.258779 , 1.817704),
Atom( 6 , 3.457479 , 0.359351 , 0.581443),
Atom( 6 , 2.585938 , -1.558925 , 1.817704),
Atom( 6 , 3.457479 , -0.359351 , -0.581443),
Atom( 6 , 3.008381 , 0.258779 , -1.817704),
Atom( 6 , 3.008381 , 1.741534 , 0.581443),
Atom( 6 , 2.585938 , 2.322977 , -0.581443),
Atom( 6 , 2.585938 , 1.558925 , -1.817704),
Atom( 6 , -0.726656 , 1.000157 , -3.300459),
Atom( 6 , 0.726656 , 1.000157 , -3.300459),
Atom( 6 , 1.410183 , 1.940951 , -2.581756),
Atom( 6 , -1.410183 , 1.940951 , -2.581756),
Atom( 6 , -3.008381 , 0.258779 , -1.817704),
Atom( 6 , -2.281725 , -0.741377 , -2.581756),
Atom( 6 , -1.175755 , -0.382026 , -3.300459),
Atom( 6 , -2.585938 , 1.558925 , -1.817704),
Atom( 6 , -3.008381 , 1.741534 , 0.581443),
Atom( 6 , -3.457479 , 0.359351 , 0.581443),
Atom( 6 , -3.457479 , -0.359351 , -0.581443),
Atom( 6 , -2.585938 , 2.322977 , -0.581443),
Atom( 6 , -0.726656 , 3.399304 , 0.581443),
Atom( 6 , -1.175755 , 2.781173 , 1.817704),
Atom( 6 , -2.281725 , 1.977639 , 1.817704),
Atom( 6 , -1.410183 , 3.177213 , -0.581443),
Atom( 6 , 0.683527 , 2.941108 , -1.817704),
Atom( 6 , 1.410183 , 3.177213 , -0.581443),
Atom( 6 , 0.726656 , 3.399304 , 0.581443),
Atom( 6 , -0.683527 , 2.941108 , -1.817704)
])

AtoBohr!(C60) # convert from Angstrom to Bohr

@time Energy,E,U = rhf(C60,verbose=true)
println("C60 Energy: $Energy")


text:
##
#
# Generate all non-zero roots of polynomials with coefficients 1 or 3 of
# degree at most 22.
#
# Warning: this may take a while (a few hours in my case) and generate almost
#          5.6 GB file (352 321 541 roots).
#

using ChristensenSet

polynomials = PolynomialIterator([1.0, 3.0], 22)

@time find_roots!(polynomials, "ex13_Float64.dat")

ext: jl


text:
function collect_branches( subsys, nodes )  #=

collects the branches in a given circuit. It is assumed that the information of the circuit is contained in the variable "subsys" given as the first argument of this function. The variable "nodes" given as the second argument of this function is an array of the independent nodes in the circuit. The independent nodes may have been collected from the circuit description using the function collect_nodes(). The branches collected by this function will be returned as the array "branches" of the struct Branch. The number "num_branches" of branches (sum of branch phases) and the number "num_dyns" of dynamic branches will also be returned. In the array "branches", the node names have been replaced by their index numbers according to the array "nodes".

Author: Taku Noda
Started on Nov. 11, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

# -- Internally used function ------------------------------------------

function find( node, nodes )
index = 0
for i = 1:length( nodes )
  if nodes[ i ] == node
    index = i
    break
  end
end
return index
end

# -- Main body starts here ---------------------------------------------

branches = Branch[]
num_branches = 0
num_dyns = 0
num_nls = 0

for elm in eachelement( subsys )
  if elm.name == "branch"
    name = ""; type = ""; nph = 0; val_str = ""
    lhs = Int[]; rhs = Int[]; p_vars = 0
    p_de = 0
    p_ns = 0; plx = []; ply = [];
    #
    type = elm[ "type" ]
    name = elm[ "name" ]
    #
    if !haskey( elm, "nph" )  # if single-phase,
      nph = 1
      if ( type == "R" )||( type == "L" )||( type == "G" )||( type == "C" )
        val_str = elm[ "val" ]
      elseif ( type == "NR" )||( type == "NL" )||( type == "NC" )||( type == "E" )||( type == "J" )||( type == "Sw" )
        val_str = elm.content
      end
    else  # if multiphase,
      nph = parse( Int, elm[ "nph" ] )
      for elm1 in eachelement( elm )
        if elm1.name == "val"
          val_str = elm1.content
          break
        end
      end
    end
    if ( type == "R" )||( type == "L" )||( type == "G" )||( type == "C" )||( type == "NR" )||( type == "NL" )||( type == "NC" )||( type == "Sw" )
      val_str = '['*replace( val_str, r"\t|\n" => "" )*']'
    elseif ( type == "E" )||( type == "J" )
      val_str = replace( val_str, r" |\t|\n" => "" )
    end
    val = Meta.parse( val_str );
    if ( type == "NR" )||( type == "NL" )||( type == "NC" )
      pldat = eval( val )
      plx = pldat[ :, 1 ]
      ply = pldat[ :, 2 ]
    elseif type == "Sw"
      plx = eval( val )
    end
    #
    if ( type == "L" )||( type == "C" )||( type == "NL" )||( type == "NC" )
      p_de = num_dyns + 1
      num_dyns = num_dyns + nph
    else
      p_de = 0
    end
    if ( type[ 1 ] == 'N' )||( type == "E" )||( type == "J" )||( type == "CSw" )||( type == "CE" )||( type == "CJ" )||( type == "VP" )||( type == "CP" )
      p_ns = num_branches + 1
    end
    if ( type[ 1 ] == 'N' )
      num_nls = num_nls + 1
      p_nl = num_nls
    end
    if ( ( type == "R" )||( type == "L" )||( type == "G" )||( type == "C" ) )&&( nph != 1 )
      for elm1 in eachelement( elm )
        nph = parse( Int, elm[ "nph" ] )
        if elm1.name == "lhs"
          node_names = split( elm1.content, ";" )
          if length( node_names ) != nph
            omegasim_error( "For the branch $name, the number of lefthand-side node names specified does not match the number of phases." )
          end
          for node_name in node_names
            push!( lhs, find( node_name, nodes ) )
          end
        elseif elm1.name == "rhs"
          node_names = split( elm1.content, ";" )
          if length( node_names ) != nph
            omega_error( "For the branch $name, the number of righthand-side node names specified does not match the number of phases." )
          end
          for node_name in node_names
            push!( rhs, find( node_name, nodes ) )
          end
        end
      end
    else
      push!( lhs, find( elm[ "lhs" ], nodes ) )
      push!( rhs, find( elm[ "rhs" ], nodes ) )
    end
    #
    p_vars = num_branches + 1
    num_branches = num_branches + nph
    #
    branch = Branch( name, type, nph, val, lhs, rhs, p_vars, p_de, p_ns, plx, ply, [ -1 ] )
    push!( branches, branch )
  end
end

return branches, num_branches, num_dyns

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
"""
Uses a Levenberg-Marquardt solver to compute the best fit of a friction model to
a predefined data set.  As the supplied data is discretely sampled but the
solver may require points between those supplied, the solver will utilize linear
interpolation to estimate values between supplied data points.

The routine returns the fitted model along with the LsqFitResults type returned
from the Levenberg-Marquardt solver.  The LsqFitResults allow exploration of
error margins, confidence intervals, etc.  Both results are returned as a
named tuple with the fitted model available as 'model' and the LsqFitResults
available as 'fit'.

Control over the ODE solver tolerances and maximum allowable time step are
exposed to the calling code via the varargs 'pargs' input parameter.  The
following inputs are currently recognized.
- reltol: The relative ODE solver tolerance. The defaults is 1e-8.
- abstol: The absolute ODE solver tolerance.  The default is 1e-6.
- dtmax: The maximum allowable step size.  The default is 1e-3.
- lower: An N-element array of lower bounds on each of the N coefficients.
- upper: An N-element array of upper bounds on each of the N coefficients.
"""
function fit_model(
    mdl::HeuristicFrictionModel,
    timedata::Array{T},
    frictiondata::Array{T},
    normaldata::Array{T},
    positiondata::Array{T},
    velocitydata::Array{T},
    z0::Array{T};
    pargs...
) where T <: Number

    # Set up the interpolation for each data set
    normal_interp = LinearInterpolation(
        timedata, normaldata,
        extrapolation_bc = Line()
    )
    position_interp = LinearInterpolation(
        timedata, positiondata,
        extrapolation_bc = Line()
    )
    velocity_interp = LinearInterpolation(
        timedata, velocitydata,
        extrapolation_bc = Line()
    )

    # Define the model.
    # - t_: The time data.
    # - p_: An array of parameter values
    function model(t_, p_) # t_ is an array of time points
        # Define the model
        m = model_from_array(mdl, p_)

        # Define functions to compute the velocity and normal force at the
        # appropriate time.  This routine will utilize the interpolation
        # objects previously defined.
        pos(ti) = position_interp(ti)
        vel(ti) = velocity_interp(ti)
        nrm(ti) = normal_interp(ti)

        # Solve the model
        rsp = friction(m, t_, nrm, pos, vel, z0; pargs)
        return rsp.f
    end

    # Convert the model to array form
    p0 = model_to_array(mdl)

    # Extract the model options
    opt = extract_options(pargs, length(p0))

    # Fit the model
    fit = curve_fit(
        model,
        timedata,
        frictiondata,
        p0,
        lower = opt.lower,
        upper = opt.upper
    )
    p = coef(fit)

    return (model = model_from_array(mdl, p), fit = fit)
end

"""
Uses a Levenberg-Marquardt solver to compute the best fit of a friction model to
a predefined data set.  As the supplied data is discretely sampled but the
solver may require points between those supplied, the solver will utilize linear
interpolation to estimate values between supplied data points.  This instance
is used when the model does not depend upon position data.

The routine returns the fitted model along with the LsqFitResults type returned
from the Levenberg-Marquardt solver.  The LsqFitResults allow exploration of
error margins, confidence intervals, etc.  Both results are returned as a
named tuple with the fitted model available as 'model' and the LsqFitResults
available as 'fit'.

Control over the ODE solver tolerances and maximum allowable time step are
exposed to the calling code via the varargs 'pargs' input parameter.  The
following inputs are currently recognized.
- reltol: The relative ODE solver tolerance. The defaults is 1e-8.
- abstol: The absolute ODE solver tolerance.  The default is 1e-6.
- dtmax: The maximum allowable step size.  The default is 1e-3.
- lower: An N-element array of lower bounds on each of the N coefficients.
- upper: An N-element array of upper bounds on each of the N coefficients.
"""
function fit_model(
    mdl::HeuristicFrictionModel,
    timedata::Array{T},
    frictiondata::Array{T},
    normaldata::Array{T},
    velocitydata::Array{T},
    z0::Array{T};
    pargs...
) where T <: Number

    pos = zeros(T, length(timedata))
    return fit_model(mdl, timedata, frictiondata, normaldata, pos, velocitydata, z0; pargs...)
end

"""
Uses a Levenberg-Marquardt solver to compute the best fit of a friction model to
a predefined data set.  As the supplied data is discretely sampled but the
solver may require points between those supplied, the solver will utilize linear
interpolation to estimate values between supplied data points.

The routine returns the fitted model along with the LsqFitResults type returned
from the Levenberg-Marquardt solver.  The LsqFitResults allow exploration of
error margins, confidence intervals, etc.  Both results are returned as a
named tuple with the fitted model available as 'model' and the LsqFitResults
available as 'fit'.

Control over the ODE solver tolerances and maximum allowable time step are
exposed to the calling code via the varargs 'pargs' input parameter.  The
following inputs are currently recognized.
- reltol: The relative ODE solver tolerance. The defaults is 1e-8.
- abstol: The absolute ODE solver tolerance.  The default is 1e-6.
- dtmax: The maximum allowable step size.  The default is 1e-3.
- lower: An N-element array of lower bounds on each of the N coefficients.
- upper: An N-element array of upper bounds on each of the N coefficients.
"""
function fit_model(
    mdl::FrictionModel,
    timedata::Array{T},
    frictiondata::Array{T},
    normaldata::Array{T},
    velocitydata::Array{T};
    pargs...
) where T <: Number

    normal_interp = LinearInterpolation(
        timedata, normaldata,
        extrapolation_bc = Line()
    )
    velocity_interp = LinearInterpolation(
        timedata, velocitydata,
        extrapolation_bc = Line()
    )

    function model(t_, p_)
        m = model_from_array(mdl, p_)
        return friction.(m, normal_interp.(t_), velocity_interp.(t_))
    end

    p0 = model_to_array(mdl)
    opt = extract_options(pargs, length(p0))

    fit = curve_fit(
        model,
        timedata,
        frictiondata,
        p0,
        lower = opt.lower,
        upper = opt.upper
    )
    p = coef(fit)
    return (model = model_from_array(mdl, p), fit = fit)
end


text:
module Stan

using Compat, Documenter, DelimitedFiles, Statistics

"""
The directory which contains the CmdStan executables such as `bin/stanc` and
`bin/stansummary`. Inferred from `Main.CMDSTAN_HOME` or `ENV["CMDSTAN_HOME"]`
when available. Use `set_cmdstan_home!` to modify.
"""
global CMDSTAN_HOME=""

function __init__()
    global CMDSTAN_HOME = if isdefined(Main, :CMDSTAN_HOME)
        eval(Main, :CMDSTAN_HOME)
    elseif haskey(ENV, "CMDSTAN_HOME")
        ENV["CMDSTAN_HOME"]
    else
        @warn("Environment variable CMDSTAN_HOME not found. Use set_cmdstan_home!.")
        ""
    end
end

"""Set the path for the `CMDSTAN_HOME` environment variable.

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`
"""
set_cmdstan_home!(path) = global CMDSTAN_HOME=path

include("main/stanmodel.jl")
include("main/stancode.jl")
include("utilities/parallel.jl")
include("utilities/create_cmd_line.jl")
include("utilities/create_r_files.jl")
include("utilities/read_stan_files.jl")
include("types/sampletype.jl")
include("types/optimizetype.jl")
include("types/diagnosetype.jl")
include("types/variationaltype.jl")

export
# From this file
set_cmdstan_home!,
CMDSTAN_HOME,
# From stanmodel.jl
Stanmodel,
# From stancode.jl
stan,
# From sampletype.jl
Sample,
# From optimizetype.jl
Optimize,
# From diagnosetype.jl
Diagnose,
# From variationaltype.jl
Variational,
# From DelimitedFiles
writedlm,
readdlm

end # module


text:
import Base:in, +, -, *, ==, convert, isempty, isequal, length, push!, sizeof, append!
import SeisIO: GeoLoc, PZResp
export FFTData

const fftfields = [:name, :loc, :fs, :gain, :freqmin, :freqmax, :cc_step,
                   :cc_len, :whitened, :time_norm, :resp]

"""
   FFTData

A structure for fourier transforms (FFT) of ambient noise data.

## Fields: FFTData
| **Field** | **Description** |
|:--------------------|:--------------------|
| :name       | Freeform channel names |
| :id         | Channel ids. use NET.STA.LOC.CHAN format when possible. |
| :loc        | Location (position) object |
| :fs         | Sampling frequency in Hz. |
| :gain       | Scalar gain; divide data by the gain to convert to units  |
| :freqmin    | Minimum frequency for whitening.  |
| :freqmax    | Maximum frequency for whitening. |
| :cc_len     | Length of each correlation in seconds. |
| :cc_step    | Spacing between correlation windows in seconds. |
| :whitened   | Whitening applied.
| :time_norm  | Apply one-bit whitening with "one_bit". |
| :resp       | Instrument response object, format [zeros poles] |
| :misc       | Dictionary for non-critical information. |
| :notes      | Timestamped notes; includes automatically-logged acquisition and |
|             | processing information. |
| :t          | Starttime of each FFT. |
| :fft        | FFTs stored in columns. |
"""
mutable struct FFTData
  name::String                                # name [Net.Sta.Loc.Chan]
  id::String                                  # id [Y-mm-dd] this is date of fft
  loc::GeoLoc                                 # loc
  fs::Float64                                 # sampling rate [Hz]
  gain::Float64                               # gain
  freqmin::Float64                            # minumum frequency [Hz]
  freqmax::Float64                            # maximum frequency [Hz]
  cc_len::Int                                 # window_length [s]
  cc_step::Int                                # step between windows [s]
  whitened::Bool                              # whitening applied
  time_norm::Union{Bool,String}               # time normaliation
  resp::PZResp                         # response poles/zeros
  misc::Dict{String,Any}                      # misc
  notes::Array{String,1}                      # notes
  t::Array{Float64,1}                         # time
  fft::Array{<:Union{Complex{Float32},Complex{Float64}},2}  # fft data

  function FFTData(
      name     ::String,
      id       ::String,
      loc      ::GeoLoc,
      fs       ::Float64,
      gain     ::Float64,
      freqmin  ::Float64,
      freqmax  ::Float64,
      cc_len   ::Int,
      cc_step  ::Int,
      whitened ::Bool,
      time_norm::Union{Bool,String},
      resp     ::PZResp,
      misc     ::Dict{String,Any},
      notes    ::Array{String,1},
      t        ::Array{Float64,1},
      fft      ::Array{<:Union{Complex{Float32},Complex{Float64}},2}
      )

      return new(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
                 whitened, time_norm, resp, misc, notes, t, fft)
    end
end

FFTData(;
          name     ::String                    = "",
          id       ::String                    = "",
          loc      ::GeoLoc                    = GeoLoc(),
          fs       ::Float64                   = zero(Float64),
          gain     ::Float64                   = one(Float64),
          freqmin  ::Float64                   = zero(Float64),
          freqmax  ::Float64                   = zero(Float64),
          cc_len   ::Int                       = zero(Int),
          cc_step  ::Int                       = zero(Int),
          whitened ::Bool                      = false,
          time_norm::Union{Bool,String}        = false,
          resp     ::PZResp                    = PZResp(),
          misc     ::Dict{String,Any}          = Dict{String,Any}(),
          notes    ::Array{String,1}           = Array{String,1}(undef, 0),
          t        ::Array{Float64,1}          = Array{Float64,1}(undef, 0),
          fft      ::Array{<:Union{Complex{Float32},Complex{Float64}},2} =
                     Array{Complex{Float32},2}(undef, 0, 2)
          ) = FFTData(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
                     whitened, time_norm, resp, misc, notes, t, fft)

FFTData(C::SeisChannel,freqmin::Float64, freqmax::Float64,cc_len::Int,
        cc_step::Int, whitened::Bool,time_norm::Union{Bool,String},
        t::Array{Float64,1},fft::Array{<:Union{Complex{Float32},Complex{Float64}},2}
       ) = FFTData(C.id, Dates.format(u2d(C.t[1,2]*1e-6),"Y-mm-dd"), C.loc, C.fs, C.gain, freqmin, freqmax, cc_len,
                   cc_step, whitened, time_norm, C.resp, C.misc, C.notes, t, fft)

in(s::String, F::FFTData) = F.id==s

isempty(F::FFTData) = min(F.name=="",F.id == "",isempty(F.loc),
                          F.fs == zero(typeof(F.fs)),
                          F.gain == one(typeof(F.gain)),
                          F.freqmin == zero(typeof(F.freqmin)),
                          F.freqmax == zero(typeof(F.freqmax)),
                          F.cc_len == zero(typeof(F.cc_len)),
                          F.cc_step == zero(typeof(F.cc_step)),
                          F.whitened == false,
                          F.time_norm == false,
                          isempty(F.resp), isempty(F.misc),isempty(F.notes),
                          isempty(F.t),isempty(F.fft))

isequal(F::FFTData, U::FFTData) = minimum([hash(getfield(F,i))==hash(getfield(U,i)) for i in fftfields]::Array{Bool,1})
==(F::FFTData, U::FFTData) = isequal(F,U)::Bool

append!(F::FFTData, U::FFTData)  = (if F == U;
  setfield!(F, :t, vcat(getfield(F,:t), getfield(U,:t)));
  setfield!(F, :fft, hcat(getfield(F,:fft), getfield(U,:fft)))
  elseif isempty(F);
  [setfield!(F, i, getfield(U,i)) for i in fftfields];
  [setfield!(F, i, getfield(U,i)) for i in [:t,:fft]];
  end;
  return F)
+(F::FFTData, U::FFTData) = (T = deepcopy(F); return append!(T, U))

function sizeof(F::FFTData)
s = sum([sizeof(getfield(F,f)) for f in fftfields])
if !isempty(F.notes)
  s += sum([sizeof(i) for i in F.notes])
end
if !isempty(F.misc)
  s += sum([sizeof(i) for i in values(F.misc)])
end
return s
end


text:
using MetacommunityDynamics
using DynamicGrids
using Distributions
using ColorSchemes
using Dispersal: DispersalKernel, OutwardsDispersal
using Plots

consumermodel =
    Eating{Tuple{:C,:R}}(
        functionalresponse=LotkaVolterra(0.2),
        dt=0.1) +
        OutwardsDispersal{:C}(DispersalKernel(radius=3)) +
    #RandomExtinction{:C}(0.1) +
    LinearMortality{:C}(0.1);

resourcemodel =
    LogisticGrowth{:R}( λ=1.1, K=200., dt=0.1) +
    OutwardsDispersal{:R}(DispersalKernel(radius=5));

model = resourcemodel + consumermodel

# helper for plotting
function plot_lv_timeseries(arrayout)
    Cs = Float64[]
    Rs = Float64[]
    for t in arrayout
        push!(Cs, sum(t[:C]))
        push!(Rs, sum(t[:R]))
    end

    plot(1:length(Cs),Cs, label="consumer")
    plot!(1:length(Rs),Rs, label="resource")
    xlabel!("time")
    ylabel!("biomass")
end

# wave initial condition
latticesize = 250

initresource = rand(Biomass, Exponential(10.) ,latticesize, latticesize)

initconsumer = zeros(Biomass, latticesize, latticesize)
rci = rand(DiscreteUniform(1, latticesize), (100, 2) )
map( i -> initconsumer[i...] = 10., rci)

arrayout = ArrayOutput((C=initconsumer, R=initresource ), tspan=1:500)
@time sim!(arrayout, model)

max(arrayout[300][:R]...)

plot_lv_timeseries(arrayout)

gifoutput = GifOutput((C=initconsumer, R=initresource ),
    filename="./docs/static/spicylv.gif",
    tspan=1:750,
    fps=25,
    minval=(0,0), maxval=(10, 30),
    scheme = (ColorSchemes.Blues_5, ColorSchemes.Greens_4),
    padval = 100
)
@time sim!(gifoutput, model)


text:
using HomotopyOpt, Test

@testset "HomotopyOpt" begin
    include("2DTests.jl")
    include("3DTests.jl")
end


text:
using StatsModels

@test terms(rcopy(R"y ~ x + (1 | g)")) == terms(@formula y ~ x + (1 | g))
@test terms(rcopy(R"y ~ x + (1 | g) + (1 | d)")) == terms(@formula y ~ x + (1 | g) + (1 | d))
@test terms(rcopy(R"y ~ (a + b) * c")) == terms(@formula y ~ (a + b) * c)
@test terms(rcopy(R"y ~ c : (a + b)")) == terms(@formula y ~ c & (a + b))

# testing association
@test terms(rcopy(R"y ~ a + b + c")) == terms(@formula y ~ a + b + c)
@test terms(rcopy(R"y ~ a + b * c * d")) == terms(@formula y ~ a + b * c * d)
@test terms(rcopy(R"y ~ a : b : c")) == terms(@formula y ~ a & b & c)
@test terms(rcopy(R"y ~ a + b : c : d")) == terms(@formula y ~ a + b & c & d)

@test rcopy(rcall(Symbol("=="), R"y ~ x + (1 | g)", RObject(@formula y ~ x + (1 | g))))
@test rcopy(rcall(Symbol("=="), R"y ~ x + (1 | g) + (1 | d)", RObject(@formula y ~ x + (1 | g) + (1 | d))))
@test rcopy(rcall(Symbol("=="), R"y ~ a + b + c + a:c + b:c", RObject(@formula y ~ (a + b) * c)))
@test rcopy(rcall(Symbol("=="), R"y ~ c + a + b + c:a + c:b", RObject(@formula y ~ c * (a + b))))

# testing association
@test rcopy(rcall(Symbol("=="), R"y ~ a + b + c", RObject(@formula y ~ a + b + c)))
@test rcopy(rcall(Symbol("=="), R"y ~ a + b + c + d + b:c + b:d + c:d + b:c:d", RObject(@formula y ~ a + b * c * d)))
@test rcopy(rcall(Symbol("=="), R"y ~ a : b : c", RObject(@formula y ~ a & b & c)))
@test rcopy(rcall(Symbol("=="), R"y ~ d + a : b : c", RObject(@formula y ~ d + a & b & c)))


text:
macro ce_curl(f, handle, args...)
  local esc_args = [esc(arg) for arg in args]
  quote
    cc = $(esc(f))($(esc(handle)), $(esc_args...))

    if cc != CURLE_OK
      err = unsafe_string(curl_easy_strerror(cc))
      error(string($f) * "() failed: " * err)
    end
  end
end

macro ce_curlm(f, handle, args...)
  local esc_args = [esc(arg) for arg in args]
  quote
    cc = $(esc(f))($(esc(handle)), $(esc_args...))

    if cc != CURLM_OK
      err = unsafe_string(curl_multi_strerror(cc))
      error(string($f) * "() failed: " * err)
    end
  end
end


text:
using Compat: cispi
# using DoubleFloats
using FourierTransforms
using Random
using Test

const FT = FourierTransforms

# `DoubleFloat` does not work: <https://github.com/JuliaMath/DoubleFloats.jl/issues/139>
const Types = [Complex{Float32}, Complex{Float64}, Complex{BigFloat}]

# # Random.rand(rng::AbstractRNG, ::Random.Random.Random.CloseOpen01{DoubleFloat}) = ???
# Random.randn(::Type{DoubleFloat}, dims::Integer...) = DoubleFloat.(randn(Float64, dims...))
# Random.randn(::Type{Complex{DoubleFloat}}, dims::Integer...) = Complex{DoubleFloat}.(randn(Complex{Float64}, dims...))
# Base.eps(::Type{DoubleFloat}) = eps(Float64)^2

# Random.rand(rng::AbstractRNG, ::Random.Random.Random.CloseOpen01{BigFloat}) = ???
Random.randn(::Type{BigFloat}, dims::Integer...) = BigFloat.(randn(Float64, dims...))
Random.randn(::Type{Complex{BigFloat}}, dims::Integer...) = Complex{BigFloat}.(randn(Complex{Float64}, dims...))

Random.seed!(100)
@testset "Direct Fourier transforms (T=$T)" for T in Types
    maxbits = T === Complex{BigFloat} ? 8 : 10
    for bits in 1:maxbits
        N = 2^bits
        x = randn(T, N)
        X = FT.direct_ft(x)
        x′ = conj(FT.direct_ft(conj(X)) / N)
        @test x′ ≈ x
    end
end

Random.seed!(100)
@testset "Radix-2 Cooley–Tukey Fast Fourier transforms (T=$T)" for T in Types
    maxbits = T === Complex{BigFloat} ? 15 : 20
    for bits in 1:maxbits
        N = 2^bits
        x = randn(T, N)
        X = FT.ditfft2(x)
        x′ = conj(FT.ditfft2(conj(X)) / N)
        @test x′ ≈ x
    end
end

Random.seed!(100)
@testset "General radix Cooley–Tukey Fast Fourier transforms (T=$T)" for T in Types
    maxbits = T === Complex{BigFloat} ? 15 : 20
    for bits in 1:maxbits
        N = 2^bits
        x = randn(T, N)
        X = FT.radix_fft(x)
        x′ = conj(FT.radix_fft(conj(X)) / N)
        @test x′ ≈ x
    end
end

Random.seed!(100)
@testset "Comparing algorithms for Fourier transforms (T=$T)" for T in Types
    maxbits = T === Complex{BigFloat} ? 8 : 10
    for bits in 1:maxbits
        N = 2^bits
        x = randn(T, N)
        X = FT.direct_ft(x)
        if N == 1
            X1 = FT.direct_ft_1(x)
            @test X1 ≈ X
        end
        if N == 2
            X2 = FT.direct_ft_2(x)
            @test X2 ≈ X
        end
        if N == 4
            X4 = FT.direct_ft_4(x)
            @test X4 ≈ X
        end
        Xct2 = FT.ditfft2(x)
        @test Xct2 ≈ X
        Xr = FT.radix_fft(x)
        @test Xr ≈ X
        for radix in 2:8
            if N % radix == 0
                Xr = FT.radix_fft(x; choose_radix=N -> radix)
                @test Xr ≈ X
            end
        end
    end
end

Random.seed!(100)
@testset "Simple Fourier transforms (T=$T)" for T in Types
    RT = typeof(real(zero(T)))
    maxbits = T === Complex{BigFloat} ? 15 : 20
    for bits in 1:maxbits
        for iter in 1:10
            N = 2^bits
            k = rand(1:N)
            x = [cispi(2 * ((k - 1) * (n - 1) % N) / RT(N)) for n in 1:N]
            X₀ = [(k′ == k) * N for k′ in 1:N]
            X = fft(x)
            @test isapprox(X, X₀; atol=N^2 * sqrt(eps()))
        end
    end
end

Random.seed!(100)
@testset "Inverse Fourier transforms (T=$T)" for T in Types
    maxbits = T === Complex{BigFloat} ? 15 : 20
    for bits in 1:maxbits
        for iter in 1:10
            N = 2^bits
            x = randn(T, N)
            X = fft(x)
            x′ = inv_fft(X)
            @test x′ ≈ x
        end
    end
end

Random.seed!(100)
@testset "Linearity of Fourier transforms (T=$T)" for T in Types
    maxbits = T === Complex{BigFloat} ? 10 : 15
    for bits in 1:maxbits
        for iter in 1:10
            N = 2^bits
            x = randn(T, N)
            y = randn(T, N)
            α = randn(T)
            z = α * x + y
            X = fft(x)
            Y = fft(y)
            Z = fft(z)
            @test α * X + Y ≈ Z
        end
    end
end


text:
module JuliaForBioinformatics

using InteractiveUtils

export show_type_tree,
	   cp_notebooks

include("type_tree.jl")

function cp_notebooks(dest=homedir())
	note_path = abspath(@__DIR__, "..", "docs", "src", "notebooks")
	for file in readdir(note_path)
		@show file
		try
			cp(joinpath(note_path, file), joinpath(dest, file), force=false)
			println("       copied!")
		catch
			println("       skipped!")
		end
	end
end

end # module

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
module BoxCoxTrans

using Optim: optimize, minimizer
using Statistics: mean, var
using StatsBase: geomean

"""
    transform(𝐱)

Transform an array using Box-Cox method.  The power parameter λ is derived
from maximizing a log-likelihood estimator.

If the array contains any non-positive values then a DomainError is thrown.
This can be avoided by providing the shift parameter α to make all values
positive.

Keyword arguments:
- α: added to all values in 𝐱 before transformation. Default = 0.
- scaled: scale transformation results.  Default = false.
"""
function transform(𝐱; α = 0, kwargs...)
    if α != 0
        𝐱 .+= α
    end
    λ, details = lambda(𝐱; kwargs...)
    #@info "estimated lambda = $λ"
    transform(𝐱, λ; kwargs...)
end

"""
    transform(𝐱, λ; α = 0)

Transform an array using Box-Cox method with the provided power parameter λ.
If the array contains any non-positive values then a DomainError is thrown.

Keyword arguments:
- α: added to all values in 𝐱 before transformation. Default = 0.
- scaled: scale transformation results.  Default = false.
"""
function transform(𝐱, λ; α = 0, scaled = false, kwargs...)
    if α != 0
        𝐱 .+= α
    end
    any(𝐱 .<= 0) && throw(DomainError("Data must be positive and ideally greater than 1.  You may specify α argument(shift). "))
    if scaled
        gm = geomean(𝐱)
        @. λ ≈ 0 ? gm * log(𝐱) : (𝐱 ^ λ - 1) / (λ * gm ^ (λ - 1))
    else
        @. λ ≈ 0 ? log(𝐱) : (𝐱 ^ λ - 1) / λ
    end
end

"""
    lambda(𝐱; interval = (-2.0, 2.0), method = :geomean)

Calculate lambda from an array using a log-likelihood estimator.

Keyword arguments:
- method: either :geomean or :normal
- any other keyword arguments accepted by Optim.optimize function e.g. abs_tol

See also: [`log_likelihood`](@ref)
"""
function lambda(𝐱 ; α = 0 , interval = (-2.0, 2.0), kwargs...)
    if α != 0
        𝐱 .+= α
    end

    i1, i2 = interval
    res = optimize(λ -> -log_likelihood(𝐱, λ; kwargs...), i1, i2)
    (value=minimizer(res), details=res)
end

"""
    log_likelihood(𝐱, λ; method = :geomean)

Return log-likelihood for the given array and lambda.

Method :geomean =>
    -N / 2.0 * log(2 * π * σ² / gm ^ (2 * (λ - 1)) + 1)

Method :normal =>
    -N / 2.0 * log(σ²) + (λ - 1) * sum(log.(𝐱))
"""
function log_likelihood(𝐱, λ; method = :geomean, kwargs...)
    N = length(𝐱)
    𝐲 = transform(float.(𝐱), λ)
    σ² = var(𝐲, corrected = false)
    gm = geomean(𝐱)
    if method == :geomean
        -N / 2.0 * log(2 * π * σ² / gm ^ (2 * (λ - 1)) + 1)
    elseif method == :normal
        -N / 2.0 * log(σ²) + (λ - 1) * sum(log.(𝐱))
    else
        throw(ArgumentError("Incorrect method. Please specify :geomean or :normal."))
    end
end

"""
    retransform(𝐱, λ; α = 0)

Retransform an array which is transformed using Box-Cox method with the provided power parameter λ and shift
argument α to the oreginal array.

Keyword arguments:
- α: added to all values in 𝐱 before transformation. Default = 0.
- scaled: scale transformation results.  Default = false.
"""

function backtransform(𝐱, λ; α = 0, scaled = false, kwargs...)
    if scaled
        gm = geomean(𝐱)
        @. λ ≈ 0 ? exp.(𝐱 / gm) - α  : (𝐱 * λ * gm ^ (λ -1) +1) ^ (1 / λ) - α
    else
        @. λ ≈ 0 ? exp.(𝐱) - α : (λ * 𝐱 + 1) ^ (1 / λ) - α
    end
end

end # module


text:
using Base: contractuser
using Base.Filesystem: path_separator

using LibGit2: LibGit2, GitCommit, GitRemote, GitRepo
using Pkg: Pkg, PackageSpec, TOML
using Random: Random, randstring
using Test: @test, @testset, @test_broken, @test_logs, @test_throws

using DeepDiffs: deepdiff
using SimpleMock: mock
using Suppressor: @suppress

using PkgTemplates
const PT = PkgTemplates

const USER = "tester"

Random.seed!(1)

# Creata a template that won't error because of a missing username.
tpl(; kwargs...) = Template(; user=USER, kwargs...)

# Generate a random package name.
pkgname() = titlecase(randstring('A':'Z', 16))

# Create a randomly named package with a template, and delete it afterwards.
function with_pkg(f::Function, t::Template, pkg::AbstractString=pkgname())
    @suppress t(pkg)
    try
        f(pkg)
    finally
        # On 1.4, this sometimes won't work, but the error is that the package isn't installed.
        # We're going to delete the package directory anyways, so just ignore any errors.
        PT.version_of(pkg) === nothing || try @suppress Pkg.rm(pkg) catch; end
        rm(joinpath(t.dir, pkg); recursive=true, force=true)
    end
end

function print_diff(a, b)
    old = Base.have_color
    @eval Base have_color = true
    try
        println(deepdiff(a, b))
    finally
        @eval Base have_color = $old
    end
end

mktempdir() do dir
    Pkg.activate(dir)
    pushfirst!(DEPOT_PATH, dir)
    try
        @testset "PkgTemplates.jl" begin
            include("template.jl")
            include("plugin.jl")
            include("show.jl")
            include("interactive.jl")

            if PT.git_is_installed()
                include("git.jl")

                # Quite a bit of output depends on the Julia version,
                # and the test fixtures are made with Julia 1.4.
                # TODO: Keep this on the latest stable Julia version.
                if VERSION.major == 1 && VERSION.minor == 4
                    include("reference.jl")
                else
                    @info "Skipping reference tests" VERSION
                end
            else
                @info "Git is not installed, skipping Git and reference tests"
            end
        end
    finally
        popfirst!(DEPOT_PATH)
    end
end


text:
module Grids

# Each grid lives on a process

using Foldable, Functor, StencilFunctor
using Cells, Defs, Norms
using Comm

import Base.show
import Norms.norm
import Cells.axpy, Cells.initial, Cells.error, Cells.rhs

export Grid
export getBoundary, axpy, norm, initial, error, rhs

immutable Grid
    imin::Int                   # spatial indices
    imax::Int
    cells::Vector{Cell}
    function Grid(imin, imax, cells)
        # println("Running on process $(Comm.comminfo.rank)")
        @assert imax>imin
        new(imin, imax, cells)
    end
end

function show(io::IO, g::Grid)
    println(io, "    grid: imin=", g.imin, " imax=", g.imax)
    for (i,c) in enumerate(g.cells)
        println(io, "        i=", i+g.imin-1, " ", c)
    end
end

imax(imin::Int) = min(imin + Defs.ncells_per_grid - 1, Defs.ncells())
x(i::Int) = Defs.xmin + (float(i) - 0.5) * Defs.dx()

getBoundary(g::Grid, face::Bool) = face ? g.cells[end] : g.cells[1]

# Linear combination
axpy(a::Float64, x::Grid, y::Grid) =
    Grid(y.imin, y.imax,
         fmap(R=Cell, (x,y) -> Cells.axpy(a,x,y), x.cells, y.cells))

# Norm
norm(g::Grid) = freduce((n,c) -> n+norm(c), Norm(), g.cells)

# Initial condition
initial(t::Float64, imin::Int) =
    Grid(imin, imax(imin), [Cells.initial(t, x(i)) for i in imin:imax(imin)])

# Error
error(g::Grid, t::Float64) =
    Grid(g.imin, g.imax, fmap(R=Cell, c -> Cells.error(c,t), g.cells))

# RHS
rhs(g::Grid, bm::Cell, bp::Cell) =
    Grid(g.imin, g.imax,
         stencil_fmap(R=Cell, Cells.rhs, (c,f)->c, g.cells, bm, bp))

end


text:
export AbstractAgent,
    get_role,
    PreExperimentStage,
    PostExperimentStage,
    PreEpisodeStage,
    PostEpisodeStage,
    PreActStage,
    PostActStage,
    PRE_EXPERIMENT_STAGE,
    POST_EXPERIMENT_STAGE,
    PRE_EPISODE_STAGE,
    POST_EPISODE_STAGE,
    PRE_ACT_STAGE,
    POST_ACT_STAGE,
    Training,
    Testing

"""
    (agent::AbstractAgent)(env) = agent(PRE_ACT_STAGE, env) -> action
    (agent::AbstractAgent)(stage::AbstractStage, env)

Similar to [`AbstractPolicy`](@ref), an agent is also a functional object which takes in an observation and returns an action.
The main difference is that, we divide an experiment into the following stages:

- `PRE_EXPERIMENT_STAGE`
- `PRE_EPISODE_STAGE`
- `PRE_ACT_STAGE`
- `POST_ACT_STAGE`
- `POST_EPISODE_STAGE`
- `POST_EXPERIMENT_STAGE`

In each stage, different types of agents may have different behaviors, like updating experience buffer, environment model or policy.
"""
abstract type AbstractAgent end

function get_role(::AbstractAgent) end

"""
                      +-----------------------------------------------------------+
                      |Episode                                                    |
                      |                                                           |
PRE_EXPERIMENT_STAGE  |            PRE_ACT_STAGE    POST_ACT_STAGE                | POST_EXPERIMENT_STAGE
         |            |                  |                |                       |          |
         v            |        +-----+   v   +-------+    v   +-----+             |          v
         --------------------->+ env +------>+ agent +------->+ env +---> ... ------->......
                      |  ^     +-----+       +-------+ action +-----+          ^  |
                      |  |                                                     |  |
                      |  +--PRE_EPISODE_STAGE            POST_EPISODE_STAGE----+  |
                      |                                                           |
                      |                                                           |
                      +-----------------------------------------------------------+
"""
abstract type AbstractStage end

struct PreExperimentStage <: AbstractStage end
struct PostExperimentStage <: AbstractStage end
struct PreEpisodeStage <: AbstractStage end
struct PostEpisodeStage <: AbstractStage end
struct PreActStage <: AbstractStage end
struct PostActStage <: AbstractStage end

const PRE_EXPERIMENT_STAGE = PreExperimentStage()
const POST_EXPERIMENT_STAGE = PostExperimentStage()
const PRE_EPISODE_STAGE = PreEpisodeStage()
const POST_EPISODE_STAGE = PostEpisodeStage()
const PRE_ACT_STAGE = PreActStage()
const POST_ACT_STAGE = PostActStage()

(agent::AbstractAgent)(env) = agent(PRE_ACT_STAGE, env)
function (agent::AbstractAgent)(stage::AbstractStage, env) end

struct Training{T<:AbstractStage} end
Training(s::T) where {T<:AbstractStage} = Training{T}()
struct Testing{T<:AbstractStage} end
Testing(s::T) where {T<:AbstractStage} = Testing{T}()

Base.show(io::IO, agent::AbstractAgent) =
    AbstractTrees.print_tree(io, StructTree(agent), get(io, :max_depth, 10))


text:
export BlockSparseMatrixCSC
@computed struct BlockSparseMatrixCSC{B,T,I} <:
                 SparseArrays.AbstractSparseMatrixCSC{T,I}
    elts::SparseMatrixCSC{SMatrix{B,B,T,B^2},I}

    function BlockSparseMatrixCSC{B,T,I}(bA::SparseMatrixCSC{SMatrix{B,B,T,L},
                                                             I}) where {B,T,
                                                                        I<:Integer,
                                                                        L}
        B::Integer
        @assert L == B^2
        @assert T ≢ Any && T ≢ Union{} # TODO
        return new{B,T,I}(bA)
    end
    function BlockSparseMatrixCSC{B,T}(bA::SparseMatrixCSC{SMatrix{B,B,T,L},I}) where {B,
                                                                                T,
                                                                                I<:Integer,
                                                                                L}
        return BlockSparseMatrixCSC{B,T,I}(bA)
    end
    function BlockSparseMatrixCSC{B}(bA::SparseMatrixCSC{SMatrix{B,B,T,L},I}) where {B,
                                                                                T,
                                                                                I<:Integer,
                                                                                L}
        return BlockSparseMatrixCSC{B,T,I}(bA)
    end
    function BlockSparseMatrixCSC(bA::SparseMatrixCSC{SMatrix{B,B,T,L},I}) where {B,
                                                                                T,
                                                                                I<:Integer,
                                                                                L}
        return BlockSparseMatrixCSC{B,T,I}(bA)
    end
end

Base.show(io::IOContext, A::BlockSparseMatrixCSC) = show(io, SparseMatrixCSC(A))
# Base.show(io::IOContext, A::BlockSparseMatrixCSC) = show(io, A.elts)

###

export blocksparse
function blocksparse(::Val{B}, Is::AbstractVector{I}, Js::AbstractVector{I},
                     Vs::AbstractVector{T}, ni::Integer,
                     nj::Integer) where {B,T,I<:Integer}
    B::Integer
    length(Is) == length(Js) == length(Vs) ||
        throw(ArgumentError("index and value vectors must have the same length"))
    ni ≥ 0 && nj ≥ 0 ||
        throw(ArgumentError("Number of elements must be non-negative"))
    @assert ni % B == 0
    @assert nj % B == 0
    ni0 = ni ÷ B
    nj0 = nj ÷ B
    Ibs = I[]
    Jbs = I[]
    Vbs = SMatrix{B,B,T,B^2}[]
    curi0 = curj0 = 0
    v0 = zero(SMatrix{B,B,T})
    for (i, j, v) in zip(Is, Js, Vs)
        @assert 1 ≤ i ≤ ni
        @assert 1 ≤ j ≤ nj
        i0, i1 = fldmod1(i, B)
        j0, j1 = fldmod1(j, B)
        if (curi0, curj0) ≠ (i0, j0)
            if curi0 ≠ 0
                push!(Ibs, curi0)
                push!(Jbs, curj0)
                push!(Vbs, v0)
            end
            curi0, curj0 = i0, j0
            v0 = zero(SMatrix{B,B,T})
        end
        v0 = setindex(v0, v, i1, j1)
    end
    if curi0 ≠ 0
        push!(Ibs, curi0)
        push!(Jbs, curj0)
        push!(Vbs, v0)
    end
    return BlockSparseMatrixCSC(sparse(Ibs, Jbs, Vbs, ni0, nj0))
end

function blocksparse(::Val{B}, A::AbstractMatrix{T}) where {B,T}
    ni, nj = size(A)
    @assert ni % B == 0
    @assert nj % B == 0
    ni0 = ni ÷ B
    nj0 = nj ÷ B
    C = Array{SMatrix{B,B,T,B^2}}(undef, ni0, nj0)
    for j0 in 1:nj0, i0 in 1:ni0
        C1 = zero(SMatrix{B,B,T})
        for j1 in 1:B, i1 in 1:B
            C1 = setindex(C1, A[B * (i0 - 1) + i1, B * (j0 - 1) + j1], i1, j1)
        end
        C[i0, j0] = C1
    end
    return BlockSparseMatrixCSC(sparse(C))
end

export blocksize
blocksize(x::BlockSparseMatrixCSC{B}) where {B} = B

for (fun, spfun) in [(:bsprand, :sprand), (:bsprandn, :sprandn)]
    @eval begin
        export $fun
        function $fun(r::AbstractRNG, ::Val{B}, ::Type{T}, ni::Integer,
                      nj::Integer, p::AbstractFloat) where {B,T}
            @assert ni % B == 0
            @assert nj % B == 0
            ni0 = ni ÷ B
            nj0 = nj ÷ B
            return BlockSparseMatrixCSC($spfun(r, SMatrix{B,B,T,B^2}, ni0, nj0,
                                               p))
        end
        function $fun(r::AbstractRNG, ::Val{B}, ni::Integer, nj::Integer,
                      p::AbstractFloat) where {B}
            return $fun(r, Val(B), Float64, ni, nj, p)
        end
        function $fun(::Val{B}, ::Type{T}, ni::Integer, nj::Integer,
                      p::AbstractFloat) where {B,T}
            return $fun(Random.default_rng(), Val(B), T, ni, nj, p)
        end
        function $fun(::Val{B}, ni::Integer, nj::Integer,
                      p::AbstractFloat) where {B}
            return $fun(Random.default_rng(), Val(B), ni, nj, p)
        end
    end
end

export bspzeros
function bspzeros(::Val{B}, ::Type{T}, ::Type{I}, ni::Integer,
                  nj::Integer) where {B,T,I<:Integer}
    B::Integer
    @assert 0 ≤ ni && ni % B == 0
    @assert 0 ≤ nj && nj % B == 0
    return BlockSparseMatrixCSC(spzeros(SMatrix{B,B,T,B^2}, fld1(ni, B),
                                        fld1(nj, B)))
end
function bspzeros(::Val{B}, ::Type{T}, ni::Integer, nj::Integer) where {B,T}
    return bspzeros(Val(B), T, Int, ni, nj)
end
function bspzeros(::Val{B}, ni::Integer, nj::Integer) where {B}
    return bspzeros(Val(B), Float64, ni, nj)
end

function (Matrix{T})(A::BlockSparseMatrixCSC{B}) where {T,B}
    R = zeros(T, size(A))
    for j0 in 1:size(A.elts, 2), k in nzrange(A.elts, j0)
        i0 = rowvals(A.elts)[k]
        v0 = nonzeros(A.elts)[k]
        for j1 in 1:B, i1 in 1:B
            R[B * (i0 - 1) + i1, B * (j0 - 1) + j1] = v0[i1, j1]
        end
    end
    return R
end
(Array{T})(A::BlockSparseMatrixCSC) where {T} = Matrix{T}(A)
Matrix(A::BlockSparseMatrixCSC{B,T}) where {B,T} = Matrix{T}(A)
Array(A::BlockSparseMatrixCSC) = Matrix(A)

### Comparison

Base.:(==)(A::BlockSparseMatrixCSC, B::BlockSparseMatrixCSC) = A.elts == B.elts

### Collection

function Base.copy!(dst::BlockSparseMatrixCSC{B},
                    src::BlockSparseMatrixCSC{B}) where {B}
    return BlockSparseMatrixCSC(copy!(dst.elts, src.elts))
end
Base.copy(A::BlockSparseMatrixCSC) = BlockSparseMatrixCSC(copy(A.elts))

Base.eltype(A::BlockSparseMatrixCSC{B,T}) where {B,T} = T
Base.eltype(::Type{<:BlockSparseMatrixCSC{B,T}}) where {B,T} = T

function Base.getindex(A::BlockSparseMatrixCSC{B}, i::Integer,
                       j::Integer) where {B}
    i0, i1 = fldmod1(i, B)
    j0, j1 = fldmod1(j, B)
    return A.elts[i0, j0][i1, j1]
end
function Base.getindex(A::BlockSparseMatrixCSC, i::CartesianIndex{2})
    return getindex(A, i[1], i[2])
end

function Base.similar(A::BlockSparseMatrixCSC; element_type)
    return BlockSparseMatrixCSC(similar(A.elts; element_type))
end

Base.size(A::BlockSparseMatrixCSC) = blocksize(A) .* size(A.elts)

### Sparse

function blocksparse(::Val{B}, A::SparseMatrixCSC{T,I}) where {B,T,I}
    # I, J, V = findnz(A)
    # TODO: Make this more efficient
    colvals = Array{I}(undef, nnz(A))
    for j in 1:size(A, 2)
        colvals[nzrange(A, j)] .= j
    end
    return blocksparse(Val(B), rowvals(A), colvals, nonzeros(A), size(A)...)
end

function (SparseMatrixCSC{T,I})(A::BlockSparseMatrixCSC{B}) where {T,B,I}
    Is = I[]
    Js = I[]
    Vs = T[]
    for j0 in 1:size(A.elts, 2), k0 in nzrange(A.elts, j0)
        i0 = rowvals(A.elts)[k0]
        v0 = nonzeros(A.elts)[k0]
        for j1 in 1:B, i1 in 1:B
            if !iszero(v0[i1, j1])
                push!(Is, B * (i0 - 1) + i1)
                push!(Js, B * (j0 - 1) + j1)
                push!(Vs, v0[i1, j1])
            end
        end
    end
    return sparse(Is, Js, Vs, size(A)...)
end
function (SparseMatrixCSC{T})(A::BlockSparseMatrixCSC{B,U,I}) where {T,B,U,I}
    return SparseMatrixCSC{T,I}(A)
end
function SparseMatrixCSC(A::BlockSparseMatrixCSC{B,T}) where {B,T}
    return SparseMatrixCSC{T}(A)
end

#TODO SparseArrays.droptol!(x::BlockSparseVector, tol) = fkeep!(x, (i, a) -> abs(a) > tol)
#TODO
#TODO SparseArrays.dropzeros!(x::BlockSparseVector) = fkeep!(x, (i, a) -> !iszero(a))
#TODO SparseArrays.dropzeros(x::BlockSparseVector) = dropzeros!(copy(x))
#TODO
#TODO function SparseArrays.fkeep!(x::BlockSparseVector{B}, f) where {B}
#TODO     fkeep!(x.elts, (i0, y) -> any(f(b * (i0 - 1) + i1, y[i1]) for i1 in 1:B))
#TODO     return x
#TODO end

SparseArrays.indtype(::BlockSparseMatrixCSC{B,T}) where {B,T,I} = I
SparseArrays.indtype(::Type{<:BlockSparseMatrixCSC{B,T,I}}) where {B,T,I} = I

SparseArrays.nnz(A::BlockSparseMatrixCSC) = blocksize(A)^2 * nnz(A.elts)

### Vector space

fixeltype(::T, A::SparseMatrixCSC{T}) where {T} = A
fixeltype(::T, A::SparseMatrixCSC{Any}) where {T} = SparseMatrixCSC{T}(A)

function Base.zero(A::BlockSparseMatrixCSC{B,T,I}) where {B,T,I}
    return BlockSparseMatrixCSC{B,T,I}(spzeros(SMatrix{B,B,T}, size(A.elts)))
end
Base.:+(A::BlockSparseMatrixCSC) = A # BlockSparseMatrixCSC(+A.elts)
Base.:-(A::BlockSparseMatrixCSC) = BlockSparseMatrixCSC(-A.elts)
function Base.:+(A::BlockSparseMatrixCSC{B},
                 C::BlockSparseMatrixCSC{B}) where {B}
    return BlockSparseMatrixCSC(A.elts + C.elts)
end
function Base.:-(A::BlockSparseMatrixCSC{B},
                 C::BlockSparseMatrixCSC{B}) where {B}
    return BlockSparseMatrixCSC(A.elts - C.elts)
end
function Base.:*(a::Number, A::BlockSparseMatrixCSC)
    return BlockSparseMatrixCSC(fixeltype(a * one(eltype(A.elts)), a * A.elts))
end
function Base.:\(a::Number, A::BlockSparseMatrixCSC)
    return BlockSparseMatrixCSC(fixeltype(a \ one(eltype(A.elts)), a \ A.elts))
end
function Base.:*(A::BlockSparseMatrixCSC, a::Number)
    return BlockSparseMatrixCSC(fixeltype(one(eltype(A.elts)) * a, A.elts * a))
end
function Base.:/(A::BlockSparseMatrixCSC, a::Number)
    return BlockSparseMatrixCSC(fixeltype(one(eltype(A.elts)) / a, A.elts / a))
end

### Category
function Base.one(A::BlockSparseMatrixCSC{B,T,I}) where {B,T,I}
    @assert size(A, 1) == size(A, 2)
    return BlockSparseMatrixCSC{B,T,I}(sparse(Diagonal([one(SMatrix{B,B,T})
                                                        for i in
                                                            1:size(A.elts, 1)])))
end
function Base.:*(A::BlockSparseMatrixCSC{B},
                 C::BlockSparseMatrixCSC{B}) where {B}
    @assert size(A, 2) == size(C, 1)
    return BlockSparseMatrixCSC(fixeltype(one(eltype(A.elts)) *
                                          one(eltype(C.elts)), A.elts * C.elts))
end

### Groupoid
function Base.:/(A::BlockSparseMatrixCSC{B},
                 C::BlockSparseMatrixCSC{B}) where {B}
    # A ⋅ inv(C)
    @assert size(A, 2) == size(C, 2)
    return BlockSparseMatrixCSC(A.elts / C.elts)
end
function Base.:\(A::BlockSparseMatrixCSC{B},
                 C::BlockSparseMatrixCSC{B}) where {B}
    # inv(A) ⋅ C
    @assert size(A, 1) == size(C, 1)
    return BlockSparseMatrixCSC(A.elts \ C.elts)
end

### Operator
function Base.:*(A::BlockSparseMatrixCSC{B}, x::BlockSparseVector{B}) where {B}
    @assert size(A, 2) == size(x, 1)
    return BlockSparseVector(A.elts * x.elts)
end
function Base.:\(A::BlockSparseMatrixCSC{B}, x::BlockSparseVector{B}) where {B}
    @assert size(A, 1) == size(x, 1)
    return BlockSparseVector(A.elts \ x.elts)
end

### Linear algebra
function LinearAlgebra.norm(A::BlockSparseMatrixCSC{B,T}, p::Real=2) where {B,T}
    if p == 0
        return zero(norm(zero(T), 0)^0 + norm(one(T), 0)^0) + length(A)
    elseif p == 1
        r = zero(norm(zero(T), 1) + norm(one(T), 1))
        for v0 in nonzeros(A.elts)
            for v in v0
                r += norm(v, 1)
            end
        end
        return r
    elseif p == 2
        r = zero(norm(zero(T), 2)^2 + norm(one(T), 2)^2)
        for v0 in nonzeros(A.elts)
            for v in v0
                r += norm(v, 2)^2
            end
        end
        return sqrt(r)
    elseif p == Inf
        r = zero(max(norm(zero(T), Inf), norm(one(T), Inf)))
        for v0 in nonzeros(A.elts)
            for v in v0
                r = max(r, norm(v, Inf))
            end
        end
        return r
    else
        r = zero(zero(T) + norm(one(T), p)^p)
        for v0 in nonzeros(A.elts)
            for v in v0
                r += norm(v, p)^p
            end
        end
        return r^inv(p)
    end
end


text:
println("Computing the amount of wrapping paper the elves need...")
# Open input
cd(raw"/Users/wapisani/Documents/Programming/Julia/Advent_of_Code/")
open("AoC_2015_2.input","r") do f
    global wrapping_dimensions_list = split(read(f,String),"\n")
end

surface_area = 0
length_ribbon = 0
for dimensions in wrapping_dimensions_list
    l = parse(Int64,split(dimensions,'x')[1])
    w = parse(Int64,split(dimensions,'x')[2])
    h = parse(Int64,split(dimensions,'x')[3])

    # Wrapping paper calculations
    side_areas = [l*w,l*h,w*h]
    SA = 2 * l * w + 2 * w * h + 2 * h * l
    SA += minimum(side_areas)
    surface_area += SA

    # Ribbon calculations
    side_1_perimeter = l + l + h + h
    side_2_perimeter = l + l + w + w
    side_3_perimeter = h + h + w + w
    volume = l * w * h
    min_perimeter = minimum([side_1_perimeter,side_2_perimeter,side_3_perimeter])
    length_ribbon += volume
    length_ribbon += min_perimeter
end

println("The elves should order $surface_area ft^2 of wrapping paper and $length_ribbon ft of ribbon")


text:
import Base.readline

macro asserteq(a,b)
    msg = "$a != $b"
    :($a == $b || throw(AssertionError("$($msg): $($a), $($b)")))
end

macro qv(v)
    s = string(v)
    :(($s, $v))
end

"""Return a dictionary containing the variables `vs` together with their values"""
macro qvs(vs...)
    v = (vs...)[1]
    if length(vs) == 1
        s = string(v)
        :(@qv $v)
    else
        :([@qv($v); @qvs($((vs...)[2:end]...))])
    end
end

"""Read line number `line` in file `filename`"""
function readline(filename::AbstractString, line::Int)
    open(filename) do io
        s = bytestring(Mmap.mmap(io))
        line <= 0 && error("negative line number")
        start = 0
        stop  = findfirst(s, '\n')
        i = 1
        while (i == 1 || start > 0) && i < line
            start = stop
            stop = findnext(s, '\n', stop+1)
            i += 1
        end
        line == i && return copy(s[start+1:stop-1])
        error("`$filename` only has $(i-1) lines")
    end
end

"""Return index of `vals` which is closest to `x`"""
function closest_index(vals, x)
    _, index = findmin([abs(v - x) for v in vals])
    return index
end
closest(vals, x) = vals[closest_index(vals, x)]

"""Fits a linear line through the values in `enumerate(ys)` and returns the maximum absolute deviation"""
lindev{T<:Number}(ys::Vector{T}) = maxabs(lindev_vec(ys))

"""Return deviation vector from a linear fitted line"""
function lindev_vec{T<:Number}(ys::Vector{T})
    a, b = linreg(collect(1.0:length(ys)), ys)
    return [y - (a + b*i) for (i, y) in enumerate(ys)]
end

"""Resample measurement S at retention times R such that the new matrix has size n*size(S,2)"""
function resample(S::Matrix, R::Vector, n::Int)
    @assert issorted(R)
    M = zeros(n, size(S, 2))
    Rmin, Rmax = R[1], R[size(S, 1)]
    for (i, t) = enumerate(linspace(Rmin, Rmax, n))
        i0 = findlast(R .<= t)
        i1 = findfirst(R .>= t)
        a = R[i1] - R[i0] <= 1e-10 ? 1.0 : (t - R[i1]) / (R[i0] - R[i1])
        for j = 1:size(S, 2)
            M[i,j] = a*S[i0, j] + (1 - a)*S[i1, j]
        end
    end
    return M
end

"""Uniform resampling"""
resample(S, n::Int) = resample(S, collect(1:n), n)

"""Resampling according to a function"""
resample(S::Matrix, f::Function, n::Integer) = resample(S, [f(i) for i=1:n], n)
export resample

"""Resample by given old and new R values"""
function resample(S::Matrix, oldR::Vector, newR::Vector)
    @assert issorted(oldR)
    @assert issorted(newR)
    @assert oldR[1] <= newR[1]
    @assert newR[end] <= oldR[end]
    @assert size(S, 1) == length(oldR)
    M = zeros(length(newR), size(S, 2))
    for (i, t) = enumerate(newR)
        i0 = max(findlast(oldR .<= t), 1)
        i1 = min(findfirst(oldR .>= t), length(oldR))
        a = oldR[i1] - oldR[i0] <= 1e-10 ? 1.0 : (t - oldR[i1]) / (oldR[i0] - oldR[i1])
        for j = 1:size(S, 2)
            M[i,j] = a*S[i0, j] + (1 - a)*S[i1, j]
        end
    end
    return M
end
resample(S::Matrix, oldR::Vector, f::Function, n::Integer) = resample(S, oldR, [f(i) for i = 1:n])

"""Compute the forward difference of a vector"""
function forward_diff{T}(x::Vector{T})
    return [x[i] - x[i-1] for i = 2:length(x)]
end

"""Plot a histogramm of an array"""
ploth(x; bins=300) = PyPlot.gca()[:hist](x[:], bins=bins)


text:
module Year2021

include("Day25.jl")

end


text:
function *(A::Matrix{Float64}, dist::MvSkewTDist)
    size(A, 2) == length(dist) || throw(ArgumentError("Matrix and distribution must have compatible dimensions"))
    ξ = A*dist.ξ
    Ω = PDMat(X_A_Xt(dist.Ω, A))
    distω, distΩz = scale_and_cor(dist.Ω)
    ω, Ωz =  scale_and_cor(Ω)
    B = distω\(dist.Ω*A')
    c = sqrt(1 + dot(dist.α,(distΩz.mat - X_invA_Xt(Ω, B))*dist.α))
    α = (ω * (Ω\(B'dist.α)))/c
    return MvSkewTDist(ξ, Ω, α, dist.df)
end

function *(x::Vector{Float64}, dist::MvSkewTDist)
    length(x) == length(dist) || throw(ArgumentError("Vector and distribution must have compatible dimensions"))
    ξ = dot(x, dist.ξ)
    ω = sqrt(quad(dist.Ω, x))
    distω, distΩz = scale_and_cor(dist.Ω)
    b = distω\(dist.Ω*x)
    c = sqrt(1.0 + dot(dist.α, (distΩz.mat - (b*b')/ω^2)*dist.α))
    α =  dot(b,dist.α)/(c*ω)
    SkewTDist(ξ, ω, α, dist.df)
end

function marginals(dist::MvSkewTDist, indices::Vector{Int})
    n = length(dist)
    (all(0 .< indices) && all(indices .<= n)) || throw(ArgumentError("Indices must be between 1 and $(n)"))
    A = zeros(Float64, length(indices), n)
    for (i, k) in enumerate(indices)
        A[i,k] = 1.0
    end
    return A*dist
end

function marginals(dist::MvSkewTDist, i::Int)
    n = length(dist)
    (i > 0 && i < n) || throw(ArgumentError("Indices must be between 1 and $(n)"))
    x = zeros(n)
    x[i] = 1.0
    return x * dist
end


text:
# plot first 5000 time points supplemental
include("../Model_scripts/Functions.jl");
# plot the number of infected & uninfected hosts
plot_population_numbers(
    hcat(labels...),
    "by=1.0\n \nNumber of infected and uninfected hosts",
    100,
    5000,
    "$figure_directory/supplemental/supp_1b.png"
)
# plot the average weighted R0
plot_reproductive_rate(6, 5000, "$figure_directory/supplemental/supp_1d.png")
# plot the average weighted V0
plot_vertical_reproductive_ratio(1, 5000, "$figure_directory/supplemental/supp_1f.png")
# plot the weighted average mortality rate through time
plot_mortality(1, 5000, "$figure_directory/supplemental/supp_1h.png")
# plot the evenness through time
plot_evenness(1, 5000, "$figure_directory/supplemental/supp_1j.png")


text:
#
# formats.jl -
#
# Definitions of OI-FITS tables.
#
#------------------------------------------------------------------------------

module Formats

"""
    def = OIFITS.Formats.FieldDefinition(...)

yields a field definition of a keyword or column, of an entry of an OI-FITS
data-block.   The following fields are available:

- `def.name` specifies the name of the keyword or column in the FITS file.

- `def.symb` specifies the symbolic field name in the corresponding Julia
  structure.

- `def.type` specifies the data type.  It is `:A` for strings, `:C` for
  complexes, `:D` or `:E` for reals (floating-points), `:I` or `:J` for
  integers, and `:L` for booleans.

- `def.rank` is `0` for keywords, `n > 0` for `n` elements per row, `n < 0` for
  cell of `-n` dimensions of length equal to the number of
  spectral channels.

- `def.optional` specifies whether the keyword / column is optional according
  to the standard.

- `def.descr` is a string describing the keyword / column.

These instances are built by the macros [`OIFITS.Formats.@header`](@ref) and
[`OIFITS.Formats.@column`](@ref) for a FITS keyword and column respectively.

Base method `ndims(def)` yields the dimensionality of `def`, that is `0` for a
FITS keyword, `1` for a column of scalar cells and `n ≥ 2` for a column of
wavelength-wise multi-dimensional cells.

"""
struct FieldDefinition
    name::String   # Keyword/column name in FITS file.
    symb::Symbol   # Keyword/column symbolic name in structures.
    type::Symbol   # Data type.
    rank::Int      # 0 for keywords, n > 0 for n elements per row, n < 0 for
                   # cell of -n dimensions of length equal to the number of
                   # spectral channels.
    optional::Bool # Optional field?
    descr::String  # Description.
    units::String  # Units.
end

Base.ndims(s::FieldDefinition) =
    (s.rank == 0   ? 0          : # keyword
     s.type === :A ? 1          : # column of strings
     s.rank == 1   ? 1          : # column of values
     s.rank < 0    ? 1 - s.rank : # column of wavelength-wise multi-dim. cells
     2)                           # column of multiple values

const FORMATS = Dict{Tuple{Symbol,Int},Vector{FieldDefinition}}()

"""
    OIFITS.get_format(ext, revn; throw_errors=false) -> A

yields a vector of `OIFITS.Formats.Field` instances corresponding to the FITS
keywords and columns of an OI-FITS data-block `ext` with revision number
`revn`.  If `ext` and `revn` do not correspond to any known definition,
`nothing` is returned if `throw_errors` is `false` (the default), otherwise an
exception is thrown.

For example:

    for def in OIFITS.get_format(:OI_VIS2, 2)
        println(
            (def.rank == 0 ? "Keyword" : "Column "), " \"",
            def.name, "\" => `", def.symb, "` -- ", def.descr)
    end

"""
function get_format(ext::Symbol, rev::Integer; throw_errors::Bool=false)
    spec = get(FORMATS, (ext,Int(rev)), nothing)
    if spec === nothing && throw_errors
        error("unknown revision ", rev, " of OI-FITS extension \"", ext, "\"")
    end
    return spec
end

"""
    @header key type descr

yields an instance of `FieldDefinition` for an OI-FITS header keyword `key` (a
string) with value type and description specified by `type` and `descr`.  The
syntax closely follows that of the tables of the papers specifying the OI-FITS
format.

References:

- G. Duvert, J. S. Young & Ch. A. Hummel: "*OIFITS 2: the 2nd version of the
  data exchange standard for optical interferometry*", Astron. & Astrophys.,
  **597**, A8 (2017).

- T. A. Pauls, J. S. Young, W. D. Cotton & J. D. Monnier: "*A Data Exchange
  Standard for Optical (Visible/IR) Interferometry*", PASP **117**, 1255-1262
  (2005).

"""
macro header(key, type, str)
    opt, sym, rank = parse_type(type)
    isa(key, String) || error("expecting a string for keyword name, got ", key)
    rank == 0 || error("invalid dimensions in `", type, "` for header keyword")
    descr, units = parse_description(str)
    quote
        FieldDefinition($key, $(QuoteNode(to_field_name(key))),
                        $(QuoteNode(sym)), $rank, $opt, $descr, $units)
    end
end

"""
    @column key type descr

yields an instance of `FieldDefinition` for an OI-FITS column `key` (a string
or symbol) with cell type and dimensions specified by `type` and description
given by `descr`.  The syntax closely follows that of the tables of the papers
specifying the OI-FITS format.

References:

- G. Duvert, J. S. Young & Ch. A. Hummel: "*OIFITS 2: the 2nd version of the
  data exchange standard for optical interferometry*", Astron. & Astrophys.,
  **597**, A8 (2017).

- T. A. Pauls, J. S. Young, W. D. Cotton & J. D. Monnier: "*A Data Exchange
  Standard for Optical (Visible/IR) Interferometry*", PASP **117**, 1255-1262
  (2005).

"""
macro column(key, type, str)
    opt, sym, rank = parse_type(type)
    isa(key, String) || error("expecting a string for column name, got ", key)
    rank != 0 || error("invalid dimensions in `", type, "` for column")
    descr, units = parse_description(str)
    quote
        FieldDefinition($key, $(QuoteNode(to_field_name(key))),
                        $(QuoteNode(sym)), $rank, $opt, $descr, $units)
    end
end

to_field_name(str::AbstractString) =
    Symbol(lowercase(replace(replace(str, r"^OI_" => ""), '-' => '_')))

"""
    parse_description(str) -> (descr, units)

yields a 2-tuple of strings extracted from string `str` which shall have the
form:

    str = "\$descr [\$units]"

with any leading/trailing spaces and where the units part is optional
(returned as an empty string if not specified).

"""
function parse_description(str::String)
    txt = strip(str)
    if endswith(txt, ']')
        i = findlast('[', txt)
        if i !== nothing
            units = String(txt[nextind(txt,i):prevind(txt,end)])
            i_first = firstindex(txt)
            if i > i_first
                return (String(strip(txt[i_first:prevind(txt,i)])), units)
            else
                return ("", units)
            end
        end
    end
    return (String(txt), "")
end

parse_description(x::Any) = error("description must be a string")

"""
     parse_type(x) -> (opt, sym, rank)

parses type `x` in OI-FITS format specification and yields a 3-tuple
`(opt,type,dims)` where `opt` is a boolean indicating whether the field is
optional, `sym` is a single letter symbol corresponding to the type of the
field and `rank` specifies the dimensions.

Argument `x` is a symbol or an expression.

""" parse_type

const APOSTROPHE = Symbol(Char(39))

# 1st pass is to figure out the optional flag (an apostrophe).
parse_type(x::Symbol) = parse_type(false, x)
parse_type(x::Expr) = begin
    if x.head === APOSTROPHE
        parse_type(true, x.args[1])
    else
        parse_type(false, x)
    end
end

# 2nd pass is to parse type letter and dimensions.
parse_type(opt::Bool, x::Symbol) = (opt, parse_type_letter(x), 0)
parse_type(opt::Bool, x::Expr) = begin
    if x.head === :call && length(x.args) ≥ 1
        symb = parse_type_letter(x.args[1])
        rank = parse_type_dims(x.args[2:end]...)
        return (opt, symb, rank)
    else
        invalid_type(opt, x)
    end
end
parse_type(opt::Bool, x::Any) = invalid_type(opt, x)

parse_type_letter(x::Symbol) = (x ∈ (:A, :C, :D, :E, :I, :J, :L) ? x :
                                error("unknown type letter/symbol: ", x))

parse_type_dims(dims...) = invalid_type_dims(dims)
parse_type_dims(dim::Integer) =
    (dim ≥ 1 ? Int(dim) : invalid_type_dims((dim,)))
parse_type_dims(dim::Symbol) =
    (dim === :W ? -1 : invalid_type_dims((dim,)))
parse_type_dims(dim1::Symbol, dim2::Symbol) =
    (dim1 === :W && dim2 === :W ? -2 : invalid_type_dims((dim1,dim2)))

@noinline invalid_type_dims(dims) = error("invalid type dimensions: ", dims)

@noinline invalid_type(opt::Bool, x) =
    error("invalid type specification: ", x, (opt ? APOSTROPHE : nothing))

function define(db::AbstractString,
                numbs::Tuple{Vararg{Integer}},
                specs::AbstractVector{FieldDefinition})
    for revn in numbs
        define(Symbol(db), revn, specs)
    end
end

define(db::AbstractString, revn::Integer, specs::AbstractVector{FieldDefinition}) =
    define(Symbol(db), revn, specs)

function define(db::Symbol, revn::Integer, specs::AbstractVector{FieldDefinition})
    FORMATS[(db,Int(revn))] = specs
end

# OI_TARGET definition (1st revision):
define("OI_TARGET", 1, [
    @header "OI_REVN"    I      "revision number of the table definition"
    #---------------------------------------------------------
    @column "TARGET_ID"  I(1)   "index number"
    @column "TARGET"     A(16)  "target name"
    @column "RAEP0"      D(1)   "RA at mean equinox [deg]"
    @column "DECEP0"     D(1)   "DEC at mean equinox [deg]"
    @column "EQUINOX"    E(1)   "equinox [yr]"
    @column "RA_ERR"     D(1)   "error in RA at mean equinox [deg]"
    @column "DEC_ERR"    D(1)   "error in DEC at mean equino [deg]"
    @column "SYSVEL"     D(1)   "systemic radial velocity [m/s]"
    @column "VELTYP"     A(8)   "reference for radial velocity"
    @column "VELDEF"     A(8)   "definition of radial velocity"
    @column "PMRA"       D(1)   "proper motion in RA [deg/yr]"
    @column "PMDEC"      D(1)   "proper motion in DEC [deg/yr]"
    @column "PMRA_ERR"   D(1)   "error of proper motion in RA [deg/yr]"
    @column "PMDEC_ERR"  D(1)   "error of proper motion in DEC [deg/yr]"
    @column "PARALLAX"   E(1)   "parallax [deg]"
    @column "PARA_ERR"   E(1)   "error in parallax [deg]"
    @column "SPECTYP"    A(16)  "spectral type"
])

# OI_TARGET definition (2nd revision):
define("OI_TARGET", 2, [
    @header "OI_REVN"    I      "revision number of the table definition"
    #---------------------------------------------------------
    @column "TARGET_ID"  I(1)   "index number"
    @column "TARGET"     A(16)  "target name"
    @column "RAEP0"      D(1)   "RA at mean equinox [deg]"
    @column "DECEP0"     D(1)   "DEC at mean equinox [deg]"
    @column "EQUINOX"    E(1)   "equinox [yr]"
    @column "RA_ERR"     D(1)   "error in RA at mean equinox [deg]"
    @column "DEC_ERR"    D(1)   "error in DEC at mean equino [deg]"
    @column "SYSVEL"     D(1)   "systemic radial velocity [m/s]"
    @column "VELTYP"     A(8)   "reference for radial velocity"
    @column "VELDEF"     A(8)   "definition of radial velocity"
    @column "PMRA"       D(1)   "proper motion in RA [deg/yr]"
    @column "PMDEC"      D(1)   "proper motion in DEC [deg/yr]"
    @column "PMRA_ERR"   D(1)   "error of proper motion in RA [deg/yr]"
    @column "PMDEC_ERR"  D(1)   "error of proper motion in DEC [deg/yr]"
    @column "PARALLAX"   E(1)   "parallax [deg]"
    @column "PARA_ERR"   E(1)   "error in parallax [deg]"
    @column "SPECTYP"    A(16)  "spectral type"
    @column "CATEGORY"   A(3)'  "CALibrator or SCIence target"
])

# OI_ARRAY definition (1st revision):
define("OI_ARRAY", 1, [
    @header "OI_REVN"    I      "revision number of the table definition"
    @header "ARRNAME"    A      "array name for cross-referencing"
    @header "FRAME"      A      "coordinate frame"
    @header "ARRAYX"     D      "array center X-coordinate [m]"
    @header "ARRAYY"     D      "array center Y-coordinate [m]"
    @header "ARRAYZ"     D      "array center Z-coordinate [m]"
    #------------------------------------------------------------
    @column "TEL_NAME"   A(16)  "telescope name"
    @column "STA_NAME"   A(16)  "station name"
    @column "STA_INDEX"  I(1)   "station index"
    @column "DIAMETER"   E(1)   "element diameter [m]"
    @column "STAXYZ"     D(3)   "station coordinates relative to array center [m]"
])

# OI_ARRAY definition (2nd revision):
define("OI_ARRAY", 2, [
    @header "OI_REVN"    I      "revision number of the table definition"
    @header "ARRNAME"    A      "array name for cross-referencing"
    @header "FRAME"      A      "coordinate frame"
    @header "ARRAYX"     D      "array center X-coordinate [m]"
    @header "ARRAYY"     D      "array center Y-coordinate [m]"
    @header "ARRAYZ"     D      "array center Z-coordinate [m]"
    #------------------------------------------------------------
    @column "TEL_NAME"   A(16)  "telescope name"
    @column "STA_NAME"   A(16)  "station name"
    @column "STA_INDEX"  I(1)   "station index"
    @column "DIAMETER"   E(1)   "element diameter [m]"
    @column "STAXYZ"     D(3)   "station coordinates relative to array center [m]"
    @column "FOV"        D(1)   "photometric field of view [arcsec]"
    @column "FOVTYPE"    A(6)   "model for FOV: 'FWHM' or 'RADIUS'"
])

# OI_WAVELENGTH definition (1st and 2nd revisions):
define("OI_WAVELENGTH", (1, 2), [
    @header "OI_REVN"    I    "revision number of the table definition"
    @header "INSNAME"    A    "name of detector for cross-referencing"
    #------------------------------------------------------------
    @column "EFF_WAVE"  E(1)  "effective wavelength of channel [m]"
    @column "EFF_BAND"  E(1)  "effective bandpass of channel [m]"
])

# OI_VIS definition (1st revision):
define("OI_VIS", 1, [
    @header "OI_REVN"    I     "revision number of the table definition"
    @header "DATE-OBS"   A     "UTC start date of observations"
    @header "ARRNAME"    A'    "name of corresponding array"
    @header "INSNAME"    A     "name of corresponding detector"
    #------------------------------------------------------------
    @column "TARGET_ID"  I(1)  "target number as index into OI_TARGET table"
    @column "TIME"       D(1)  "UTC time of observation [s]"
    @column "MJD"        D(1)  "modified Julian Day [day]"
    @column "INT_TIME"   D(1)  "integration time [s]"
    @column "VISAMP"     D(W)  "visibility amplitude"
    @column "VISAMPERR"  D(W)  "error in visibility amplitude"
    @column "VISPHI"     D(W)  "visibility phase [deg]"
    @column "VISPHIERR"  D(W)  "error in visibility phase [deg]"
    @column "UCOORD"     D(1)  "U coordinate of the data [m]"
    @column "VCOORD"     D(1)  "V coordinate of the data [m]"
    @column "STA_INDEX"  I(2)  "station numbers contributing to the data"
    @column "FLAG"       L(W)  "flag"
])

# FIXME: In OI-FITS Rev. 2 only MJD and DATE-OBS must be used to express
#        time. The TIME column is retained only for backwards compatibility
#        and must contain zeros.

# OI_VIS definition (2nd revision):
define("OI_VIS", 2, [
    @header "OI_REVN"          I        "revision number of the table definition"
    @header "DATE-OBS"         A        "UTC start date of observations"
    @header "ARRNAME"          A        "name of corresponding array"
    @header "INSNAME"          A        "name of corresponding detector"
    @header "CORRNAME"         A'       "name of corresponding correlation table"
    @header "AMPTYP"           A'       "'ABSOLUTE', 'DIFFERENTIAL' or 'CORRELATED FLUX'"
    @header "PHITYP"           A'       "'ABSOLUTE', or 'DIFFERENTIAL'"
    @header "AMPORDER"         I'       "polynomial fit order for differential amplitudes"
    @header "PHIORDER"         I'       "polynomial fit order for differential phases"
    #------------------------------------------------------------------------
    @column "TARGET_ID"        I(1)     "target number as index into OI_TARGET table"
    @column "TIME"             D(1)     "UTC time of observation [s]"
    @column "MJD"              D(1)     "modified Julian Day [day]"
    @column "INT_TIME"         D(1)     "integration time [s]"
    @column "VISAMP"           D(W)     "visibility amplitude"
    @column "VISAMPERR"        D(W)     "error in visibility amplitude"
    @column "CORRINDX_VISAMP"  J(1)'    "index into correlation matrix for 1st VISAMP element"
    @column "VISPHI"           D(W)     "visibility phase [deg]"
    @column "VISPHIERR"        D(W)     "error in visibility phase [deg]"
    @column "CORRINDX_VISPHI"  J(1)'    "index into correlation matrix for 1st VISPHI element"
    @column "VISREFMAP"        L(W,W)'  "true where spectral channels were taken as reference for differential visibility computation"
    @column "RVIS"             D(W)'    "real part of complex coherent flux"
    @column "RVISERR"          D(W)'    "error on RVIS"
    @column "CORRINDX_RVIS"    J(1)'    "index into correlation matrix for 1st RVIS element"
    @column "IVIS"             D(W)'    "imaginary part of complex coherent flux"
    @column "IVISERR"          D(W)'    "error on IVIS"
    @column "CORRINDX_IVIS"    J(1)'    "index into correlation matrix for 1st IVIS element"
    @column "UCOORD"           D(1)     "U coordinate of the data [m]"
    @column "VCOORD"           D(1)     "V coordinate of the data [m]"
    @column "STA_INDEX"        I(2)     "station numbers contributing to the data"
    @column "FLAG"             L(W)     "flag"
])

# OI_VIS2 definition (1st revision):
define("OI_VIS2", 1, [
    @header "OI_REVN"    I     "revision number of the table definition"
    @header "DATE-OBS"   A     "UTC start date of observations"
    @header "ARRNAME"    A'    "name of corresponding array"
    @header "INSNAME"    A     "name of corresponding detector"
    #------------------------------------------------------------
    @column "TARGET_ID"  I(1)  "target number as index into OI_TARGET table"
    @column "TIME"       D(1)  "UTC time of observation [s]"
    @column "MJD"        D(1)  "modified Julian Day [day]"
    @column "INT_TIME"   D(1)  "integration time [s]"
    @column "VIS2DATA"   D(W)  "squared visibility"
    @column "VIS2ERR"    D(W)  "error in squared visibility"
    @column "UCOORD"     D(1)  "U coordinate of the data [m]"
    @column "VCOORD"     D(1)  "V coordinate of the data [m]"
    @column "STA_INDEX"  I(2)  "station numbers contributing to the data"
    @column "FLAG"       L(W)  "flag"
])

# OI_VIS2 definition (2nd revision):
define("OI_VIS2", 2, [
    @header "OI_REVN"            I      "revision number of the table definition"
    @header "DATE-OBS"           A      "UTC start date of observations"
    @header "ARRNAME"            A      "name of corresponding array"
    @header "INSNAME"            A      "name of corresponding detector"
    @header "CORRNAME"           A'     "name of corresponding correlation table"
    #------------------------------------------------------------
    @column "TARGET_ID"          I(1)   "target number as index into OI_TARGET table"
    @column "TIME"               D(1)   "UTC time of observation [s]"
    @column "MJD"                D(1)   "modified Julian Day [day]"
    @column "INT_TIME"           D(1)   "integration time [s]"
    @column "VIS2DATA"           D(W)   "squared visibility"
    @column "VIS2ERR"            D(W)   "error in squared visibility"
    @column "CORRINDX_VIS2DATA"  J(1)'  "index into correlation matrix for 1st VIS2DATA element"
    @column "UCOORD"             D(1)   "U coordinate of the data [m]"
    @column "VCOORD"             D(1)   "V coordinate of the data [m]"
    @column "STA_INDEX"          I(2)   "station numbers contributing to the data"
    @column "FLAG"               L(W)   "flag"
])

# OI_T3 definition (1st revision):
define("OI_T3", 1, [
    @header "OI_REVN"    I     "revision number of the table definition"
    @header "DATE-OBS"   A     "UTC start date of observations"
    @header "ARRNAME"    A'    "name of corresponding array"
    @header "INSNAME"    A     "name of corresponding detector"
    #------------------------------------------------------------
    @column "TARGET_ID"  I(1)  "target number as index into OI_TARGET table"
    @column "TIME"       D(1)  "UTC time of observation [s]"
    @column "MJD"        D(1)  "modified Julian Day [day]"
    @column "INT_TIME"   D(1)  "integration time [s]"
    @column "T3AMP"      D(W)  "triple product amplitude"
    @column "T3AMPERR"   D(W)  "error in triple product amplitude"
    @column "T3PHI"      D(W)  "triple product phase [deg]"
    @column "T3PHIERR"   D(W)  "error in triple product phase [deg]"
    @column "U1COORD"    D(1)  "U coordinate of baseline AB of the triangle [m]"
    @column "V1COORD"    D(1)  "V coordinate of baseline AB of the triangle [m]"
    @column "U2COORD"    D(1)  "U coordinate of baseline BC of the triangle [m]"
    @column "V2COORD"    D(1)  "V coordinate of baseline BC of the triangle [m]"
    @column "STA_INDEX"  I(3)  "station numbers contributing to the data"
    @column "FLAG"       L(W)  "flag"
])

# OI_T3 definition (2nd revision):
define("OI_T3", 2, [
    @header "OI_REVN"         I      "revision number of the table definition"
    @header "DATE-OBS"        A      "UTC start date of observations"
    @header "ARRNAME"         A      "name of corresponding array"
    @header "INSNAME"         A      "name of corresponding detector"
    @header "CORRNAME"        A'     "name of corresponding correlation table"
    #------------------------------------------------------------
    @column "TARGET_ID"       I(1)   "target number as index into OI_TARGET table"
    @column "TIME"            D(1)   "UTC time of observation [s]"
    @column "MJD"             D(1)   "modified Julian Day [day]"
    @column "INT_TIME"        D(1)   "integration time [s]"
    @column "T3AMP"           D(W)   "triple product amplitude"
    @column "T3AMPERR"        D(W)   "error in triple product amplitude"
    @column "CORRINDX_T3AMP"  J(1)'  "index into correlation matrix for 1st T3AMP element"
    @column "T3PHI"           D(W)   "triple product phase [deg]"
    @column "T3PHIERR"        D(W)   "error in triple product phase [deg]"
    @column "CORRINDX_T3PHI"  J(1)'  "index into correlation matrix for 1st T3PHI element"
    @column "U1COORD"         D(1)   "U coordinate of baseline AB of the triangle [m]"
    @column "V1COORD"         D(1)   "V coordinate of baseline AB of the triangle [m]"
    @column "U2COORD"         D(1)   "U coordinate of baseline BC of the triangle [m]"
    @column "V2COORD"         D(1)   "V coordinate of baseline BC of the triangle [m]"
    @column "STA_INDEX"       I(3)   "station numbers contributing to the data"
    @column "FLAG"            L(W)   "flag"
])

# OI_FLUX definition (1st revision):
define("OI_FLUX", 1, [
    @header "OI_REVN"            I      "revision number of the table definition"
    @header "DATE-OBS"           A      "UTC start date of observations"
    @header "INSNAME"            A      "name of corresponding detector"
    @header "ARRNAME"            A'     "name of corresponding array"
    @header "CORRNAME"           A'     "name of corresponding correlation table"
    @header "FOV"                D'     "area of sky over which flux is integrated [arcsec]"
    @header "FOVTYPE"            A'     "model for FOV: 'FWHM' or 'RADIUS'"
    @header "CALSTAT"            A      "'C': spectrum is calibrated, 'U': uncalibrated"
    #------------------------------------------------------------
    @column "TARGET_ID"          I(1)   "target number as index into OI_TARGET table"
    @column "MJD"                D(1)   "modified Julian Day [day]"
    @column "INT_TIME"           D(1)   "integration time [s]"
    @column "FLUXDATA"           D(W)   "flux"
    @column "FLUXERR"            D(W)   "flux error"
    @column "CORRINDX_FLUXDATA"  J(1)'  "index into correlation matrix for 1st FLUXDATA element"
    @column "STA_INDEX"          I(1)'  "station number contributing to the data"
    @column "FLAG"               L(W)   "flag"
])

# OI_CORR definition (1st revision):
define("OI_CORR", 1, [
    @header "OI_REVN"   I     "revision number of the table definition"
    @header "CORRNAME"  A     "name of correlation data set"
    @header "NDATA"     I     "number of correlated data"
    #------------------------------------------------------------
    @column "IINDX"     J(1)  "1st index of correlation matrix element"
    @column "JINDX"     J(1)  "2nd index of correlation matrix element"
    @column "CORR"      D(1)  "matrix element"
])

# OI_INSPOL definition (1st revision):
define("OI_INSPOL", 1, [
    @header "OI_REVN"    I      "revision number of the table definition"
    @header "DATE-OBS"   A      "UTC start date of observations"
    @header "NPOL"       I      "number of polarization types in this table"
    @header "ARRNAME"    A      "identifies corresponding OI_ARRAY"
    @header "ORIENT"     A      "orientation of the Jones Matrix: 'NORTH' (for on-sky orientation), or 'LABORATORY'"
    @header "MODEL"      A      "describe the way the Jones matrix is estimated"
    #---------------------------------------------------------------------------------------
    @column "TARGET_ID"  I(1)   "target number as index into OI_TARGET table"
    @column "INSNAME"    A(16)  "INSNAME of this polarization"
    @column "MJD_OBS"    D(1)   "modified Julian day, start of time lapse"
    @column "MJD_END"    D(1)   "modified Julian day, end of time lapse"
    @column "JXX"        C(W)   "complex Jones Matrix component along X axis"
    @column "JYY"        C(W)   "complex Jones Matrix component along Y axis"
    @column "JXY"        C(W)   "complex Jones Matrix component between Y and X axis"
    @column "JYX"        C(W)   "complex Jones Matrix component between Y and X axis"
    @column "STA_INDEX"  I(1)'  "station number for the above matrices"
])

end # module


text:
using Plots;
using Random;
using LinearAlgebra;
pyplot()

#===============================================#
#                                               #
#                Plots Settings                 #
#                                               #
#===============================================#

upscale = 2                             # 8x upscaling in resolution
fntsm = Plots.font("sans-serif", 10.0*upscale)
fntlg = Plots.font("sans-serif", 14.0*upscale)
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
default(size=(800*upscale,600*upscale)) # Plot canvas size
default(dpi=300)                        # Only for PyPlot - presently broken

#===============================================#
#                                               #
#                  Definations                  #
#                                               #
#===============================================#

CORE_BIAS = 3
SAMPLE_SIZE = 50
FEATURE_DIM = 2
NUM_CLASSES = 2
λ = 0.1             # normalize
α = 0.001           # learning rate

#===============================================#
#                                               #
#                Tools functions                #
#                                               #
#===============================================#

function generate_dataset(sample_size::Int, feature_dim::Int=FEATURE_DIM, num_classes::Int=NUM_CLASSES)
    # random get
    Y = rand(1:num_classes,sample_size)
    # make features
    X = randn(sample_size,feature_dim) .+ (CORE_BIAS .* Y)
    # make label
    Y = Int.([y == class for y in Y, class = 1:num_classes])
    return hcat(ones(SAMPLE_SIZE,1), X) ,Y
end

# Gradient Descent
# π: predict(θ,X)
# ∇: gradient(π(θ,X),X,Y)
# loss: loss(π(θ,X),X,Y)
function gradient_descent(X,Y,π::Function,∇::Function,loss::Function,θ::Array,λ = 0)
    LR_DECAY = 0.1                  # lr_decay
    MOMENTUM = 0.9                  # Momentum
    i = 0

    P = π(θ,X)
    losses = Float64[loss(P,X,Y)]
    v = zeros(size(θ))
    while losses[end] > 0.1
        p = π(θ,X)
        Δ = ∇(p,X,Y) + λ * θ ./ length(Y)        # normalize, if λ = 0, equal to no normalization
        lr_i = α / (1 + LR_DECAY * i)
        v = MOMENTUM * v - lr_i * Δ
        θ += v
        push!(losses,loss(p,X,Y) + λ * sum(θ .* θ) ./ 2length(Y))
        println(losses[end])
    end
    return θ,losses
end

function gradient_descent(X,Y,π::Function,∇::Function,loss::Function,λ = 0)
    θ = ones(size(X,2),size(Y,2))
    return gradient_descent(X,Y,π,∇,loss,θ,λ)
end

function newtown(X,Y,π::Function,loss::Function,θ::Array,λ = 0)
    LR_DECAY = 0                      # lr_decay
    MOMENTUM = 0                      # Momentum
    i = 0

    P = π(θ,X)
    losses = Float64[loss(P,X,Y)]
    v = zeros(size(θ))
    while losses[end] > 0.1
        P = π(θ,X)
        A = Diagonal((P .* (1 .- P))[:,])
        U = X' * (Y .- P)
        H = X' * A * X
        lr_i = α / (1 + LR_DECAY * i)
        v = MOMENTUM * v + lr_i * (H^-1 * U - λ * θ ./ length(Y))
        θ = θ + v
        i = i + 1
        push!(losses,loss(P,X,Y) + λ * sum(θ .* θ) ./ 2length(Y))
        println(losses[end])
    end
    return θ,losses
end

function newtown(X,Y,π::Function,loss::Function,λ = 0)
    θ = ones(size(X,2),size(Y,2))
    return newtown(X,Y,π,loss,θ,λ)
end

function draw(θ,X,Y,losses,name::String)
    # draw the line and points
    x1 = Array(1:1:7)
    x2 = (- θ[1,1] .- θ[2,1] .* x1) ./ θ[3,1]
    colors = [[:red, :green][l+1] for l in Y[:,1]]
    scatter(X[:,2], X[:,3],  markercolor = colors)
    plot!(x1, x2)
    savefig("$name-plot")
    # losses plot
    plot(1:size(losses,1),losses[1:end])
    savefig("$name-loss")
end

#===============================================#
#                                               #
#               Logistic Regression             #
#                                               #
#===============================================#

# generate data
X,Y = generate_dataset(SAMPLE_SIZE,FEATURE_DIM,2)
Y = Y[:,1]

# class model
sigmoid(θ,X) = 1 ./ (1 .+ exp.(- X * θ))

# gradient
# ∇(sigmoid(θ,X),X,Y)
∇(π,X,Y) = X' * (π - Y) ./ size(Y,1)

# loss function
# loss(sigmoid(θ,X),X,Y)
loss(π,X,Y) = - sum(Y .* log.(π) + (1 .- Y) .* log.(1 .- π)) ./ length(Y)

θ, losses = gradient_descent(X,Y,sigmoid,∇,loss)
draw(θ,X,Y,losses,"grad")

θ, losses = newtown(X,Y,sigmoid,loss)
draw(θ,X,Y,losses,"newtown")

#===============================================#
#                                               #
#                    Test  UCI                  #
#                                               #
#===============================================#

using DataFrames, CSV;

function toLabel(labels)
    classes = unique(labels)
    function λ(label)
        for (index, value) in enumerate(classes)
            if label == value
                return index
            end
        end
    end
    return Int.([y == class for y in λ.(labels), class = 1:size(classes,1)]),classes
end

datasets = CSV.read("data/iris-cut.csv");
labels = Array{String}(copy(datasets[:Species]))

X = Matrix{Float64}(datasets[:,1:end-1])
X = hcat(ones(size(X,1),1), X)
Y,classes = toLabel(labels)
Y = Y[:,1]

θ, losses = gradient_descent(X,Y,sigmoid,∇,loss)
draw(θ,X,Y,losses,"iris")

#===============================================#
#                                               #
#                Softmax Regression             #
#                                               #
#===============================================#
using DataFrames, CSV;

function toLabel(labels)
    classes = unique(labels)
    function λ(label)
        for (index, value) in enumerate(classes)
            if label == value
                return index
            end
        end
    end
    return Int.([y == class for y in λ.(labels), class = 1:size(classes,1)]),classes
end

# class model
function softmax(θ,X)
    Y = exp.(X * θ)
    total = sum(Y,dims = 2)
    return Y ./ total
end

# gradient
# ∇(softmax(θ,X),X,Y)
∇(π,X,Y) = X' * (π - Y) ./ length(Y)

# loss function
# loss(softmax(θ,X),X,Y)
loss(π,X,Y) = - sum(Y .* log.(π)) ./ length(Y)

datasets = CSV.read("data/iris.csv");
labels = Array{String}(copy(datasets[:Species]))

# generate data
X = Matrix{Float64}(datasets[:,1:end-1])
# set bias row
X = hcat(ones(size(X,1),1), X)
Y,classes = toLabel(labels)

θ, losses = gradient_descent(X,Y,softmax,∇,loss)
θ, losses = gradient_descent(X,Y,softmax,∇,loss,θ)

# softmax
softmax(θ,X)

# losses plot
plot(1:size(losses,1),losses)

max_stars_repo_name: AlongWY/HIT-WorkShop

max_issues_repo_name: AlongWY/HIT-WorkShop

max_forks_repo_name: AlongWY/HIT-WorkShop


text:
using QXTns
using QXTools
using QXTools.Circuits

@testset "Test circuit to Tensor Network Circuit conversion" begin
    circ = create_test_circuit()

    # check tensornetwork size when no input or output
    tnc = convert_to_tnc(circ, no_input=true, no_output=true, decompose=false)
    @test QXTns.qubits(tnc) == 3
    @test length(tnc) == 3
    # make sure there are gates on all qubits as expected
    @test all(tnc.input_indices .!= tnc.output_indices)

    tnc = convert_to_tnc(circ, no_input=true, no_output=true, decompose=true)
    @test qubits(tnc) == 3
    @test length(tnc) == 5
    # make sure there are gates on all qubits as expected
    @test all(tnc.input_indices .!= tnc.output_indices)

    # convert again adding input
    tnc = convert_to_tnc(circ, no_input=false, no_output=true, decompose=false)
    @test length(tnc) == 6

    # check data of first tensor matches that of first gate
    @test tensor_data(tnc, first(keys(tnc))) == gate_matrix(first(circ.circ_ops))
end

@testset "Test simple circuit contraction to verify conversion to circuit working" begin
    circ = create_test_circuit()
    # we add input but no output to get full output vector
    tnc = convert_to_tnc(circ, no_input=false, no_output=true)
    tnc.tn.tensor_map
    # disable_hyperindices!(tnc.tn)
    output = QXTns.simple_contraction(tnc.tn); output = reshape(output, prod(size(output)))

    ref = zeros(8)
    ref[[1,8]] .= 1/sqrt(2)
    @test all(output .≈ ref)

    tnc = convert_to_tnc(circ, no_input=false, no_output=true, decompose=false)

    output = simple_contraction(tnc)
    output = reshape(output, prod(size(output)))
    ref = zeros(8)
    ref[[1,8]] .= 1/sqrt(2)
    @test all(output .≈ ref)
end


text:
#####
##### Validation
#####

# Check invariants between the CacheManager and the managed heaps.
function check(manager::CacheManager)
    @spinlock alloc_lock(manager) maybe_cleanup!(manager)

    passed = true
    if !check(manager.remote_heap)
        println("Remote Heap failed!")
        passed = false
    end

    if !check(manager.local_heap)
        println("Local Heap failed!")
        passed = false
    end

    #####
    ##### Blocks allocated by heaps
    #####

    remote_ids = Set{UInt64}()
    remote_size_allocated = 0
    for block in manager.remote_heap
        if !isfree(block)
            push!(remote_ids, CachedArrays.getid(block))
            remote_size_allocated += length(block)
        end
    end

    local_ids = Set{UInt64}()
    local_size_allocated = 0
    for block in manager.local_heap
        if !isfree(block)
            push!(local_ids, CachedArrays.getid(block))
            local_size_allocated += length(block)
        end
    end

    #####
    ##### Blocks visible from the manager
    #####

    manager_local_ids = Set{UInt64}()
    manager_remote_ids = Set{UInt64}()

    for (id, (_, backedge)) in manager.map.dict
        block = unsafe_block(unsafe_load(backedge))

        id = getid(block)
        pool = getpool(block)

        set = pool == Local ? manager_local_ids : manager_remote_ids
        other = pool == Local ? manager_remote_ids : manager_local_ids
        push!(set, id)
        sibling = getsibling(block)
        if sibling !== nothing
            # Sanity checks.
            if getsibling(sibling) !== block
                println("Non-symmetrical linking between the following siblings:")
                println(block)
                println(sibling)
                passed = false
            end

            if getid(sibling) != id
                println("Non-symmetrical id between the following siblings:")
                println(block)
                println(sibling)
                passed = false
            end

            if getpool(sibling) == pool
                println("Non-symmetrical pool between the following siblings:")
                println(block)
                println(sibling)
                passed = false
            end

            if length(sibling) != length(block)
                println("Non-symmetrical length between the following siblings:")
                println(block)
                println(sibling)
                passed = false
            end
            push!(other, id)
        end
    end

    if !issubset(manager_local_ids, local_ids)
        println("Heap local IDs not a subset of Manager local ids!")
        passed = false
    end

    if !issubset(local_ids, manager_local_ids)
        println("Manager local ids not a subset of heap ids!")
        passed = false
    end

    if !issubset(manager_remote_ids, remote_ids)
        println("Heap remote IDs not a subset of Manager remote ids!")
        passed = false
    end

    if !issubset(remote_ids, manager_remote_ids)
        println("Manager remote ids not a subset of heap ids!")
        @show setdiff(remote_ids, manager_remote_ids)
        @show setdiff(manager_remote_ids, remote_ids)
        passed = false
    end

    return passed
end