text:
"""
$(DocStringExtensions.README)
"""
module ModelingToolkit
using DocStringExtensions
using AbstractTrees
using DiffEqBase, SciMLBase, Reexport
using Distributed
using StaticArrays, LinearAlgebra, SparseArrays, LabelledArrays
using InteractiveUtils
using Latexify, Unitful, ArrayInterface
using MacroTools
@reexport using UnPack
using Setfield, ConstructionBase
using DiffEqJump
using DataStructures
using SpecialFunctions, NaNMath
using RuntimeGeneratedFunctions
using Base.Threads
using DiffEqCallbacks
import MacroTools: splitdef, combinedef, postwalk, striplines
import Libdl
using DocStringExtensions
using Base: RefValue
import IfElse

import Distributions

RuntimeGeneratedFunctions.init(@__MODULE__)

using RecursiveArrayTools

import SymbolicUtils
import SymbolicUtils: istree, arguments, operation, similarterm, promote_symtype,
                      Symbolic, Term, Add, Mul, Pow, Sym, FnType,
                      @rule, Rewriters, substitute, metadata
using SymbolicUtils.Code
import SymbolicUtils.Code: toexpr
import SymbolicUtils.Rewriters: Chain, Postwalk, Prewalk, Fixpoint
import JuliaFormatter

using Reexport
@reexport using Symbolics
export @derivatives
using Symbolics: _parse_vars, value, @derivatives, get_variables,
                 exprs_occur_in, solve_for, build_expr, unwrap, wrap,
                 VariableSource, getname, variable
import Symbolics: rename, get_variables!, _solve, hessian_sparsity,
                  jacobian_sparsity, islinear, _iszero, _isone,
                  tosymbol, lower_varname, diff2term, var_from_nested_derivative,
                  BuildTargets, JuliaTarget, StanTarget, CTarget, MATLABTarget,
                  ParallelForm, SerialForm, MultithreadedForm, build_function,
                  rhss, lhss, prettify_expr, gradient,
                  jacobian, hessian, derivative, sparsejacobian, sparsehessian,
                  substituter, scalarize, getparent

import DiffEqBase: @add_kwonly

import LightGraphs: SimpleDiGraph, add_edge!, incidence_matrix

using Requires

for fun in [:toexpr]
    @eval begin
        function $fun(eq::Equation; kw...)
            Expr(:(=), $fun(eq.lhs; kw...), $fun(eq.rhs; kw...))
        end

        $fun(eqs::AbstractArray; kw...) = map(eq->$fun(eq; kw...), eqs)
        $fun(x::Integer; kw...) = x
        $fun(x::AbstractFloat; kw...) = x
    end
end

"""
$(TYPEDEF)

TODO
"""
abstract type AbstractSystem end
abstract type AbstractTimeDependentSystem <: AbstractSystem end
abstract type AbstractTimeIndependentSystem <: AbstractSystem end
abstract type AbstractODESystem <: AbstractTimeDependentSystem end
abstract type AbstractMultivariateSystem <: AbstractSystem end

"""
$(TYPEDSIGNATURES)

Get the set of independent variables for the given system.
"""
function independent_variables end

function independent_variable end

"""
$(TYPEDSIGNATURES)

Get the set of states for the given system.
"""
function states end

"""
$(TYPEDSIGNATURES)

Get the set of parameters variables for the given system.
"""
function parameters end

include("bipartite_graph.jl")
using .BipartiteGraphs

include("variables.jl")
include("parameters.jl")

include("utils.jl")
include("domains.jl")

include("systems/abstractsystem.jl")

include("systems/diffeqs/odesystem.jl")
include("systems/diffeqs/sdesystem.jl")
include("systems/diffeqs/abstractodesystem.jl")
include("systems/diffeqs/first_order_transform.jl")
include("systems/diffeqs/modelingtoolkitize.jl")
include("systems/diffeqs/basic_transformations.jl")

include("systems/jumps/jumpsystem.jl")

include("systems/nonlinear/nonlinearsystem.jl")

include("systems/optimization/optimizationsystem.jl")

include("systems/control/controlsystem.jl")

include("systems/pde/pdesystem.jl")

include("systems/discrete_system/discrete_system.jl")
include("systems/validation.jl")
include("systems/dependency_graphs.jl")
include("systems/systemstructure.jl")
using .SystemStructures

include("systems/alias_elimination.jl")
include("structural_transformation/StructuralTransformations.jl")
@reexport using .StructuralTransformations

for S in subtypes(ModelingToolkit.AbstractSystem)
    S = nameof(S)
    @eval convert_system(::Type{<:$S}, sys::$S) = sys
end

struct Flow end

export AbstractTimeDependentSystem, AbstractTimeIndependentSystem, AbstractMultivariateSystem
export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system
export DAEFunctionExpr, DAEProblemExpr
export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
export SystemStructure
export JumpSystem
export ODEProblem, SDEProblem
export NonlinearProblem, NonlinearProblemExpr
export OptimizationProblem, OptimizationProblemExpr
export AutoModelingToolkit
export SteadyStateProblem, SteadyStateProblemExpr
export JumpProblem, DiscreteProblem
export NonlinearSystem, OptimizationSystem
export ControlSystem
export alias_elimination, flatten, connect, @connector
export ode_order_lowering, liouville_transform
export runge_kutta_discretize
export PDESystem
export Differential, expand_derivatives, @derivatives
export Equation, ConstrainedEquation
export Term, Sym
export SymScope, LocalScope, ParentScope, GlobalScope
export independent_variables, independent_variable, states, parameters, equations, controls, observed, structure
export structural_simplify
export DiscreteSystem, DiscreteProblem

export calculate_jacobian, generate_jacobian, generate_function
export calculate_control_jacobian, generate_control_jacobian
export calculate_tgrad, generate_tgrad
export calculate_gradient, generate_gradient
export calculate_factorized_W, generate_factorized_W
export calculate_hessian, generate_hessian
export calculate_massmatrix, generate_diffusion_function
export stochastic_integral_transform
export initialize_system_structure
export generate_difference_cb

export BipartiteGraph, equation_dependencies, variable_dependencies
export eqeq_dependencies, varvar_dependencies
export asgraph, asdigraph

export toexpr, get_variables
export simplify, substitute
export build_function
export modelingtoolkitize
export @variables, @parameters, Flow
export @named, @nonamespace, @namespace, extend, compose

end # module


text:
function convolve_gauss(xs::AA{T,1}, ys::AA{T,1}; new_res::T=1.17e5) where T<:AbstractFloat
    # TODO: input spectrum is not infinitely sharp
    # new res is actually convolution res
    σ = (5.0/7.0) * 5434.5232 / new_res / 2.354
    g(x, n) = (one(T)/(σ * sqrt(2.0 * π))) * exp(-0.5 * ((x - n)/σ)^2)

    # pad x array to deal with edges
    # TODO: make padding smarter
    xstep = step(xs)
    newxs = range(xs[1]-100*xstep, xs[end]+100*xstep, step=xstep)

    # find matching indices
    ind1 = searchsortedfirst(newxs, xs[1])
    ind2 = searchsortedfirst(newxs, xs[end])

    # pad y array
    newys = similar(newxs)
    newys[1:ind1-1] .= 1.0
    newys[ind1:ind2] .= ys
    newys[ind2+1:end] .= 1.0

    # perform convolution
    # TODO: do convolution on different grid than newxs
    conv = similar(newxs)
    for i in eachindex(conv)
        conv[i] = sum(newys .* g.(newxs, newxs[i]))
    end

    # normalize and return sliced convolution
    conv ./= maximum(conv)
    return conv[ind1:ind2]
end

function convolve_gauss(xs::AA{T,1}, ys::AA{T,2}; new_res::T=1.17e5) where T<:AbstractFloat
    conv = zeros(length(xs), size(ys,2))
    for i in 1:size(ys,2)
        conv[:, i] = convolve_gauss(xs, ys[:,i], new_res=new_res)
    end
    return conv
end

function degrade_resolution(wave::AA{T,1}, spec::AA{T,1}; R_new::T=700000.0, snr::T=Inf) where T<:AF
    # get new resolution element and wavelength grid
    dλ = wave[1]/R_new
    λs = range(wave[1], wave[end], step=dλ)

    # interpolate the spectrum onto new wavelength grid
    spl = LinearInterpolation(wave, spec, extrapolation_bc=Flat())
    new_spec = spl(λs)

    # add noise and return
    return λs, add_noise(new_spec, snr)
end

function degrade_resolution(wave::AA{T,1}, spec::AA{T,2}; R_new::T=700000.0) where T<:AF
    return [degrade_resolution(wave, spec[:,i]) for i in 1:size(spec,2)]
    # return map(x -> degrade_resolution(wave, spec[:,i]), [])
end

ext: jl

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
include("GO_read.jl");
include("GO_read2.jl");

mutable struct PowersenseData
    nbus::Int  # Number of buses
    ngen::Int  # Number of generators
    nbr::Int  # Number of branches
    nbss::Int  # Number of controllable shunt
    Vmax::Vector{Float64}
    Vmin::Vector{Float64}
    Pmax::Vector{Float64}
    Pmin::Vector{Float64}
    Qmax::Vector{Float64}
    Qmin::Vector{Float64}
    bmax::Vector{Float64}
    bmin::Vector{Float64}
    Imax::Vector{Float64}

    Pg::Vector{Float64}
    Qg::Vector{Float64}
    Bss::Vector{Float64}
    V::Vector{Float64}
    θ::Vector{Float64}
    Gf::Vector{Float64}
    Bf::Vector{Float64}
    Gt::Vector{Float64}
    Bt::Vector{Float64}
    gf::Vector{Float64}
    bf::Vector{Float64}
    gt::Vector{Float64}
    bt::Vector{Float64}

    vr::Vector{Float64}
    vi::Vector{Float64}
    br::Vector{Tuple{Int64,Int64}}
    node::Dict{String,Any}
    Pd::Vector{Float64}
    Qd::Vector{Float64}
    Gl::Vector{Float64}
    Bl::Vector{Float64}

    Iiij::Vector{Float64}
    Iiji::Vector{Float64}
    Irij::Vector{Float64}
    Irji::Vector{Float64}
    Pij::Vector{Float64}
    Qij::Vector{Float64}
    Pji::Vector{Float64}
    Qji::Vector{Float64}

    P::Vector{Float64}
    Q::Vector{Float64}
    lps::Vector{Int}
    pgh::Vector{Vector{Float64}}
    cgh::Vector{Vector{Float64}}

    cost::Float64
    c0::Vector{Float64}
    c1::Vector{Float64}
    c2::Vector{Float64}
    cost_order::Int
    source_type::String

    cg::Vector{Float64}
    slack::Int
    Ag::SparseMatrixCSC{Float64,Int64}
    Abss::SparseMatrixCSC{Float64,Int64}
    Aij::SparseMatrixCSC{Float64,Int64}
    Aji::SparseMatrixCSC{Float64,Int64}

    model

    PowersenseData(nbus::Int,ngen::Int,nbr::Int,nbss::Int) = new(
        nbus,ngen,nbr,nbss,zeros(nbus),zeros(nbus),zeros(ngen),zeros(ngen),zeros(ngen),zeros(ngen),zeros(nbss),zeros(nbss),zeros(nbr),
        zeros(ngen),zeros(ngen),zeros(nbss),ones(nbus),zeros(nbus), zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),
        ones(nbus),zeros(nbus),Array{Tuple{Int64,Int64},1}(),Dict(),zeros(nbus),zeros(nbus),zeros(nbus),zeros(nbus),
        zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),zeros(nbr),
        zeros(nbus),zeros(nbus),zeros(ngen),Vector{Vector{Float64}}(),Vector{Vector{Float64}}(),
        0.0,zeros(ngen),zeros(ngen),zeros(ngen),0, " ",
        zeros(ngen),1,spzeros(nbus, ngen),spzeros(nbus, nbss),spzeros(nbus, nbr),spzeros(nbus, nbr)
    )
end

function create_PowersenseData(path)
    network = process_data(path)
    m = PowersenseData(length(network["bus"]),length(network["gen"]),length(network["branch"]),length(network["sshunt"])-1);
    cost_segments = 30;
    "Assigning bss parameters"
    @simd for i=1:m.nbss sh=network["sshunt"][string(i)]; st=sh["status"]; bus=sh["shunt_bus_index"];
        @inbounds m.bmax[i] = st * sh["shmax"];
        @inbounds m.bmin[i] = st * sh["shmin"];
        m.Abss[bus,i] = 1.0;
    end
    m.slack = network["slack"];
    m.source_type=network["source_type"];

    m.cost_order = (m.source_type == "matpower") ? network["gen_cost_order"] : 0;
    "Assigning generator parameters"
    #@show m.ngen
    @simd for i=1:m.ngen gen=network["gen"][string(i)]; st=gen["gen_status"]; bus=gen["gen_bus_index"];
    	m.Ag[bus,i] = 1.0;
        m.Pmin[i] = st * gen["pmin"];
        m.Pmax[i] = st * gen["pmax"];
        m.Qmin[i] = st * gen["qmin"];
        m.Qmax[i] = st * gen["qmax"];

        if network["source_type"] == "matpower"
        	gen = merge(gen, Dict("LinCost"=>Dict()))
        	if length(network["gen"][string(i)]["cost"]) == 3
        		m.c2[i] = gen["cost"][1];
        		m.c1[i] = gen["cost"][2];
        		m.c0[i] = gen["cost"][3];
        		gen["LinPoints"] = cost_segments;
        		l_seg = (m.Pmax[i] - m.Pmin[i])/ (gen["LinPoints"] - 1)
        		for l=1:gen["LinPoints"]
        			p = m.Pmin[i] + (l-1) * l_seg
        			gen["LinCost"] = merge(gen["LinCost"], Dict(string(l)=>[p, m.c2[i] * p^2 + m.c1[i] * p + m.c0[i]]))
        		end
        	elseif length(network["gen"][string(i)]["cost"]) == 2
        		m.c1[i] = gen["cost"][1];
        		m.c0[i] = gen["cost"][2];
        		gen["LinPoints"] = 2;
        		gen["LinCost"] = merge(gen["LinCost"], Dict(string(1)=>[m.Pmin[i],m.c0[i]]))
        		gen["LinCost"] = merge(gen["LinCost"], Dict(string(2)=>[m.Pmax[i],m.c1[i]]))
        	elseif length(network["gen"][string(i)]["cost"]) == 1
        		m.c0[i] = gen["cost"][1];
        		gen["LinPoints"] = 1;
        		gen["LinCost"] = merge(gen["LinCost"], Dict(string(1)=>[m.Pmin[i],m.c0[i]]))
        	elseif length(network["gen"][string(i)]["cost"]) < 1
        		gen["LinPoints"] = 1;
        		gen["LinCost"] = merge(gen["LinCost"], Dict(string(1)=>[m.Pmin[i],m.c0[i]]))
        	end
        end

        push!(m.cgh, Float64[]); push!(m.pgh, Float64[]);
        #@show i
        m.lps[i] = gen["LinPoints"];
        genCost = gen["LinCost"];
        if (m.Pmin[i] <= m.Pmax[i] && !(m.Pmax[i] == 0 && m.Pmin[i] == 0))
            for l=1:m.lps[i]
                if genCost[string(l)][2] < 10e-6
                    genCost[string(l)][2] = 10e-6
                end
                push!(m.cgh[i],genCost[string(l)][2]);
                push!(m.pgh[i],genCost[string(l)][1]);
            end
            if (m.Pmin[i] < m.pgh[i][1])
                m.pgh[i][1] = m.Pmin[i];
            end
            if (m.Pmax[i] > m.pgh[i][m.lps[i]])
                m.pgh[i][m.lps[i]] = m.Pmax[i];
            end
        elseif (m.Pmax[i] == 0 && m.Pmin[i] == 0)
            m.lps[i] = 1;
            push!(m.cgh[i],0);
            push!(m.pgh[i],0);
        end
    end
     "Assigning bus parameters"
    for key in keys(network["bus"]) bus=network["bus"][key]; i = bus["index"];
        m.Vmax[i] = bus["vmax"];
        m.Vmin[i] = bus["vmin"];
        m.Pd[i] = bus["Pd"];
        m.Qd[i] = bus["Qd"];
        m.Gl[i] = bus["GL"];
        m.Bl[i] = bus["BL"];
    end
    "Assigning branch parameters"
    for i=1:m.nbr br=network["branch"][string(i)]; st=br["br_status"]; f = br["f_bus_index"]; t = br["t_bus_index"];
        g = real(1 / (br["br_r"] + br["br_x"]im));
        b = imag(1 / (br["br_r"] + br["br_x"]im))
        τ = 1/br["tap"] * cos(-br["shift"]) + (1/br["tap"] * sin(-br["shift"]))im
        α = real(τ);
        β = imag(τ);
        m.Gf[i] = st * (-g*α - b*β);
        m.Bf[i] = st * (g*β - b*α);
        m.Gt[i] = st * (-g*α + b*β);
        m.Bt[i] = st * (-g*β - b*α);
        m.gf[i] = st * (br["gfM"] + (g + br["g_fr"]) / (br["tap"] * br["tap"]));
        m.bf[i] = st * (br["bfM"] + (b + br["b_fr"]) / (br["tap"] * br["tap"]));
        m.gt[i] = st * (g + br["g_fr"]);
        m.bt[i] = st * (b + br["b_fr"]);
        m.Imax[i] = haskey(br,"rate_a") ? st * br["rate_a"] : 10000 ;
        m.Imax[i] = m.Imax[i] == 0 ? 10000 : m.Imax[i] ;
        m.Aij[f,i] = 1.0;
        m.Aji[t,i] = 1.0;
        push!(m.br, (f,t));
    end
    m.node = network["nbal"];
    return m
end

function display_preprocess_info(m::PowersenseData)
	println("---> System Information");
	println(" Buses#: ", m.nbus);
	println(" Generators#: ", m.ngen);
	println(" Branches#: ", m.nbr);
	println(" Controllable Shunts#: ", m.nbss);
	println(" Threads: ", Threads.nthreads());
	#println(" ... time_stamp: ", time()-m.start_time, " seconds");
	println("\n");
end


text:
struct DPVectorialTerm{A,K}
    alpha::A
    gamma::K
end

struct DPScalarTerm{A,K}
    alpha::A
    gamma::K
end

const DPTerm = Union{DPVectorialTerm, DPScalarTerm}

quaddata(pot::DPTerm, local_space, charts) = quadpoints(local_space, charts, (2,))
quadrule(pot::DPTerm, local_space, test_index, test_point, trial_index, trial_chart, quad_data) = quad_data[1,trial_index]

function kernelvals(pot::DPTerm, test_point, trial_neighborhood)

    γ = pot.gamma
    α = pot.alpha
    r = test_point - cartesian(trial_neighborhood)
    R = norm(r)

    γR = γ*R
    expn = exp(-γR)
    green = expn / (4π*R)
    gradgreen = -(γ + 1/R) * green / R * r

    return (vector=r, distance=R, green=green, gradgreen=gradgreen)
end

function integrand(pot::DPVectorialTerm, kernel_values, test_point, trial_values, trial_neighborhood)

    α = pot.alpha
    G = kernel_values.green
    f = trial_values.value

    return α*G*f
end

function integrand(pot::DPScalarTerm, kernel_values, test_point, trial_values, trial_neighborhood)

    α = pot.alpha
    ∇G = kernel_values.gradgreen
    f = trial_values.value

    return α*∇G*f
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of basic streams
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael de P. Soares
#* $Id$
#*---------------------------------------------------------------------
type source2
	source2()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Protected=>true
			)),(NComp)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Valid Phases for Flash Calculation",
				:Valid=>["Vapour-Only", "Liquid-Only","Vapour-Liquid"],
				:Default=>"Vapour-Liquid"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Molar or Mass Composition",
				:Valid=>["Molar", "Mass"],
				:Default=>"Molar"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>0,
				:PosY=>0.5256,
				:Symbol=>"_{out}",
				:Protected=>true
			)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Stream Composition"
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Sum of Stream Composition",
				:Protected=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Stream Molar Flow Rate"
			)),
			flow_mass(Dict{Symbol,Any}(
				:Brief=>"Stream Mass Flow"
			)),
			flow_vol(Dict{Symbol,Any}(
				:Brief=>"Volumetric Flow"
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Stream Temperature"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Temperature in  C",
				:Lower=>-250,
				:Upper=>5000
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Stream Pressure"
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Fraction",
				:Hidden=>true
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar Fraction",
				:Hidden=>true
			)),(NComp)),
			molweight(Dict{Symbol,Any}(
				:Brief=>"Average Mol Weight",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Molar Volume",
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Stream Mass Density",
				:Protected=>true
			)),
			dens_mol(Dict{Symbol,Any}(
				:Brief=>"Stream Molar Density",
				:Protected=>true
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Mass Fraction",
				:Protected=>true
			)),(NComp)),
			[
				:(Outlet.z = Composition/sum(Composition)),
				:(zmass = M*Outlet.z / Mw),
				:(zmass = Composition/sum(Composition)),
				:(Outlet.z*sum(zmass/M) = zmass/M),
				:(Outlet.v = 0),
				:(x = Outlet.z),
				:(y = Outlet.z),
				:(Outlet.h = PP.LiquidEnthalpy(Outlet.T, Outlet.P, x)),
				:(vm = PP.LiquidVolume(Outlet.T, Outlet.P, x)),
				:(Outlet.v = 1),
				:(x = Outlet.z),
				:(y = Outlet.z),
				:(Outlet.h = PP.VapourEnthalpy(Outlet.T, Outlet.P, y)),
				:(vm = PP.VapourVolume(Outlet.T, Outlet.P, y)),
				:([Outlet.v, x, y] = PP.Flash(Outlet.T, Outlet.P, Outlet.z)),
				:(Outlet.h = (1-Outlet.v)*PP.LiquidEnthalpy(Outlet.T, Outlet.P, x) + Outlet.v*PP.VapourEnthalpy(Outlet.T, Outlet.P, y)),
				:(vm = (1-Outlet.v)*PP.LiquidVolume(Outlet.T, Outlet.P, x) + Outlet.v*PP.VapourVolume(Outlet.T,Outlet.P,y)),
				:(SumOfComposition = sum(Composition)),
				:(rhom * vm = 1),
				:(Mw = sum(M*Outlet.z)),
				:(rhom * Mw = rho),
				:(Fw = Mw*Outlet.F),
				:(Fvol = Outlet.F*vm),
				:(T_Cdeg = Outlet.T/"K" - 273.15),
				:(Outlet.F = F),
				:(Outlet.P = P),
				:(Outlet.T = T),
			],
			[
				"Stream Molar Composition","Stream Mass Composition","Stream Mass Composition","Stream Molar Composition","Vapour Fraction","Liquid Composition","Vapour Composition","Overall Enthalpy","Molar Volume","Vapor Fraction","Liquid Composition","Vapour Composition","Overall Enthalpy","Molar Volume","Flash Calculation","Overall Enthalpy","Molar Volume","Sum of Composition","Molar Density","Average Molecular Weight","Mass or Molar Density","Flow Mass","Volumetric Flow","Temperature in  C","Equate Flow","Equate Pressures","Equate Temperatures",
			],
			[:PP,:NComp,:M,:ValidPhases,:CompositionBasis,],
			[:Outlet,:Composition,:SumOfComposition,:F,:Fw,:Fvol,:T,:T_Cdeg,:P,:x,:y,:Mw,:vm,:rho,:rhom,:zmass,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	M::Array{molweight}
	ValidPhases::DanaSwitcher
	CompositionBasis::DanaSwitcher
	Outlet::stream
	Composition::Array{positive}
	SumOfComposition::positive
	F::flow_mol
	Fw::flow_mass
	Fvol::flow_vol
	T::temperature
	T_Cdeg::DanaReal
	P::pressure
	x::Array{fraction}
	y::Array{fraction}
	Mw::molweight
	vm::volume_mol
	rho::dens_mass
	rhom::dens_mol
	zmass::Array{fraction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export source2
function set(in::source2)
	M = PP.MolecularWeight()

end
function setEquationFlow(in::source2)
	let switch=CompositionBasis
		if switch=="Molar"
			addEquation(1)
			addEquation(2)
		elseif switch=="Mass"
			addEquation(3)
			addEquation(4)
		end
	end
	let switch=ValidPhases
		if switch=="Liquid-Only"
			addEquation(5)
			addEquation(6)
			addEquation(7)
			addEquation(8)
			addEquation(9)
		elseif switch=="Vapour-Only"
			addEquation(10)
			addEquation(11)
			addEquation(12)
			addEquation(13)
			addEquation(14)
		elseif switch=="Vapour-Liquid"
			addEquation(15)
			addEquation(16)
			addEquation(17)
		end
	end
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	addEquation(23)
	addEquation(24)
	addEquation(25)
	addEquation(26)
	addEquation(27)
end
function atributes(in::source2,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Source2"
	fields[:Brief]="Material stream source"
	fields[:Info]="
	This model should be used for boundary streams.
	Usually these streams are known and come from another process
	units.

	The user should specify:
	 * Total molar (mass or volumetric) flow
	 * Temperature
	 * Pressure
	 * Molar or mass composition

	No matter the specification set, the model will calculate some
	additional properties:
	 * Mass density
	 * Mass flow
	 * Mass compostions
	 * Specific volume
	 * Vapour fraction
	 * Volumetric flow
	 * Liquid and Vapour compositions
	"
	drive!(fields,_)
	return fields
end
source2(_::Dict{Symbol,Any})=begin
	newModel=source2()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
immutable DataArray{T<:Number,D<:Number,A<:AbstractArray} <: AbstractData
    index::Vector
    values::A
    DataArray(i::AbstractVector{T},d::AbstractArray{D}) = new(collect(i),d)
end

function _data(::Type{Val{:array}},i::AbstractVector,d::AbstractArray)
    @assert length(i) == size(d,1)
    DataArray{eltype(i),eltype(d),typeof(d)}(i,d)
end

numvalues(d::DataArray) = size(d.values,1)
numvars(d::DataArray) = size(d.values,2)
eltype(d::DataArray) = eltype(d.values)
generate!(d::DataArray) = d.values
dataindex(d::DataArray) = d.index
datavalues(d::DataArray) = d.values
datatypename(d::DataArray) = "DataArray"


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.ContactDetectionMPR (Modia3D/contactDetection/ContactDetectionMPR/_module.jl)
#

const Point = SVector{3,Float64}

const i16max = Int64(typemax(Int16))
const i32max = Int64(typemax(Int32))

pack16(i1::Integer, i2::Integer) = Int64(i1) + i16max*Int64(i2)

function pack(i1::Integer, i2::Integer,i3::Integer,i4::Integer)
   @assert(i1 >= 0 && i1 <= typemax(Int16))
   @assert(i2 >= 0 && i2 <= typemax(Int16))
   @assert(i3 >= 0 && i3 <= typemax(Int16))
   @assert(i4 >= 0 && i4 <= typemax(Int16))
   return pack16(i1,i2) + i32max*pack16(i3,i4)
end

function unpack16(i::Int64)
    i1 = rem(i,i16max)
    i2 = div(i-i1, i16max)
    return (i1,i2)
end

function unpack32(i::Int64)
    i1 = rem(i,i32max)
    i2 = div(i-i1, i32max)
    return (i1,i2)
end

function unpack(i::Int64)
    @assert(i >= 0 && i <= typemax(Int64))
    (tmp1,tmp2) = unpack32(i)
    (i1,i2) = unpack16(tmp1)
    (i3,i4) = unpack16(tmp2)
    return (i1,i2,i3,i4)
end

### -------------------computation of pairID -----------------------------------
### it returns a unique ID
function orderPositions(is,i,js,j)
    if is < js
        return pack(is,i,js,j)
    elseif is > js
        return pack(js,j,is,i)
    else
        error("from orderPositions: is == js.")
end; end

getPositionsOfObj(scene::Composition.Scene, obj::Composition.Object3D,
                  movablePos::Nothing) = (false, 0, 0)
getPositionsOfObj(scene::Composition.Scene, obj::Composition.Object3D,
                  movablePos::Int64) = (true, movablePos, findall(x->x==obj, scene.superObjs[movablePos].superObjMovable.superObj)[1] )

function computePairID(scene::Composition.Scene,
        actObj::Composition.Object3D, nextObj::Composition.Object3D,
        is, i, js, j)
    # is: actual super - object
    # js: subsequent super - object
    # i: Object3D of is_th super - object
    # j: Object3D of js_th super - object
    (isSetAct, isPosActSuper, iPosActObj)    = getPositionsOfObj(scene, actObj,  actObj.interactionManner.movablePos)
    (isSetNext, jsPosNextSuper, jPosNextObj) = getPositionsOfObj(scene, nextObj, nextObj.interactionManner.movablePos)
    if isSetAct && isSetNext
        return orderPositions(isPosActSuper,iPosActObj,jsPosNextSuper,jPosNextObj)
    elseif isSetAct && !isSetNext
        return orderPositions(isPosActSuper,iPosActObj,js,j)
    elseif isSetNext && !isSetAct
        return orderPositions(is,i,jsPosNextSuper,jPosNextObj)
    else
        return orderPositions(is,i,js,j)
end; end


text:
struct Color
    c::Union{Symbol, Int}
end
Color() = Color(:normal)

struct Token
    s::String
    c::Color
end

function parsecolor(s, i)
    i >= ncodeunits(s) && return i, Color(), false
    i = nextind(s, i)
    i1 = i
    while i1 <= ncodeunits(s)
        c = s[i1]
        isspace(c) && return i, Color(), false
        if c == '{'
            i1 == i && return i, Color(), false
            val = s[i:prevind(s, i1)]
            v1 = tryparse(Int, val)
            i1 = nextind(s, i1)
            if v1 !== nothing
                return i1, Color(v1), true
            else
                return i1, Color(Symbol(val)), true
            end
        end
        i1 = nextind(s, i1)
    end

    return i, Color(), false
end

function parse!(s, i, res, col, depth)
    isescape = false
    i1 = i
    i2 = i
    while i2 <= ncodeunits(s)
        c = s[i2]
        if c == '\\'
            isescape = !isescape
        elseif c == ':' && !isescape
            i0, col2, iscol = parsecolor(s, i2)
            if iscol
                i1 < i2 && push!(res, Token(s[i1:prevind(s, i2)], col))
                i2 = parse!(s, i0, res, col2, 1)
                i1 = nextind(s, i2)
            end
        elseif c == '{' && !isescape
            depth += (depth > 0) # we do not calculate open and closed brackets for initial parse
        elseif c == '}' && !isescape
            depth -= (depth > 0)
            if depth == 0
                push!(res, Token(s[i1:prevind(s, i2)], col))
                return i2
            end
        else
            isescape = false
        end
        i2 = nextind(s, i2)
    end

    i1 < i2 && push!(res, Token(s[i1:prevind(s, i2)], col))
    return i2
end

function parse(s)
    res = Token[]
    col = Color(:normal)
    parse!(s, 1, res, col, -1)

    return res
end

macro c_str(s)
    quote
        s2 = @interpolate $s
        parse(s2)
    end
end


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf, JLD2, FileIO

# noise1=convert(Array{BigFloat},(rand(n).-1/2)./BigFloat("1e1000"))
# f=sin.(25x)
# f=f+noise1

p=50;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("10"); n=BigInt(100); iN=convert(Int64,n);
dx=BigFloat("2")*L/(n-1); x=convert(Array{Float64},-L:dx:L)
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end

V=TaylorBigF.calc_Vandermonde(x,n,n)
xi=(x.+dx/BigFloat("2"))[1:end-1]
Vi=TaylorBigF.calc_Vandermonde(xi,n-1,n)

# f=convert(Array{BigFloat},(rand(n).-1/2))
# tansf1=BigFloat("1e100")
# f_transform=exp.(-f.^2)
f=sin.(x)+sin.(x./BigFloat("2"))

# fi=(f[2:end]+f[1:end-1])/2
fi=sin.(xi)+sin.(xi./BigFloat("2"))
# fi_transform=exp.(-fi.^2)

# f[75]+=BigFloat("1e-20")
f+=convert(Array{BigFloat},(rand(n).-1/2)./1)
# plot(x,f)

a=V\f
# a.*=(1 .+(rand(n).-1/2)/1e3)
r_root=Array{Float64}(undef,0); for i=1:Int64(n-1) push!(r_root,abs(a[i])^(1/(n-1))) end; r_rootF=convert(Array{Float64},r_root);
# plot(r_rootF)
# plot(sign.(a).*log10.(abs.(a)))
error1=maximum(abs.(f-V*a))
fi_comp=Vi*a
error2=maximum(abs.(fi_comp-fi))

plot(x,f)
scatter!(x,V*a,markersize=1)
scatter!(xi,fi,markersize=1)
scatter!(xi,fi_comp,markersize=1)
# fi_comp=sqrt.(log.(-fi_comp_transform))
# fi_comp=fi_comp_transform*tansf1
# fi_analyt=sin.(25xi);

# dV=TaylorBigF.calc_Vandermonde_dot(x,n,n)

# Eq1=V[1,:]
# q1=f[1]
# Eq2=dV[1,:]
# q2=(f[2]-f[1])/dx
# Eq3=V[end,:]
# q3=f[end]
# Eq4=dV[end,:]
# q4=(f[end]-f[end-1])/dx

# VV=[V;Eq1';Eq2';Eq3';Eq4']
# q_all=[f;q1;q2;q3;q4]
# a=VV[1:end,:]\q_all[1:end]
# r_root=Array{Float64}(undef,0); for i=1:Int64(n-1) push!(r_root,abs(a[i])^(1/(n-1))) end; r_rootF=convert(Array{Float64},r_root);
# plot(r_rootF)
# error1=maximum(abs.(f-V*a))
# plot(x,f)
# scatter!(x,V*a)
# Vi=TaylorBigF.calc_Vandermonde(xi,n-1,n);
# fi_comp=Vi*a
# error2=maximum(abs.(fi_comp-fi))
# scatter!(xi,Vi*a)


text:
# coding: utf-8

# In[1]:

#i would recommend you to read algorithm and data structure before messing with hashtable
# http://interactivepython.org/runestone/static/pythonds/SortSearch/Hashing.html
#this script features the hash function described in the book, which is mod eleven %11
#we apply hash function to all the values in a given list
#we shall obtain the output value for each element
#and assuming this is a perfect hash function
#each value from the list should have a unique hash value
#we can create a dictionary based on the unique hash value
#and assign all the values from the list into dictionary
#when we wanna search for something
#we just need to apply hash function to the target
#and make query on the dictionary for that particular hash value
#O(1), simple AF, not really...
#we could encounter the same hash value for different values from the list
#this is so called hash collision
#there are three ways of solving hash collision in the book
#which are chaining, linear probing and quadratic probing
#lets deliberately create an imperfect hash function to see how it goes

# In[2]:

arr=[21,55,89,67,12,12]

# In[3]:

#hash chaining
#the easiest one
#basically stack all the collision together
#form a list under that hash value
#does it make search easier?
#not really, if the list under one hash value gets too large
#it would slow down the search
#the good thing about this method is that you always have all the values in the dictionary
#for other methods, you either have to increase hash table size or drop values

# In[4]:

#the first function is to create a dictionary
#assign values from the list to the dictionary based on hash value
function create_hashtable(arr,default=11)

    #create hashtable
    hashtable=Dict()

    #initialize
    for i in 0:(default-1)
        hashtable[i]=[]
    end

    #add value into hashtable
    #using folding method
    for i in arr
        push!(hashtable[i%default],i)
    end

    return hashtable
end

# In[5]:

#now its the search part
#we just apply hash function on target and get hash value
#we look up the hash value in dictionary
function hash_search(target,arr,default=11)

    #get hashtable
    hashtable=create_hashtable(arr,default)

    #find value through hash value
    if target in hashtable[target%default]
        return true
    else
        return false
    end
end

# In[6]:

hash_search(12,arr)

# In[7]:

#linear probing
#when collision occurs, we try to find the next empty slot to store it
#it sounds reasonable but it is also trouble
#what if we run through the whole dictionary
#and there is no slot?
#we can choose to drop the values
#or we can reset the hash function or expand the dictionary
#in the best case, it is faster than chaining
#in the worst case, it is slower

# In[8]:

#the first function is to create a dictionary
#assign values from the list to the dictionary based on hash value
function create_hashtable(arr,default=11)

    #create hashtable
    hashtable=Dict()

    #unwanted is a temporary list to append collision items
    unwanted=[]

    #keep track of values with no allocation
    badhash=[]

    #initialize
    for i in 0:(default-1)
        hashtable[i]=""
    end

    #add value into hashtable
    for i in arr
        if hashtable[i%default]==""
            hashtable[i%default]=i
        else
            #make sure every collision will be tracked
            push!(unwanted,i)
        end
    end

    for i in unwanted

        #c is a counter
        #in the inner loop
        #c is to determine whether we have gone through the entire list
        c=0
        modnumber=i%default
        stop=false

        while !stop && c<=default-1

            #when modnumber exceeds ten, we return it to 0
            #alternatively we can use mod eleven %11
            modnumber+=1
            if modnumber>default-1
                modnumber=0
            end

            #when the next slot isnt empty, we keep iterating
            if hashtable[modnumber]==""
                hashtable[modnumber]=i
                stop=true
            end

            c+=1
        end

        #make sure that we will print out those items which didnt get assigned
        if !stop
            push!(badhash,i)
        end
    end

    #if the hashing is imperfect, we print out badhash list
    if !isempty(badhash)
        println(badhash)
    end

    return hashtable
end

# In[9]:

#the search part is very similar to the chaining one
function hash_search(target,arr,default=11)

    #get hashtable
    hashtable=create_hashtable(arr,default)

    #find value through hash value
    modnumber=target%default
    if target in hashtable[modnumber]
        return true
    else

        #when we cannot find the value at hash value
        #we begin our linear probing
        #its the same process as the hash function
        #except we only need to return T/F
        counter=0

        while counter<=default-1
            modnumber+=1

            if modnumber>default-1
                modnumber=0
            end

            if target==hashtable[modnumber]
                return true
            end
            counter+=1
        end
    end
    return false
end

# In[10]:

hash_search(89,arr)

# In[11]:

#quadratic probing
#it sounds math intensive with the word quadratic
#as a matter of fact, it is simple AF
#we just replace the add one method with add quadratic values

# In[12]:

#the first function is to create a dictionary
#assign values from the list to the dictionary based on hash value
function create_hashtable(arr,default=11)

    #create hashtable
    hashtable=Dict()

    unwanted=[]
    badhash=[]

    #initialize
    for i in 0:(default-1)
        hashtable[i]=""
    end

    #add value into hashtable
    for i in arr
        if hashtable[i%default]==""
            hashtable[i%default]=i
        else
            push!(unwanted,i)
        end
    end

    for i in unwanted
        c=1
        modnumber=i%default
        stop=false

        #the loop is basically the same as linear probing
        #except we add quadratic value
        #note that its quite difficult
        #to determine whether we have been through the entire list
        #so i still set counter limit at 10
        while !stop && c<=default

            modnumber+=c^2

            #note that i use mod eleven %11 when iteration exceeds hash table size
            if modnumber>default-1
                modnumber=modnumber%default
            end

            if hashtable[modnumber]==""
                hashtable[modnumber]=i
                stop=true
            end

            c+=1
        end

        if !stop
            push!(badhash,i)
        end
    end

    if !isempty(badhash)
        println(badhash)
    end

    return hashtable
end

# In[13]:

#the search is basically the same as linear probing
#except linear part is substituted with quadratic
function hash_search(target,arr,default=11)

    #get hashtable
    hashtable=create_hashtable(arr,default)

    #find value through hash value
    modnumber=target%default
    if target in hashtable[modnumber]
        return true
    else
        counter=1

        while counter<=default
            modnumber+=counter^2

            if modnumber>default-1
                modnumber=modnumber%default
            end

            if target==hashtable[modnumber]
                return true
            end
            counter+=1
        end
    end
    return false
end

# In[14]:

hash_search(78,[21,55,89,67,12,12,12,12,12,12,12,12,12,12,78])

# In[15]:

#we get False in the end
#its quite interesting that for the same hash value 67,12,78
#we can store 67 in hash table but not 78
#67 and 12 are processed earlier than 78
#quadratic probing doesnt iterate through all slots
#all empty slots we can iterate have been occupied by the time we reach 78

max_stars_repo_path: hash search.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: hash search.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: hash search.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
## This code implements algorithms in
## Rudolfo Gorenflo, Joulia Loutchko and Yuri Loutchko,
## *Computation of the Mittag-Leffler function and its derivative*,  Fract. Calc. Appl. Anal, **(2002)**

# When support for v0.6 is dropped, use fully qualified name
if VERSION >= v"0.7-"
    import SpecialFunctions: gamma, erfc
end

macro br(n)
#    esc(:(println(stderr, "branch ", $n)))
    nothing
end

function P(α,β,ϵ,ϕ,z)
    ω = ϕ * (1+(1-β)/α) + ϵ^(1/α) * sin(ϕ/α)
    return (1/(2*α*pi)) * ϵ^(1+(1-β)/α)*exp(ϵ^(1/α) * cos(ϕ/α)) * (cos(ω) + im * sin(ω))/(ϵ*exp(im*ϕ)-z)
end

if VERSION >= v"0.5-"
    import QuadGK: quadgk
end

ourquadgk(f,a,b) = quadgk(f,a,b; order=7)[1]

function Kint(α,β,a,χ0,z)
    return ourquadgk(χ -> K(α,β,χ,z), a, χ0)
end

function K(α,β,χ,z)
    den = (χ^2-2*χ*z*cos(α*pi)+z^2)
    return (1/(α*pi)) * χ^((1-β)/α) * exp(-χ^(1/α))*( χ * sin(pi*(1-β)) - z * sin(pi*(1-β+α)))/den
end

Kint(α,β,χ0,z) = Kint(α,β,0,χ0,z)

mpow(x::Complex, y) = x^y
mpow(x::Real, y) = x >= 0 ? x^y : Complex(x,0)^y

# The following gymnastics are to get around the type-unstable mpow
function sum2(α, β, z::Real, k0)
    return z > 0 ? sum2_pos(α, β, z, k0) : sum2_neg(α, β, z, k0)
end
sum2(α, β, z::Complex, k0) = sum2_neg(α, β, z, k0)

for funcname in (:sum2_neg, :sum2_pos)
    local zarg
    local stype
    if funcname == :sum2_neg
        zarg = :( (z + Complex(0,0) ))
        stype = :( typeof(z + Complex(0,0) ) )
    else
        zarg = :z
        stype = :( typeof(z) )
    end
    @eval function ($funcname)(α, β, z, k0)
        s::($stype) = zero($stype)
        for k=1:k0
            arg = β - α * k
            if !( round(arg) == arg && arg < 0)
                s += ($zarg)^(-k) / gamma(arg)
            end
        end
        return s
    end
end

function choosesum(α,β,z,ρ)
    k0 = floor(Int, -log(ρ)/log(abs(z)))
    if abs(angle(z)) < pi*α/4 + 1//2 * min(pi,pi*α)
        @br 3
        return 1/α * z^((1-β)/α) * exp(z^(1/α)) - sum2(α,β,z,k0)
    else
        @br 4
        return - sum2(α,β,z,k0)
    end
end

function mittleffsum(α, β, z)
    @br 1
    k0::Int = floor(Int, α) + 1
    s = zero(z)
    for k=0:(k0 - 1)
        s += mittleff(α / k0, β, mpow(z, (1 // k0)) * exp(2pi * im * k / k0))
    end
    return s / k0
end

# mittleffsum sometimes passes complex z with magnitude nearly equal to 1
# Then k0 is astronomically high. Because the gamma function is evaluated by a call to gsl,
# the loop is not interruptible.
# So, we compute k0 using only the real part of k0.
# The effect on the accuracy is not controlled.
function mittleffsum2(α,β,z,ρ)
    @br 2
    zr = real(z)
    k0 = max(ceil(Int,(1-β)/α), ceil(Int, log(ρ*(1-abs(zr)))/log(abs(zr))))
    s = zero(z)
    for k=0:k0
        s += z^k/gamma(β+α*k)
    end
    return s
end

Pint(α,β,ϵ,z) = ourquadgk(ϕ -> P(α,β,ϵ,ϕ,z), -α*pi, α*pi)
Pint(α,β,z) = Pint(α,β,1,z)

function mittleffints(α,β,z,ρ)
    az = abs(z)
    ab = abs(β)
    χ0 = β >= 0 ?
      max(1,2*az,(-log(pi*ρ/6))^α) :
      max((ab+1)^α, 2*az,(-2*log( pi*ρ/(6*(ab+2)*(2*ab)^ab)))^α)
    aaz = abs(angle(z))
    if aaz > α * pi
        if β <= 1
            @br 5
            return Kint(α,β,χ0,z)
        else
            @br 6
            return Kint(α,β,1,χ0,z) + Pint(α,β,z)
        end
    elseif aaz < pi*α
        if β <= 1
            @br 7
            return Kint(α,β,χ0,z) + (1/α)*z^((1-β)/α) * exp(z^(1/α))
        else
            @br 8
            return Kint(α,β,az/2,χ0,z) + Pint(α,β,(az)/2,z) + (1/α)*z^((1-β)/α) * exp(z^(1/α))
        end
    else
        @br 9
        return Kint(α,β,(az+1)/2,χ0,z) + Pint(α,β,(az+1)/2,z)
    end
end

"""
    mittlefferr(α,z,ρ)

Compute mittlefferr(α,1,z,ρ).
"""
mittlefferr(α,z,ρ) = mittlefferr(α,1,z,ρ)

"""
    mittlefferr(α,β,z,ρ)

Compute the Mittag-Leffler function at `z` for parameters `α,β` with
accuracy `ρ`.
"""
function mittlefferr(α,β,z,ρ::Real)
    ρ > 0 || throw(DomainError(ρ))
    _mittlefferr(α,β,z,ρ)
end

_mittlefferr(α::Real, β::Real, z::Real, ρ::Real) = real(_mittleff(α, β, z, ρ))
_mittlefferr(α::Real, β::Real, z::Complex, ρ::Real) = _mittleff(α, β, z, ρ)

# The second definition would work for both complex and real
myeps(x) = x |> one |> float |> eps
myeps(x::Complex) =  x |> real |> myeps

"""
    mittleff(α,β,z)

Compute the Mittag-Leffler function at `z` for parameters `α,β`.
"""
mittleff(α, β, z) = _mittleff(α, β, float(z))
mittleff(α, β, z::Union{Integer,Complex{T}}) where {T<:Integer} = mittleff(α, β, float(z))

#mittleff(α, β, z) = _mittlefferr(α,β,z,myeps(z))

"""
    mittleff(α,z)

Compute `mittleff(α,1,z)`.
"""
mittleff(α, z) = _mittleff(α,1,z)
#mittleff(α, z::Union{Integer,Complex{T}}) where {T<:Integer} = mittleff(α, float(z))

function _mittleff_special_beta_one(α,z)
    z == 0 && return one(z)
    if α == 1/2 && abs(z) < 10
        return exp(z^2) * erfc(-z)
    end
    α == 0 && return 1/(1-z)  # FIXME needs domain check
    α == 1 && return exp(z)
    zc = isa(z, Real) && z < 0 ? complex(z) : z # Julia sometimes requires explicit complex numbers for efficiency
    α == 2 && return cosh(sqrt(zc))
    α == 3 && return (1//3) * (exp(zc^(1//3)) + 2 * exp(-zc^(1//3)/2) * cos(sqrt(convert(typeof(zc),3))/2 * zc^(1//3)))
    α == 4 && return (1//2) * (cosh(zc^(1//4)) + cos(zc^(1//4)))
    return nothing
end

function _mittleff_slow_with_eps(α, β, z, ρ)
    az = abs(z)
    az < 1 && return mittleffsum2(α,β,z,ρ)
    az > floor(10+5*α) && return choosesum(α,β,z,ρ)
    return mittleffints(α,β,z,ρ)
end

_mittleff(α::Real, β::Real, z::Real) = real(_mittleff0(α, β, z))
_mittleff(α, β, z) = _mittleff0(α, β, z)

function _mittleff0(α, β, z)
    if β == 1
        res = _mittleff_special_beta_one(α,z)
        res != nothing && return res
    end
    z == 0 && return 1/gamma(β)
    if β == 2
        α == 1 && return (exp(z) - 1)/z
        zc = isreal(z) && z < 0 ? complex(z) : z
        α == 2 && return sinh(sqrt(zc))/sqrt(zc)
    end
    1 < α && return mittleffsum(α,β,z)
    ρ = myeps(z)
    return _mittleff_slow_with_eps(α, β, z, ρ)
end

function _mittleff(α,β,z,ρ)
    if β == 1
        res = _mittleff_special_beta_one(α,z)
        res != nothing && return res
    end
    z == 0 && return 1/gamma(β)
    α == 1 && β == 1 && return(exp(z))
    1 < α && return mittleffsum(α,β,z)
    return _mittleff_slow_with_eps(α, β, z, ρ)
end

"""
    mittleffderiv(α,β,z)

Compute the derivative of the Mittag-Leffler function at `z` for parameters `α,β`.
"""
function mittleffderiv(α, β, z)
    #derivative of Mittag Leffler function WRT to main argument Z.
    #take q = 0.5. Paper requires |z| <= q < 1.
    q = 1//2
    #case 1, small z
    if abs(z) <= q
        ω = α + β - 3//2
        D = α^2 - 4*α*β + 6*α + 1
        #k1
        if α>1
            k₁ = ((2-α-β)/(α-1)) + 1
        elseif α>0 && α<=1 && D<=0
            k₁ = ((3-α-β)/α) + 1
        else
            k₁ = maximum([((3-α-β)/α) + 1, ((1-2*ω*α+sqrt(D))/(2*(α^2)))+1])
        end
        k₀ = maximum([k₁, log(myeps(z)*(1-abs(z)))/log(abs(z))])
        k₀ = ceil(Int,k₀) #take ceiling (not specified in paper whether floor or ceiling)
        out = zero(z)
        for k in 0:k₀
            out = out + ((k+1)*z^k)/(gamma(α+β+α*k))
        end
    #case 2, larger z
    else
        out = (mittleff(α,β-1,z) - (β-1)*mittleff(α,β,z))/(α*z)
    end
    return out
end

"""
    mittleffderiv(α,z)

Compute mittleffderiv(α,1,z)
"""
mittleffderiv(α, z) = mittleffderiv(α,1,z)


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

# 12 steps to Navier–Stokes
# =====
# ***

# In the previous step, we solved the [2D Burgers' equation](./10_Step_8.ipynb): an important equation in the study of fluid mechanics because it contains the full convective nonlinearity of the flow equations. With that exercise, we also build the experience to incrementatlly code a Navier–Stokes solver.

# In the next two steps, we will solve Laplace and then Poisson equation. We will then put it all together!

# Step 9: 2D Laplace Equation
# ----
# ***

# Here is Laplace's equation in 2D:

# $$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0$$

# We know how to discretize a 2nd order derivative. But think about this for a minute — Laplace's equation has the features typical of diffusion phenomena. For this reason, it has to be discretized with *central differences*, so that the discretization is consistent with the physics we want to simulate.

# The discretized equation is:

# $$\frac{p_{i+1, j}^n - 2p_{i,j}^n + p_{i-1,j}^n}{\Delta x^2} + \frac{p_{i,j+1}^n - 2p_{i,j}^n + p_{i, j-1}^n}{\Delta y^2} = 0$$

# Notice that the Laplace Equation does not have a time dependence — there is no $p^{n+1}$.  Instead of tracking a wave through time (like in the previous steps), the Laplace equation calculates the equilibrium state of a system under the supplied boundary conditions.

# If you have taken coursework in Heat Transfer, you will recognize the Laplace Equation as the steady-state heat equation.

# Instead of calculating where the system will be at some time $t$, we will iteratively solve for $p_{i,j}^n$ until it meets a condition that we specify.  The system will reach equilibrium only as the number of iterations tends to $\infty$, but we can approximate the equilibrium state by iterating until the change between one iteration and the next is *very* small.

# Let's rearrange the discretized equation, solving for $p_{i,j}^n$:

# $$p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}$$

# Using second-order central-difference schemes in both directions is the most widely applied method for the Laplace operator. It is also known as the **five-point difference operator**, alluding to its stencil.

# We are going to solve Laplace's equation numerically by assuming an initial state of $p=0$ everywhere. Then we add boundary conditions as follows:

# $p=0$ at $x=0$

# $p=y$ at $x=2$

# $\frac{\partial p}{\partial y}=0$ at $y=0, \ 1$

# Under these conditions, there is an analytical solution for Laplace's equation:

# $$p(x,y)=\frac{x}{4}-4\sum_{n=1,odd}^{\infty}\frac{1}{(n\pi)^2\sinh2n\pi}\sinh n\pi x\cos n\pi y$$

# ### Exercise

# Write your own code to solve Poisson's equation using loops, in the style of coding used in our first lessons.
# Then, consider the demonstration of how to write it using functions (below) and modify your code in that style.
# Can you think of reasons to abandon the old style and adopt modular coding?

# Other tips:

# + Visualize each step of the iterative process
# + Think about what the boundary conditions are doing
# + Think about what the PDE is doing

# ### Using functions

# Remember the lesson on writing [functions with Julia](./11_Defining_Functions_in_Julia.jl)? We will use that style of code in this exercise.

# We're going to define two functions: one that plots our data in a 3D projection plot and the other that iterates to solve
# for $p$ until the change in the [L1 Norm][1] of $p$ is less than a specified value.

# [1]: http://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm

using Plots

# TODO PLOT
function  plot2D(x, y, p)
    fig = pyplot.figure(figsize=(11, 7), dpi=100)
    ax = fig.gca(projection='3d')
    X, Y = numpy.meshgrid(x, y)
    surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,
            linewidth=0, antialiased=False)
    ax.set_xlim(0, 2)
    ax.set_ylim(0, 1)
    ax.view_init(30, 225)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
end

# The function `plot2D` takes three arguments, an x-vector, a y-vector and our p matrix.
#     Given these three values, it produces a 3D projection plot,
#     sets the plot limits and gives us a nice viewing angle.

# $$p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}$$

function laplace2d(p, y, dx, dy, l1norm_target)
    l1norm = 1
    pn = empty(p)

    while l1norm > l1norm_target
        pn = copy(p)
        p[2:end, 2:end] = ((dy^2 * (pn[2:end, 2:] + pn[2:end, 0:-2]) +
                         dx^2 * (pn[2:, 2:end] + pn[0:-2, 2:end])) /
                        (2 * (dx^2 + dy^2)))

        p[:, begin] = 0  # p = 0 @ x = 0
        p[:,   end] = y  # p = y @ x = 2
        p[begin, :] = p[2, :]  # dp/dy = 0 @ y = 0
        p[end,   :] = p[end-1, :]  # dp/dy = 0 @ y = 1
        # TODO Port
        l1norm = (sum(abs.(p[:]) - abs.(pn[:])) /
                sum(abs.(pn[:]))
    end

    return p
end

# `laplace2d` takes five arguments, the `p` matrix, the `y`-vector, `dx`, `dy` and the value `l1norm_target`.
# This last value defines how close the `p` matrix should be in two consecutive iterations before the loop breaks and returns the calculated `p` value.

# Note that when executing the cells above in your own notebook, there will be no output.
# You have *defined* the function but you have not yet *called* the function.
# It is now available for you to use, the same as `numpy.linspace` or any other function in our namespace.

## Variable declarations
nx = 31
ny = 31
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)

# initial conditions
p = zeros((ny, nx))  # create a XxY vector of 0's

# plotting aids
x = LinRange(0, 2, nx)
y = LinRange(0, 1, ny)

##boundary conditions
p[:, begin] = 0  # p = 0 @ x = 0
p[:, end]   = y  # p = y @ x = 2
p[begin, :] = p[2, :]  # dp/dy = 0 @ y = 0
p[end, :]   = p[end-1, :]  # dp/dy = 0 @ y = 1

# Now let's try using our `plot2D` function to look at our initial conditions.
# If the function has been correctly defined,
# you should be able to begin typing `plot2D`
# and hit the **tab** key for auto-complete options.

# TODO PLOT
plot2D(x, y, p)

# It worked!  This is the initial state of our problem, where the value of `p`
# is zero everywhere except for along $x=2$ where $p=y$.
# Now let's try to run our `laplace2d` function with a specified L1 target of .01

# [Hint: if you are having trouble remembering the order in which variables are sent to a function,
# you can just type `laplace2d(` and the iPython Notebook will put up a little popup box to remind you]

p = laplace2d(p, y, dx, dy, 1e-4)

# Now try plotting this new value of `p` with our plot function.

# TODO PLOT
plot2D(x, y, p)

# ***

## Learn More

# The [next step](./13_Step_10.ipynb) will be to solve Poisson's equation. Watch **Video Lesson 11** on You Tube to understand why we need Poisson's equation in CFD.
#src="https://www.youtube.com/embed/ZjfxA3qq2Lg"

# And for a detailed walk-through of the discretization of Laplace and Poisson equations (steps 9 and 10), watch **Video Lesson 12** on You Tube:
#src="https://www.youtube.com/embed/iwL8ashXhWU"

max_stars_repo_name: miguelraz/CFDPython

max_issues_repo_name: miguelraz/CFDJulia

max_forks_repo_name: miguelraz/CFDJulia


text:
function read_analog(t::NIDAQ.AITask, precision::DataType, num_samples_per_chan::Integer = -1; timeout::Float64=-1.0)
    data = allocate_sampbuf(t, precision, num_samples_per_chan)
    read_analog!(data, t, precision, num_samples_per_chan; timeout = timeout)
end

function allocate_sampbuf(t::NIDAQ.AITask, precision::DataType, num_samples_per_chan::Integer = -1)
    num_channels = getproperties(t)["NumChans"][1]
    buffer_size = num_samples_per_chan==-1 ? 1024 : num_samples_per_chan
    return Array{precision}(undef, num_channels, buffer_size)
end

#modified from NIDAQ.read
#idea:  could modify this to be mutating (!) function that writes to a view of a bigger array.
#       but warning! the input data vector must be the correct size and contiguous in memory
function read_analog!(data, t::NIDAQ.AITask, precision::DataType, num_samples_per_chan::Integer = -1; timeout::Float64=-1.0)
    num_channels = getproperties(t)["NumChans"][1]
    num_samples_per_chan_read = Int32[0]
    buffer_size = num_samples_per_chan==-1 ? 1024 : num_samples_per_chan
    NIDAQ.catch_error( NIDAQ.read_analog_cfunctions[precision](t.th,
                                                        convert(Int32,num_samples_per_chan),
                                                        timeout,
                                                        reinterpret(Bool32,NIDAQ.Val_GroupByScanNumber),
                                                        Ref(data,1),
                                                        convert(UInt32,buffer_size*num_channels),
                                                        Ref(num_samples_per_chan_read,1),
                                                        reinterpret(Ptr{Bool32},C_NULL)) )
    if num_samples_per_chan_read[1] != num_samples_per_chan
        error("Read $num_samples_per_chan_read samples instead of the $num_samples_per_chan samples requested")
    end
    #data = data[1:num_samples_per_chan_read[1]*num_channels]
    #num_channels==1 ? data : reshape(data, (convert(Int64, num_channels), div(length(data),num_channels)))
    return data
end

function read_digital(t::NIDAQ.DITask, precision::DataType, num_samples_per_chan::Integer = -1; timeout::Float64=-1.0)
    num_channels = getproperties(t)["NumChans"][1]
    num_samples_per_chan_read = Int32[0]
    buffer_size = num_samples_per_chan==-1 ? 1024 : num_samples_per_chan
    data = Array{precision}(buffer_size*num_channels)
    NIDAQ.catch_error( NIDAQ.read_digital_cfunctions[precision](t.th,
        convert(Int32,num_samples_per_chan),
        timeout,
        reinterpret(Bool32,NIDAQ.Val_GroupByScanNumber),
        Ref(data, 1),
        convert(UInt32,buffer_size*num_channels),
        Ref(num_samples_per_chan_read, 1),
        reinterpret(Ptr{Bool32},C_NULL)) )
    if num_samples_per_chan_read[1] != num_samples_per_chan
        error("Read $num_samples_per_chan_read samples instead of the $num_samples_per_chan samples requested")
    end
    #data = data[1:num_samples_per_chan_read[1]*num_channels]
    #num_channels==1 ? data : reshape(data, (convert(Int64, num_channels, div(length(data),num_channels))))
    return data
end

function prepare_ai(coms, nsamps::Integer, bufsz::Int, trigger_terminal::String; terminal_config = "referenced single-ended", clock_source::AbstractString = "")
    print("preparing DAQ for ai...\n")
    rig = rig_name(coms[1])
    usb_req_size = Ref{UInt32}(0)
    chns = map(daq_channel, coms)
    nchans = length(coms)
    tsk = []
    if !in(split(DEVICE_PREFIX, "/")[1], NIDAQ.devices())
        print("trying to throw error\n") #this error isn't visible.  TODO:fix
        error("Device $DEVICE_PREFIX not detected")
    end
    for i = 1:nchans
        if i == 1
            tsk = analog_input(DEVICE_PREFIX * chns[1], terminal_config = terminal_config)
            if nchans == 1
                NIDAQ.setproperty!(tsk, DEVICE_PREFIX * chns[i], "Max", 10.0) #"/ao0:1", "Max", 10.0)
            end
        else
            analog_input(tsk, DEVICE_PREFIX * chns[i])
            NIDAQ.setproperty!(tsk, DEVICE_PREFIX * chns[i], "Max", 10.0) #"/ao0:1", "Max", 10.0)
        end
    end
    if rig == "dummy-6002"
        #note: default Xfer size seems to be 8000 bytes (4000 samples)
        NIDAQ.catch_error(NIDAQ.GetAIUsbXferReqSize(tsk.th, Vector{UInt8}(DEVICE_PREFIX * chns[1]), usb_req_size))
        #See article here http://digital.ni.com/public.nsf/allkb/B7B47F50F9813DFD862575890054EF7C
        usb_req_size = Ref{UInt32}(usb_req_size[] * UInt32(2))
        for i = 1:nchans
            print("    Set Usb transfer request size for AI\n")
            NIDAQ.catch_error(NIDAQ.SetAIUsbXferReqSize(tsk.th, Vector{UInt8}(DEVICE_PREFIX * chns[i]), usb_req_size[]))
        end
    else
        print("Setting DMA Xfer for AI\n")
        for i = 1:nchans
            NIDAQ.catch_error(NIDAQ.SetAIDataXferMech(tsk.th, convert(Ref{UInt8}, Vector{UInt8}(DEVICE_PREFIX * chns[i])), Int32(NIDAQ.DAQmx_Val_DMA)))
        end
    end
    clk_b = UInt8[]
    if isempty(clock_source)
        clk_b = C_NULL
    elseif clock_source[1:2] == "ai"
        clk_b = convert(Ref{UInt8}, b"OnboardClock")
    else
        for ch in clock_source
            push!(clk_b, ch)
        end
        clk_b = convert(Ref{UInt8}, clk_b)
    end
    NIDAQ.catch_error(NIDAQ.CfgSampClkTiming(tsk.th,
            clk_b,
            Float64(ustrip(samprate(first(coms)))),
            NIDAQ.Val_Rising,
            NIDAQ.Val_FiniteSamps,
            UInt64(nsamps)))
    NIDAQ.catch_error(NIDAQ.DAQmxCfgInputBuffer(tsk.th, UInt32(bufsz)))
    print("AI buffer size (nsamps): $bufsz\n")
    if trigger_terminal != "disabled"
        print("    setting up digital start trigger for analog recording...\n")
        props = NIDAQ.getproperties(String(split(DEVICE_PREFIX, "/")[1])) #; warning=true)
        terms = props["Terminals"][1]
        if !in("/" * DEVICE_PREFIX * trigger_terminal, terms)
            error("Terminal $trigger_terminal is does not exist.  Check valid t erminals with NIDAQ.getproperties(dev)[\"Terminals\"]")
        end
        NIDAQ.catch_error(NIDAQ.DAQmxSetStartTrigType(tsk.th, NIDAQ.Val_DigEdge))
        NIDAQ.catch_error(NIDAQ.CfgDigEdgeStartTrig(tsk.th, Ref(codeunits("/" * DEVICE_PREFIX * trigger_terminal),1), NIDAQ.Val_Rising))
    end
    return tsk
end

function prepare_di(coms, nsamps::Integer, bufsz::Int, trigger_terminal::String; clock_source = "")
    error("Not yet implemented")
end

function _record_analog_signals(ai_name::AbstractString, coms::AbstractVector{T}, nsamps::Integer, samps_per_read::Int, trigger_terminal::String, ready_chan::RemoteChannel, clock::AbstractString) where T<:ImagineSignal
    if any(map(isdigital, coms))
        error("Only analog signals are allowed")
    end
    nchans = length(coms)
    #device quirk?
    if isequal(rig_name(coms[1]), "dummy-6002")
        global autostart = isequal(trigger_terminal, "disabled") ? true : false
    else
        global autostart = true
    end
    tsk = prepare_ai(coms, nsamps, 2*samps_per_read, trigger_terminal, clock_source=clock)
    output_array = -1
    fid = -1
    if !isempty(ai_name)
        if isfile(ai_name)
            rm(ai_name)
        end
        open(ai_name, "w+") do fid
            output_array = Mmap.mmap(fid, Matrix{rawtype(first(coms))}, (nchans, nsamps))
        end
    else
        output_array = Matrix{rawtype(first(coms))}(undef, nchans, nsamps)
    end
    record_loop!(output_array, mapper(first(coms)), tsk, nsamps, samps_per_read, autostart, ready_chan)
    sigs = parse_ai(output_array, map(daq_channel, coms), rig_name(first(coms)), samprate(first(coms)))
    return sigs
end

function _record_digital_signals(di_name::AbstractString, coms::AbstractVector{T}, nsamps::Integer, samps_per_read::Int, trigger_terminal::String, clock::AbstractString) where T<:ImagineSignal
    if !all(map(isdigital, coms))
        error("Only digital signals are allowed")
    end
    nchans = length(coms)
    autostart = isequal(trigger_terminal, "disabled") ? true : false
    tsk = prepare_di(coms, nsamps, 2*samps_per_read, trigger_terminal)
    fid = open(di_name, "w+")
    output_array = Mmap.mmap(fid, BitArray, (nchans, nsamps))
    try
        record_loop!(output_array, mapper(first(coms)), tsk, nsamps, samps_per_read, autostart)
    finally
        output_array = 0
        close(fid)
        gc()
    end
    return myid()
end

#Should work for both analog and digital
function record_loop!(output::Matrix{Traw}, m::ImagineInterface.SampleMapper{Traw, TV, TW}, tsk, nsamps::Integer, samps_per_read::Integer, autostart::Bool, ready_chan::RemoteChannel) where {Traw, TV, TW}
    nsamps = size(output,2)
    nchans = size(output,1)
    finished = false
    curstart = 1
    read_incr = samps_per_read - 1
    is_digi = isa(tsk, NIDAQ.DITask)
    if is_digi
        pp = x->Traw(x)
    else #it's an analog voltage
        #TODO: remove max() after handling piezo signal issue
        pp = x-> try
                    ImagineInterface.volts2raw(m)(max(0.0, x)*Unitful.V)
                catch
                    print("error with this sample value: $x \n")
                end
    end
    try
        if autostart
	    print("Auto-starting task\n")
            start(tsk)
        end
        sampbuf = -1
        if !is_digi
            sampbuf = allocate_sampbuf(tsk, Float64, samps_per_read)
        else
            sampbuf = allocate_sampbuf(tsk, UInt8, samps_per_read)
        end
        put!(ready_chan, myid())
        while !finished
            curstop = curstart + read_incr
            if curstop >= nsamps
                finished = true
                curstop = min(nsamps, curstop)
            end
            nsamps_to_read = curstop - curstart + 1
            #now read from buffer
            if !is_digi
                read_analog!(sampbuf, tsk, Float64, nsamps_to_read, timeout=-1.0); #takes almost all the time
                cur_samps = view(sampbuf, :, 1:nsamps_to_read);
                output[:, curstart:curstop] .= pp.(cur_samps);
            else
                output[:, curstart:curstop] .= pp.(read_digital!(sampbuf, tsk, UInt8, nsamps_to_read, timeout=-1.0))
            end
            print("So far $curstop of $nsamps samples have been read...\n")
            curstart = curstop + 1
        end
    catch
        print("some error\n")
        clear(tsk)
        rethrow()
    end
    print("waiting until task is done...\n")
    NIDAQ.catch_error(NIDAQ.WaitUntilTaskDone(tsk.th, 0.0)) #should be done already, throw error if not
    stop(tsk)
    clear(tsk)
    print("Finished.\n")
    return true
end


text:
#=
Type and methods to simplify data sourcing and management of the universe of tradable assets
=#

using ProgressMeter

import Base: show

const SEPARATORS = ['/', '_', '.']

# function guess_tickers(assets::Vector{String})::Vector{Symbol}
#     tickers = Symbol.([Temporal.namefix(split(asset, SEPARATORS)[end]) for asset in assets])
#     @assert tickers == unique(tickers)  "Non-unique ticker symbols found in universe"
#     return tickers
# end

mutable struct Universe
    assets::Vector{String}
    # tickers::Vector{Symbol}
    data::Dict{String,TS}
    from::TimeType
    thru::TimeType
    function Universe(assets::Vector{String}, from::TimeType=Dates.Date(0), thru::TimeType=Dates.today())
        @assert assets == unique(assets)
        # tickers = guess_tickers(assets)
        data = Dict{String,TS}()
        @inbounds for asset in assets
            data[asset] = TS()
        end
        return new(assets, data, from, thru)
    end
end

#TODO: ensure type compatibility across variables (specifically with regard to TimeTypes)
function gather!(universe::Universe; source::Function=Temporal.quandl, verbose::Bool=true)::Nothing
    t0 = Vector{Dates.Date}()
    tN = Vector{Dates.Date}()
    verbose ? progress = Progress(length(universe.assets), 1, "Gathering Universe Data") : nothing
    @inbounds for asset in universe.assets
        verbose ? next!(progress) : nothing
        indata = source(asset)
        push!(t0, indata.index[1])
        push!(tN, indata.index[end])
        universe.data[asset] = indata
    end
    universe.from = max(minimum(t0), universe.from)
    universe.thru = min(maximum(tN), universe.thru)
    return nothing
end

#FIXME: make robust to other time types
function get_overall_index(universe::Universe)::Vector{Date}
    idx = Vector{Date}()
    for asset in universe.assets
        idx = union(idx, universe.data[asset].index)
    end
    return idx
end

function show(io::IO, universe::Universe)
    println(io, "# Universe:")
    for (i, asset) in enumerate(universe.assets  )
        println(io, TAB, "Asset $i:", TAB, asset)
        data = universe.data[asset]
        if isempty(data)
            println(io, TAB, TAB, "(No Data Gathered)")
        else
            println(io, TAB, TAB, "Range:", TAB, data.index[1], " to ", data.index[end])
            println(io, TAB, TAB, "Fields:", TAB, join(String.(data.fields), "  "))
        end
    end
end


text:
module CHSF

using JuLIP
using LinearAlgebra

export chsf, chsf_RADF, chsf_desc_RADF, chsf_desc, chsf_RADF_weights, c_RADFnnl_template, chsf_ijR_template

chsfPI = 3.14159265358979

function chebyshev(x, n)
    T = zeros((n+1, size(x,1)))
    T[1,:] = ones(size(x,1))
    if n > 0
        T[2,:] = x
        for i=3:n+1
            T[i,:] = 2.0 .* x .* T[i-1,:] .- T[i-2,:]
        end
    end
    return T
end

function ch_T(x, n)
    T = zeros(n+1)
    T[1] = 1
    if n > 0
        T[2] = x
        for i=3:n+1
            T[i] = 2.0 * x * T[i-1] - T[i-2]
        end
    end
    return T
end

function cutoff_func(dij, Rc)
    fc = zeros(size(dij))
    dij_in_Rc = dij .< Rc
    fc[dij_in_Rc] = 0.5 .* ( cos.(dij[dij_in_Rc] .* chsfPI/Rc) .+ 1.0)
    return fc
end

function c_RDF(dj, n, Rc)
    wtj = 1 # Single specie case!
    fc_ij = cutoff_func(dj, Rc)
    scaled_r = (2.0 .* dj .- Rc) ./ Rc
    phi_ij = chebyshev(scaled_r, n-1)
    return phi_ij * fc_ij .* wtj
end

function c_ADF(Rj, dj, l, Rc)
    fc_ij = cutoff_func(dj, Rc)
    D = zeros(l)
    for j = 1:length(Rj)-1, k = j+1:length(Rj)
       cos_theta_ijk = dot(Rj[j], Rj[k]) / (dj[j]*dj[k])
       scaled_theta = (2.0 .* cos_theta_ijk .- chsfPI) ./ chsfPI
       Aijk = ch_T(scaled_theta, l-1)
       D .+= Aijk .* fc_ij[j] .* fc_ij[k]
    end
    return D
end

function c_RADFn(Rj, nmax, lmax, Rc)
    dj = norm.(Rj)
    fc_ij = cutoff_func(dj, Rc)
    D = zeros(nmax * lmax)
    for j = 1:length(Rj)-1, k = j+1:length(Rj)
       scaled_rj = (2.0 * dj[j] - Rc) / Rc
       scaled_rk = (2.0 * dj[k] - Rc) / Rc
       Pj = ch_T(scaled_rj, nmax-1) * fc_ij[j]
       Pk = ch_T(scaled_rk, nmax-1) * fc_ij[k]
       cos_theta_ijk = dot(Rj[j], Rj[k]) / (dj[j]*dj[k])
       scaled_theta = (2.0 * cos_theta_ijk - chsfPI) / chsfPI
       Aijk = ch_T(scaled_theta, lmax-1)
       for n = 1:nmax, l = 1:lmax
           D[n + (l-1) * nmax] += Aijk[l] * Pj[n] * Pk[n]
       end
    end
    return D
end

function c_RADFnnl_template(Rj, nmax, lmax, Rc)
    dj = norm.(Rj)
    fc_ij = cutoff_func(dj, Rc)
    rcount = 0
    for j = 1:length(Rj)-1, k = j+1:length(Rj)
        rcount += 1
    end
    D = zeros(rcount, nmax, nmax, lmax)
    rcount = 1
    for j = 1:length(Rj)-1, k = j+1:length(Rj)
       scaled_rj = (2.0 * dj[j] - Rc) / Rc
       scaled_rk = (2.0 * dj[k] - Rc) / Rc
       Pj = ch_T(scaled_rj, nmax-1) * fc_ij[j]
       Pk = ch_T(scaled_rk, nmax-1) * fc_ij[k]
       cos_theta_ijk = dot(Rj[j], Rj[k]) / (dj[j]*dj[k])
       scaled_theta = (2.0 * cos_theta_ijk - chsfPI) / chsfPI
       Aijk = ch_T(scaled_theta, lmax-1)
       for n = 1:nmax, np = 1:nmax, l = 1:lmax
           D[rcount, n, np, l] = Aijk[l] * Pj[n] * Pk[np]
       end
       rcount += 1
    end
    return D
end

function c_RADFnnl(Rj, nmax, lmax, Rc; w=nothing)
    dj = norm.(Rj)
    fc_ij = cutoff_func(dj, Rc)
    D = zeros(nmax, nmax, lmax)
    for j = 1:length(Rj)-1, k = j+1:length(Rj)
       wj = 1.
       wk = 1.
       if w != nothing
          wj = w[j]
          wk = w[k]
       end
       scaled_rj = (2.0 * dj[j] - Rc) / Rc
       scaled_rk = (2.0 * dj[k] - Rc) / Rc
       Pj = ch_T(scaled_rj, nmax-1) * fc_ij[j]
       Pk = ch_T(scaled_rk, nmax-1) * fc_ij[k]
       cos_theta_ijk = dot(Rj[j], Rj[k]) / (dj[j]*dj[k])
       scaled_theta = (2.0 * cos_theta_ijk - chsfPI) / chsfPI
       Aijk = ch_T(scaled_theta, lmax-1)
       for n = 1:nmax, np = 1:nmax, l = 1:lmax
           D[n, np, l] += Aijk[l] * Pj[n] * Pk[np] .* wj .* wk
       end
    end
    return D[:]
end

function c_RADFnnl_weights(DT, W, lenRj, nmax, lmax)
    D = zeros(size(W,1), nmax, nmax, lmax)
    for wi = 1:size(W,1)
       c = 0
       for j = 1:lenRj-1, k = j+1:lenRj
          for n = 1:nmax, np = 1:nmax, l = 1:lmax
             D[wi, n, np, l] += DT[c, n, np, l] .* W[wi, j] .* W[wi, k]
          end
          c += 1
       end
    end
    return D[:]
end

function chsf_desc_RADF(Rs, Rc; n=nothing, l=nothing, np=false, w=nothing)
    if n == nothing
        n = 0
    end
    if l == nothing
        l = 0
    end
    if np
        return c_RADFnnl(Rs, n+1, l+1, Rc, w=w)
    else
        return c_RADFn(Rs, n+1, l+1, Rc)
    end
end

function chsf_desc(Rs, Rc; n=nothing, l=nothing)
    ds = norm.(Rs)
    descriptor = []
    if n != nothing
        if n>-1
            # Concatenate radial basis functions
            descriptor = vcat(descriptor,c_RDF(ds, n+1, Rc))
        end
    end
    if l != nothing
        if l>-1
            # Concatenate angular basis functions
            descriptor = vcat(descriptor,c_ADF(Rs, ds, l+1, Rc))
        end
    end
    return descriptor
end

function chsf_ijR_template(at, Rc, n, l)
    DT = []
    ni = []
    nj = []
    nR = []
    for (i, j, R) in pairs(neighbourlist(at, Rc))
        push!(ni, i)
        push!(nj, j)
        push!(nR, R)
    end
    for i=1:length(at)
        println(i)
        Rs = nR[(ni .== i)]
        push!(DT, c_RADFnnl_template(Rs, n, l, Rc))
    end
    return DT, ni, nj, nR
end

function chsf_RADF_weights(DT, w, at, ni, nj, lenR, n, l)
    representation = []
    for i=1:length(at)
        ns = nj[(ni .== i)]
        ws = [w[j,:] for j in ns].T
        push!(representation, w[i,:] .* c_RADFnnl_weights(DT, ws, lenR, n, l))
    end
    return representation
end

function chsf_RADF(at, Rc; n=nothing, l=nothing, np=true, w=nothing)
    representation = []
    ni = []
    nR = []
    if w != nothing
       nw = []
    end
    for (i, j, R) in pairs(neighbourlist(at, Rc))
        push!(ni, i)
        push!(nR, R)
        if w != nothing
           push!(nw, w[j])
        end
    end
    for i=1:length(at)
        Rs = nR[(ni .== i)]
        if w != nothing
           ws = nw[(ni .== i)]
           push!(representation, w[i] .* chsf_desc_RADF(Rs, Rc, n=n, l=l, np=np, w=ws))
        else
           push!(representation, chsf_desc_RADF(Rs, Rc, n=n, l=l, np=np))
        end
    end
    return representation
end

function chsf(at, Rc; n=nothing, l=nothing)
    representation = []
    ni = []
    nj = []
    nR = []
    nd = []
    for (i, j, R) in pairs(neighbourlist(at, Rc))
        push!(ni, i)
        push!(nj, j)
        push!(nR, R)
        push!(nd, norm(R))
    end
    for i=1:length(at)
        Rs = nR[(ni .== i)]
        push!(representation,chsf_desc(Rs, Rc, n=n, l=l))
    end
    return representation
end

end # module


text:
using DataFrames
using Gadfly
using HDF5, JLD

finetune = load("finetune/statistics.jld")["statistics"]
randinit = load("randinit/statistics.jld")["statistics"]
stats = [finetune, randinit]
category = ["pre-train", "random init"]

plot_vars = (("test-accuracy-accuracy", "Accuracy (Test set)"),
             ("obj-val", "Objective Function (Training set)"))

for (var, title) in plot_vars
  dfs = [begin
    data = stats[i][var]
    x = collect(keys(data)); y = [data[xx] for xx in x]
    DataFrame(x=x, y=y, category=category[i])
  end for i = 1:length(stats)]

  df = vcat(dfs...)
  the_plot = plot(df, x="x", y="y", color="category",
      Geom.line, Guide.xlabel("Iteration"), Guide.title(title))

  draw(SVG("$var.svg", 18cm, 9cm), the_plot)
end


text:
using EHANK

### Scalar error
# Given dt, what size growth rate Q is needed such that
# a kth-order approximation of power series is no longer okay
dt = 1. / 3 # monthly for quarterly
tol = 1e-5
Q_low = 0 # run bisection search
Q_high = 1
kvec = [1; 2; 3]
Q_final = zeros(length(kvec))

for ind = 1:length(kvec)
    Q = 0.
    for i = 1:100000
        Q += .01
        truth = exp(Q * dt)/Q - 1 / Q
        power = kvec[ind]
        approx = dt
        for j = 1:power
            approx += Q^j * dt^(j+1) / factorial(j+1)
        end
        if abs(truth - approx) >= tol
            Q_final[ind] = Q
            break
        end
    end
end

### Matrix error w/KrusellSmith
# Given T, for what size dt is
# a kth-order approximation of power series no longer okay
tol = 1e-8
dt_vec = collect(linspace(1/90, .5, 10))
power_vec = zeros(length(dt_vec))

# Solve model
m = KrusellSmith()
T, ~, ~, ~, ~ = solve(m)

# Bisection search
for ind = 1:length(dt_vec)
    dt = dt_vec[ind]
    power = 0
    while power <= 10
        power += 1
        approx = dt
        for j = 1:power
            approx += T^j * dt^(j+1) / factorial(j+1)
        end
        if maximum(abs.(T^power * dt^(power + 1) / factorial(power + 1))) <= tol
            power_vec[ind] = power-1
            break
        end
        if power == 10
            power_vec[ind] = NaN
        end
    end
end


text:
mutable struct RegressionTable
    headerString::String
    header::Array{String, 2}
    bodies::Vector{Array{String, 2}}
    footerString::String

    function RegressionTable(columns::Int64, headerString::String, header::Array{String, 2}, bodies::Vector{Array{String, 2}}, footerString::String)
        this = new()
        if any([size(body,2) for body in bodies] .!= columns)
            error("Incorrect number of columns in table")
        end
        if size(header, 2) != columns
            error("Header has wrong number of columns")
        end
        if (size(bodies,1) == 0 ) || (size(bodies[1],1)==0)
            error("Table must contain at least one body, and at least one row in the first body.")
        end
        this.headerString = headerString
        this.header = header
        this.bodies = bodies
        this.footerString = footerString
        return this
    end

end

# some helper functions
columns(tab::RegressionTable) = size(tab.bodies[1],2)


text:
using Base.Test
using Distributions

@test convert(Binomial, Bernoulli(0.75)) == Binomial(1, 0.75)

@test convert(Gamma, Exponential(3.0)) == Gamma(1.0, 3.0)
@test convert(Gamma, Erlang(5, 2.0)) == Gamma(5.0, 2.0)


text:
# gen/fetcher.jl
# https://developer.riotgames.com/docs/lol#data-dragon

include(normpath(@__DIR__, "common.jl")) # data_version locales resources

rootdir = normpath(@__DIR__, data_version)
for locale in locales
    datadir = normpath(rootdir, "data", locale)
    !isdir(datadir) && mkpath(datadir)
    for (k, v) in pairs(resources)
        path = normpath(datadir, string(k, ".json"))
        url = "http://ddragon.leagueoflegends.com/cdn/$data_version/data/$locale/$k.json"
        if !isfile(path)
            @info relpath(path, @__DIR__)
            download(url, path)
        end
    end
end


text:
const ESC = Char(0x1B)
const HIDE_CURSOR = ESC * "[?25l"
const SHOW_CURSOR = ESC * "[?25h"
const CLEAR_SCREEN = ESC * "[2J"
const MOVE_CURSOR_TO_ORIGIN = ESC * "[H"
const CLEAR_SCREEN_BEFORE_CURSOR = ESC * "[1J"
const EMPTY_SCREEN = CLEAR_SCREEN_BEFORE_CURSOR * MOVE_CURSOR_TO_ORIGIN

open_maybe(file_name::AbstractString) = open(file_name, "w")
open_maybe(::Nothing) = nothing

close_maybe(io::IO) = close(io)
close_maybe(io::Nothing) = nothing

write_maybe(io::IO, content) = write(io, content)
write_maybe(io::Nothing, content) = 0
write_maybe(io1, io2, content) = write_maybe(io1, content) + write_maybe(io2, content)

function play!(terminal::REPL.Terminals.UnixTerminal, env::AbstractGridWorld, file_name::Union{Nothing, AbstractString}, frame_start_delimiter::Union{Nothing, AbstractString})
    if !isnothing(file_name) && isnothing(frame_start_delimiter)
        error("Please additionally specify a string for the keyword argument `frame_start_delimiter`, for example `frame_start_delimiter = \"FRAME_START_DELIMITER\"`. This is needed when you want to replay the recording.")
    end
    if isnothing(file_name) && !isnothing(frame_start_delimiter)
        error("Please additionally specify a string for the keyword argument `file_name`, for example `file_name = \"recording.txt\"`. This is needed to save the recording in a text file and replay it later.")
    end

    terminal_out = terminal.out_stream
    terminal_in = terminal.in_stream
    file = open_maybe(file_name)

    action_keys = get_action_keys(env)
    key_bindings = "Key bindings to play: 'q': quit, 'r': GW.reset!, $(action_keys): GW.act!"

    char = nothing

    write(terminal_out, CLEAR_SCREEN)
    write(terminal_out, MOVE_CURSOR_TO_ORIGIN)
    write(terminal_out, HIDE_CURSOR)

    REPL.Terminals.raw!(terminal, true)

    try
        while true
            write_maybe(file, frame_start_delimiter)

            frame = key_bindings
            frame = frame * "\n" * "Last play character read: $(char)"
            frame = frame * "\n" * repr(MIME"text/plain"(), env)

            write_maybe(terminal_out, file, frame)

            char = read(terminal_in, Char)

            if char == 'q'
                write(terminal_out, SHOW_CURSOR)
                close_maybe(file)
                REPL.Terminals.raw!(terminal, false)
                return nothing
            elseif char == 'r'
                reset!(env)
            elseif char in action_keys
                act!(env, findfirst(==(char), action_keys))
            end

            write(terminal_out, EMPTY_SCREEN)
        end
    finally
        write(terminal_out, SHOW_CURSOR)
        close_maybe(file)
        REPL.Terminals.raw!(terminal, false)
    end

    return nothing
end

play!(env::AbstractGridWorld; file_name = nothing, frame_start_delimiter = nothing) = play!(REPL.TerminalMenus.terminal, env, file_name, frame_start_delimiter)

function replay(terminal::REPL.Terminals.UnixTerminal, file_name::AbstractString, frame_start_delimiter::AbstractString, frame_rate = Union{Nothing, Real})
    terminal_out = terminal.out_stream
    strings = split(read(file_name, String), frame_start_delimiter)
    frames = @view strings[2:end]
    num_frames = length(frames)

    write(terminal_out, CLEAR_SCREEN)
    write(terminal_out, MOVE_CURSOR_TO_ORIGIN)
    write(terminal_out, HIDE_CURSOR)

    if isnothing(frame_rate)
        terminal_in = terminal.in_stream
        replay_key_bindings = "Key bindings to replay: 'q': quit, 'f': first frame, 'n': next frame, 'p': previous frame"
        current_frame = 1
        char = nothing

        REPL.Terminals.raw!(terminal, true)

        try
            while true
                replay_frame = replay_key_bindings
                replay_frame = replay_frame * "\n" * "Last replay character read: $(char)"
                replay_frame = replay_frame * "\n" * "------------FRAME_START------------"
                replay_frame = replay_frame * "\n" * frames[current_frame]

                write(terminal_out, replay_frame)

                char = read(terminal_in, Char)

                if char == 'q'
                    write(terminal_out, SHOW_CURSOR)
                    REPL.Terminals.raw!(terminal, false)
                    return nothing
                elseif char == 'f'
                    current_frame = 1
                elseif char == 'n'
                    current_frame = mod1(current_frame + 1, num_frames)
                elseif char == 'p'
                    current_frame = mod1(current_frame - 1, num_frames)
                end

                write(terminal_out, EMPTY_SCREEN)
            end
        finally
            write(terminal_out, SHOW_CURSOR)
            REPL.Terminals.raw!(terminal, false)
            return nothing
        end
    else
        write(terminal_out, CLEAR_SCREEN)
        write(terminal_out, MOVE_CURSOR_TO_ORIGIN)
        write(terminal_out, HIDE_CURSOR)

        for frame in frames
            write(terminal_out, frame)
            sleep(1 / frame_rate)
            write(terminal_out, EMPTY_SCREEN)
        end

        write(terminal_out, SHOW_CURSOR)
        return nothing
    end
end

replay(; file_name, frame_start_delimiter, frame_rate = nothing) = replay(REPL.TerminalMenus.terminal, file_name, frame_start_delimiter, frame_rate)


text:
using ReferenceTests, Test
include("../scripts-part2/part2.jl")

opt = SimIn_t()
opt.nx, opt.ny = 257, 65
opt.tol = 1.0e-12
opt.W_init_strategy = W_from_file
sim_results = navier_stokes_2D(; opt=opt, verbose=false, do_vis=false, testmode=true)

comp_size(lhs, rhs) = size(lhs) == size(rhs)

function comp_val(nx, ny; atol=1e-8)
    comp_val_inner(lhs::Matrix{Float64}, rhs::Matrix{Float64}) =
        all(abs.(lhs[2:nx-1, 2:ny-1] - rhs[2:nx-1, 2:ny-1]) .< atol)
    return comp_val_inner
end

function my_test_reference(binfile, arr; comp=(==))
    bin_data = open(binfile, "r") do io
        return load(io)
    end
    @test comp(bin_data, arr)
end

# Add unit and reference tests
@testset "Test against FORTRAN reference implementation: array sizes." begin
    my_test_reference(joinpath(@__DIR__, "reftest-files/fortran/T.bin"), sim_results.T, comp=comp_size)
    my_test_reference(joinpath(@__DIR__, "reftest-files/fortran/W.bin"), sim_results.W, comp=comp_size)
    my_test_reference(joinpath(@__DIR__, "reftest-files/fortran/S.bin"), sim_results.S, comp=comp_size)
end

@testset "Test against FORTRAN reference implementation: array values." begin
    atol = 1e-8
    my_test_reference(joinpath(@__DIR__, "reftest-files/fortran/T.bin"), sim_results.T, comp=comp_val(opt.nx, opt.ny, atol=atol))
    my_test_reference(joinpath(@__DIR__, "reftest-files/fortran/W.bin"), sim_results.W, comp=comp_val(opt.nx, opt.ny, atol=atol))
    my_test_reference(joinpath(@__DIR__, "reftest-files/fortran/S.bin"), sim_results.S, comp=comp_val(opt.nx, opt.ny, atol=atol))
end

@ParallelStencil.reset_parallel_stencil()


text:
using Base.Test
using MPFI

include("MPFI_test.jl")


text:
function get_size(obj::BoundarySizeVisitor)
    return jcall(obj, "getSize", jdouble, ())
end

function visit_internal_node(obj::BoundarySizeVisitor, arg0::BSPTree)
    return jcall(obj, "visitInternalNode", void, (BSPTree,), arg0)
end

function visit_leaf_node(obj::BoundarySizeVisitor, arg0::BSPTree)
    return jcall(obj, "visitLeafNode", void, (BSPTree,), arg0)
end

function visit_order(obj::BoundarySizeVisitor, arg0::BSPTree)
    return jcall(obj, "visitOrder", BSPTreeVisitor_Order, (BSPTree,), arg0)
end


text:
immutable OVCallbacks
    read_func::Ptr{Void}
    seek_func::Ptr{Void}
    close_func::Ptr{Void}
    tell_func::Ptr{Void}
end

immutable OggVorbisFile
    datasource::Ptr{Void}
    seekable::Cint
    offset::Int64
    nd::Int64
    oy::Ogg.SyncState
    links::Cint
    offsets::Ptr{Int64}
    dataoffsets::Ptr{Int64}
    serialnos::Ptr{Clong}
    pcmlengths::Int64
    vi::Ptr{VorbisInfo}
    vc::Ptr{VorbisComment}
    pcm_offset::Int64
    ready_state::Cint
    current_serialno::Clong
    current_link::Cint
    bittrack::Cdouble
    samptrack::Cdouble
    os::Ogg.StreamState
    vd::VorbisDSPState
    vb::VorbisBlock
    callbacks::OVCallbacks
end

clear(vf::OggVorbisFile) =
    ccall((:ov_clear,vorbisfile),Cint,(Ref{OggVorbisFile},),vf)
fopen(path::AbstractString,vf::OggVorbisFile) =
    ccall((:ov_fopen,vorbisfile),Cint,(Cstring,Ref{OggVorbisFile}),path,vf)


text:
module Users

using SearchLight, SearchLight.Validation, UsersValidator

export User

mutable struct User <: AbstractModel
  ### INTERNALS
  _table_name::String
  _id::String
  _serializable::Vector{Symbol}

  ### FIELDS
  id::DbId
  username::String
  password::String
  name::String
  email::String

  ### VALIDATION
  validator::ModelValidator

  ### CALLBACKS
  # before_save::Function
  # after_save::Function
  # on_save::Function
  # on_find::Function
  # after_find::Function

  ### SCOPES
  # scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

  ### constructor
  User(;
    ### FIELDS
    id = DbId(),
    username = "",
    password = "",
    name = "",
    email = "",

    ### VALIDATION
    validator = ModelValidator([
      ValidationRule(:username, UsersValidator.not_empty)
    ]),

    ### CALLBACKS
    # before_save = (m::Todo) -> begin
    #   @info "Before save"
    # end,
    # after_save = (m::Todo) -> begin
    #   @info "After save"
    # end,
    # on_save = (m::Todo, field::Symbol, value::Any) -> begin
    #   @info "On save"
    # end,
    # on_find = (m::Todo, field::Symbol, value::Any) -> begin
    #   @info "On find"
    # end,
    # after_find = (m::Todo) -> begin
    #   @info "After find"
    # end,

    ### SCOPES
    # scopes = Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}()

  ) = new("users", "id", Symbol[],                                                      ### INTERNALS
          id, username, password, name, email,                                          ### FIELDS
          validator,                                                                    ### VALIDATION
          # before_save, after_save, on_save, on_find, after_find                       ### CALLBACKS
          # scopes                                                                      ### SCOPES
          )
end

end

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text: # no tests


text:
# include('gridworld.jl')
# todo, remove vec2, makes moves really annoying to perform
struct Vec2
	x::UInt8
	y::UInt8
end

function grid_init(dims::Vec2)
	grid = zeros(UInt8, (dims.x, dims.y))
	return grid
end

function player_init(grid::Array{UInt8, 2}, pos::Vec2)
	grid[pos.x, pos.y] = UInt8(1)
	return grid
end

function move(grid::Array{UInt8, 2}, pos::Array{UInt8, 1}, m::Char)
	# todo valid move checker
	#new_x = MOVES[m][1] + pos.x
	#new_y = MOVES[m][2] + pos.y
	#new_pos = UInt8[new_x, new_y]
	new_pos = MOVES[m] +
	if
	return (new_x, new_y)
end

function legal_moves(grid_dims::Tuple, pos::Array{UInt8, 1}, m::Char)
	# this is garbage logic code. could be way cleaner
	if pos[1] == 1 & m == 'u'
		#invalid
		print("invalid move")
	else if pos[1] == grid_dims[1] & m == 'd'
		#invalid
		print("invalid move")
	else if pos[2] == 1 & m == 'l'
		print("invalid move")
	else if pos[2] == grid_dims[2] & m == 'r'
		print("invalid move")
	else
		print("valid")
	end
end

# function is_valid(grid_dims::Tuple, pos::Array{UInt8, 1}

dims = Vec2(10, 20)
pos = Vec2(5, 5)

# fukvec2
player_pos = UInt8[5, 5]

inv_case1 = UInt8[1, 2]
inv_case2 = UInt8[10, 2]
inv_case3 = UInt8[2, 1]
inv_case4 = UInt8[2, 20]

test_cases = [inv_case1, inv_case2, inv_case3, inv_case4]
test_moves = ['u', 'd', 'l', 'r']

# tuple? grid_dims = UInt8[
grid_dims = (10, 20)

g = grid_init(dims)
print(g)

g = player_init(g, pos)
println(g)

# 'd' = [1, 0], 'u' = [-1, 0]
# 'l' = [0, -1], 'r' = [0, 1]

MOVES = Dict('d' => Int8[1, 0], 'u' => Int8[-1, 0], 'l' => Int8[0, -1], 'r'=>Int8[0, 1])
println(MOVES)

after_m = move(g, pos, 'd')
println(after_m)
function test_legal_moves(grid_dims, cases, moves)
	if length(cases) != length(moves)
		println("cases len must == moves len")
		return 0
	end
	for i in 1:length(cases)
		legal_moves(grid_dims,
legal_moves(grid_dims, inv_case1


text:
using HelloWorld
using Test

@testset "HelloWorld.jl" begin
    # Write your own tests here.
end


text:
struct DataSet1
    x::Float64
    y::Float64
    z::Float64
end

mutable struct DataSet2
    a::Array{Float64}
    b::Array{Float64}
    c::Array{Float64}
end

struct InitValues
    init_a::Vector
    init_b::Vector
    init_c::Vector
end

mutable struct DekaiStruct

    data1::DataSet1
    data2::DataSet2

end

DekaiStruct(param::Tuple{Float64, Float64, Float64}, initvalues::InitValues, datanum::Int64) = begin

    data1 = DataSet1(param[1], param[2], param[3])

    a_array = zeros(10, datanum)
    a_array[:, begin] = initvalues.init_a

    b_array = zeros(100, datanum)
    b_array[:, begin] = initvalues.init_b

    c_array = zeros(100, datanum)
    c_array[:, begin] = initvalues.init_c

    data2 = DataSet2(a_array, b_array, c_array)

    return DekaiStruct(data1, data2)
end

initvalues = InitValues(rand(10), ones(100), randn(100))

dekaiyatu = DekaiStruct((1.0,2.0,3.0), initvalues, 1000)