text:
import Calculus: simplify
import DataStructures: OrderedSet, OrderedDict
using Distributions
# import Parameters: with_kw

# =================================
# main codegen functions and macros
# =================================

function sde_struct(typename::Symbol, supertype::Symbol, d::Integer, m::Integer, parameter_vars, docstring)
  block = Expr(:block, [:($p::Float64) for p in parameter_vars]...)
  typedef =  Expr(:struct, false, :($typename <: SDEModels.$supertype{$d,$m}), block)
  # if length(parameter_vars) > 0
     # typedef = with_kw(typedef, current_module())
  # end
  quote
    # @doc $docstring ->
    $typedef
    export $typename
  end
end

# eltype_promote(::Type{T1}, val::T2)  where {T1<:Number,T2<:Number} =
#   convert(promote_type(T1,T2), val)
#
# eltype_promote(::Type{T1}, val::SVector{S,T2}) where {T1<:Number,T2<:Number,S} =
#   convert(SVector{S,promote_type(T1,T2)}, val)
#
# eltype_promote(::Type{T1}, val::SMatrix{S1,S2,T2}) where {T1<:Number,T2<:Number,S1,S2} =
#   convert(SMatrix{S1,S2,promote_type(T1,T2)}, val)

function marked_sde_state_function(typename::Symbol, functionname::Symbol, model_vars, parameter_vars, mark_vars, ex)
  docstring = "$typename: $ex"
  replacements = Dict()
  # if length(model_vars) == 1
    # push!(replacements, first(model_vars) => :x)
  # else
    merge!(replacements, Dict(j => :(x[$i]) for (i,j) in enumerate(model_vars)))
  # end
  merge!(replacements, Dict(broadcast(s -> s => :(model.$s), parameter_vars)))
  if length(mark_vars) == 1
    push!(replacements, first(mark_vars) => :ξ)
  else
    merge!(replacements, Dict(j => :(ξ[$i]) for (i,j) in enumerate(mark_vars)))
  end
  m = length(model_vars)
  ex = replace_symbols(ex, replacements)
  quote
    # @doc $docstring ->
    function (SDEModels.$functionname)(model::$typename, t::Number, x::S, ξ) where S
      # SDEModels.eltype_promote(eltype(S), $ex)
      $ex
    end
  end
end

function sde_state_function(typename::Symbol, functionname::Symbol, model_vars, parameter_vars, ex)
  docstring = "$typename: $ex"
  replacements = Dict()
  # if length(model_vars) == 1
    # push!(replacements, first(model_vars) => :x)
  # else
    merge!(replacements, Dict(j => :(x[$i]) for (i,j) in enumerate(model_vars)))
  # end
  merge!(replacements, Dict(broadcast(s -> s => :(model.$s), parameter_vars)))
  m = length(model_vars)
  ex = replace_symbols(ex, replacements)
  quote
    # @doc $docstring ->
    function (SDEModels.$functionname)(model::$typename, t::Number, x::S) where S
      # SDEModels.eltype_promote(eltype(S), $ex)
      $ex
    end
  end
end

function sde_model_function(typename::Symbol, functionname::Symbol, parameter_vars, ex)
  replacements = Dict(broadcast(s -> s => :(model.$s), parameter_vars))
  ex = replace_symbols(ex, replacements)
  quote
    function (SDEModels.$functionname)(model::$typename)
      $ex
    end
  end
end

function sde_model(typename::Symbol, ex::Expr, parameter_vars...)
  if ex.head == :block
    equations = filter(e -> e isa Expr && (e.head == :(=) || (e.head == :call && e.args[1] == :~)), ex.args)
  elseif ex.head == :(=)
    equations = [ex]
  else
    error("Expression must be block or assignment $(ex.head)")
  end

  correlation_equations = filter(eq -> eq.head == :(=) && isa(eq.args[1], Expr), equations)
  mark_equations = filter(eq -> eq.head == :call && eq.args[1] == :~, equations)
  equations = filter(eq -> eq.head == :(=) && isa(eq.args[1], Symbol), equations)

  model_vars = foldl(merge!, matchdict(r"(?<=^d).*", e.args[1]) for e in equations)
  time_vars = OrderedDict([:dt => :t])
  process_vars = foldl(merge!, matchdict(r"(?<=^d)[wW].*", e.args[2]) for e in equations)
  jump_vars = foldl(merge!, matchdict(r"(?<=^d)[N]", e.args[2]) for e in equations)
  differentials = union(keys(time_vars), keys(process_vars), keys(jump_vars))
  correlation_matrix = correlation_matrix_extract(correlation_equations, process_vars)
  correlation_chol = symbolic_chol(correlation_matrix)
  drift_expressions = [factor_extract(e.args[2], :dt, differentials) for e in equations]
  drift = cat_expressions(drift_expressions)
  # drift_jacobian_expressions = [differentiate(f_i, x_j) for f_i in drift_expressions, x_j in values(model_vars)]
  # drift_jacobian = isempty(drift_jacobian_expressions) ? 0 : cat_expressions(drift_jacobian_expressions)
  diffusion_expressions = [factor_extract(e.args[2], dw, differentials) for e in equations, dw in keys(process_vars)]
  diffusion = isempty(diffusion_expressions) ? 0 : cat_expressions(symbolic_mul(diffusion_expressions, correlation_chol))
  jump_expressions = [factor_extract(e.args[2], dn, differentials) for e in equations, dn in keys(jump_vars)]
  jump = isempty(jump_expressions) ? 0 : cat_expressions(jump_expressions[:,1]) # currently only allows one jump
  mark_vars = isempty(mark_equations) ? Symbol[] : foldl(merge!, [e.args[2]] for e in mark_equations)
  if isempty(parameter_vars)
    parameter_vars = setdiff(union(symbols(drift), symbols(diffusion), symbols(jump)), union(values(model_vars), [:t], values(mark_vars)))
  end

  if length(jump_vars) > 0
    supertype = :JumpSDE
  elseif isempty(intersect(symbols(diffusion), values(model_vars)))
    supertype = :StateIndependentDiffusion
  else
    supertype = :AbstractSDE
  end

  docstring = """
    Model variables: $(join(values(model_vars), ", "))

    Process variables: $(join(values(process_vars), ", "))

    Parameter variables: $(join(parameter_vars, ", "))

    Jump variables: $(join(values(jump_vars), ", "))

    Definition:

    $(join(string.(equations), "\n\n"))
  """
  blk = Expr(:block)
  append!(blk.args, sde_struct(typename, supertype, length(equations), length(process_vars), parameter_vars, docstring).args)
  append!(blk.args, sde_state_function(typename, :drift, values(model_vars), parameter_vars, drift).args)
  # append!(blk.args, corrected_drift_function(typename, collect(values(model_vars)), parameter_vars, drift_expressions, diffusion_expressions).args)
  # append!(blk.args, sde_state_function(typename, :drift_jacobian, values(model_vars), parameter_vars, drift_jacobian).args)
  append!(blk.args, sde_state_function(typename, :diffusion, values(model_vars), parameter_vars, diffusion).args)
  append!(blk.args, marked_sde_state_function(typename, :jump, values(model_vars), parameter_vars, mark_vars, jump).args)
  # append!(blk.args, sde_model_function(typename, :variables, parameter_vars, :($(values(model_vars)...))).args)
  append!(blk.args, sde_model_function(typename, :variables, parameter_vars, Tuple(values(model_vars))).args)
  if supertype == :JumpSDE
    mark_dist = mark_equations[1].args[3]
    append!(blk.args, sde_model_function(typename, :mark_distribution, parameter_vars, mark_dist).args)
  end
  blk
end

macro sde_model(typename, ex, parameter_vars...)
  esc(sde_model(typename, ex, parameter_vars...))
end

# ================================
# generic codegen helper functions
# ================================

function matchdict(r::Regex, ex)
  syms = symbols(ex)
  dict = OrderedDict()
  for sym in syms
    str = string(sym)
    m = match(r, str)
    if m != nothing
      push!(dict, sym => Symbol(m.match))
    end
  end
  dict
end

function factor_extract(ex, one_sym::Symbol, zero_syms)
  replacements = Dict(s => 0 for s in zero_syms)
  replacements[one_sym] = 1
  factor = simplify(replace_symbols(ex, replacements))
  # isa(factor, Number) ? Float64(factor) : factor
end

cat_expressions(x::Array{T,1}) where T =
  length(x) == 1 ? x[1] : :(StaticArrays.SVector{$(size(x)...)}($(x...)))

cat_expressions(x::Array{T,2}) where T =
  length(x) == 1 ? x[1] : :(StaticArrays.SMatrix{$(size(x)...)}($(x...)))

replace_symbols(ex, dict) = replace_symbols!(copy(ex), dict)
replace_symbols(sym::Symbol, dict) = replace_symbols!(sym, dict)
replace_symbols(tup::Tuple, dict) = Tuple(replace_symbols(t, dict) for t in tup)
replace_symbols!(ex, dict::AbstractDict) = haskey(dict, ex) ? dict[ex] : ex
function replace_symbols!(ex::Expr, dict::AbstractDict)
  for i in 1:length(ex.args)
    ex.args[i] = replace_symbols!(ex.args[i], dict)
  end
  ex
end

symbols(ex) = symbols!(OrderedSet{Symbol}(), ex)
symbols!(set::OrderedSet{Symbol}, ex) = set
symbols!(set::OrderedSet{Symbol}, sym::Symbol) = push!(set, sym)
function symbols!(set::OrderedSet{Symbol}, ex::Expr)
  start = ex.head == :call ? 2 : 1
  for arg in ex.args[start:end]
    symbols!(set, arg)
  end
  set
end

function symbolic_mul(A, B)
  C = Array{Any}(undef, size(A, 1), size(B, 2))
  for i=1:size(A, 1), j=1:size(B, 2)
    C[i,j] = simplify(Expr(:call, :+, [:($(A[i,k])*$(B[k,j])) for k=1:size(A, 2)]...))
  end
  C
end

function symbolic_chol(A)
  # note: this assumes A is symmetric
  n = size(A, 1)
  B = Array{Any}(undef, n, n)
  for k=1:n
    tmp = simplify(Expr(:call, :-, A[k,k], Expr(:call, :+, 0, [:($(B[k,j])^2) for j=1:k-1]...)))

    B[k,k] = tmp == 1 ? 1 : simplify(Expr(:call, :sqrt, tmp))
    for i=1:n
      if i == k
        continue
      end
      B[i,k] = simplify(Expr(:call, :/, Expr(:call, :-, A[i,k], Expr(:call, :+, 0, [:($(B[i,j])*$(B[k,j])) for j=1:k-1]...)), B[k,k]))
    end
  end
  B
end

# =============================
# specialized codegen functions
# =============================

function correlation_matrix_extract(equations, process_vars)
  process_var_keys = collect(keys(process_vars))
  n = length(process_var_keys)

  S = Any[i==j ? 1 : 0 for i=1:n, j=1:n]
  for e in equations
    lhs = e.args[1]
    rhs = e.args[2]
    if rhs.head == :block # if lhs is not a symbol, rhs appears to be a block. parsed as infix function definition?
      rhs = rhs.args[end]
    end
    if lhs.args[1] == :* && length(lhs.args) == 3
      i = findfirst(isequal(lhs.args[2]), process_var_keys)
      j = findfirst(isequal(lhs.args[3]), process_var_keys)
      ρ = factor_extract(rhs, :dt, [])
      S[i,j] = S[j,i] = simplify(ρ)
    end
  end
  S # note: S is symmetric
end

# function corrected_drift_function(typename::Symbol, model_vars, parameter_vars, drift_equations, diffusion_equations)
#   D = length(model_vars)
#   corrected = Array{Any}(undef, D)
#   for k in 1:D
#     terms = [Expr(:call, :*, f_i, differentiate(f_i, model_vars[k])) for f_i in diffusion_equations[k,:]]
#     sum_terms = Expr(:call, :+, terms...)
#     corrected[k] = simplify(:($(drift_equations[k]) - correction * $sum_terms))
#   end
#   ex = cat_expressions(corrected)
#   replacements = Dict()
#   # if length(model_vars) == 1
#     # push!(replacements, first(model_vars) => :x)
#   # else
#     merge!(replacements, Dict(j => :(x[$i]) for (i,j) in enumerate(model_vars)))
#   # end
#   merge!(replacements, Dict(broadcast(s -> s => :(model.$s), parameter_vars)))
#   ex = replace_symbols(ex, replacements)
#   quote
#     function SDEModels.corrected_drift(model::$typename, t::Number, x::S, correction::Number) where S
#       # SDEModels.eltype_promote(eltype(S), $ex)
#       $ex
#     end
#   end
# end


text:
"""
    AlleeExtinction(minfounders)
    AlleeExtinction(; minfounders=5.0)
    AlleeExtinction{R,W}(minfounders)

Causes extinction in a cell when a population is below the
minimum number of individuals required to maintain it.

Pass grid name `Symbol`s to `R` and `W` type parameters to use specific grids.

$(FIELDDOCTABLE)
"""
@columns struct AlleeExtinction{R,W,MF} <: CellRule{R,W}
    # Field         | Default | Flatten | Limits       | Description
    minfounders::MF | 5.0     | true    | (1.0, 200.0) | "Minimum founding individuals required to to start an ongoing population"
end

@inline applyrule(data, model::AlleeExtinction, state, args...) =
    (state >= model.minfounders ? state : zero(state))

ext: jl


text:
######### CmdStan sample example  ###########

using StanSample

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Keep tmpdir to prevent re-compilation
stanmodel = SampleModel(
  "bernoulli", bernoulli_model;
  tmpdir = tmpdir,
  method = StanSample.Sample(adapt=StanSample.Adapt(delta=0.85))
)


text:
using Pkg

Pkg.add(
    url="https://github.com/terasakisatoshi/libhello_jll.jl.git",
    rev="libhello-v0.2.0+0"
)


text:
include(joinpath(@__DIR__, "../../src/sketches/FindModel.jl"))
include(joinpath(@__DIR__, "FLS.jl"))  # where the sketches are defined

using Test
using Catlab.CategoricalAlgebra
using Catlab.Graphs.BasicGraphs: ReflexiveGraph, TheoryReflexiveGraph

met = catpres_to_linear(TheoryReflexiveGraph);
r1, r2 = ReflexiveGraph(2), ReflexiveGraph(2)
set_subpart!(r2,:src,[1,1])
@test length(check_diagram(met, r1, :V)) == 2
@test length(check_diagram(met, r2, :V)) == 1

@test check_diagrams(met, r1)
@test_throws AssertionError("Equations for V not satisfied by indices [2]"
                           ) check_diagrams(met, r2)

#############################
# LEFT INVERSE / INVOLUTION #
#############################

srch_leftinv = c->[crels_to_cset(leftinv, r) for r in search(leftinv, c)]

for (x, y) in [[2,2]=>2, # inv can be swap or id, f is forced to be an iso
               [2,1]=>0, # injective function to smaller card
               [1,2]=>2, # inv can be swap or id, f's choices are symmetric.
               [2,3]=>4] # see note below
    @test length(search(leftinv, x)) == y
end

# There are 4 nonisomorphic options for A=2,B=3:
# inv can be id OR swap 2 of the three elements in B.
# There is only one automorph where inv is id, and
# if B swaps, then f(1) and f(2) can map to the swapped
# elements, in which case there is just one automorph.
# There are two automorphs in the case where only one out of
# f(1),f(2) is in the swapped pair in B. These are
# distinguished by whether or not the element of g that is
# not mapped to maps to the A whose image in in the swap or not.

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

#############
# EQUALIZER #
#############

srch_equalizer = c->[crels_to_cset(equalizer, r)
                     for r in search(equalizer, c)]

# The pullback is empty because f=[1,1] and g=[2,2]
@test length(search(equalizer,[2,2,0])) == 1

# There are 4 models.
# f=[1,2], g=[1,2] means the pullback is {(1,1),(2,2)}
# f=[1,2], g=[1,1] means the pullback is {(1,1),(1,2)}
# f=[1,1], g=[1,2] means the pullback is {(1,1),(2,1)}
# f=[2,1], g=[1,2] means the pullback is {(1,2),(2,1)}
@test length(search(equalizer,[2,2,2]))==4

# If f,g map to same element, then pullback is A×A
@test length(search(equalizer,[2,2,4]))==1

#################################
# SET PERMUTATIONS / PARTITIONS #
#################################

srch_perm = c -> [crels_to_cset(perm, r) for r in search(perm, [c])]

# Isomorphism classes of permutations AKA PARTITIONS
# [1]=id, [2]= adds swap, [3] adds a 3-cycle,
# [4] adds a 4-cycle + having two swaps
iso_perms = [1,2,3,5]
for (i, n_perms) in enumerate(iso_perms)
    @assert length(search(perm, [i])) == n_perms
end

####################
# REFLEXIVE GRAPHS #
####################

srch_reflg = cs -> [crels_to_cset(reflg, r) for r in search(reflg, cs)]

for (consts, expected) in [
    [1,1]=>1,   # Single node with self loop
    [0,1]=>0,   # arrow needs a vertex
    [1,0]=>0,   # vertex needs an id arrow
    [2,2]=>1,   # both arrows must be id
    [2,3]=>2,   # the non-id arrow can bridge or not
    [2,4]=>(2+  # 2 loops (either same vertex or not),
            2+  # 1 loop (either ↺← or ↺→)
            2)] # 0 loops (either ⇆ or ⇉)
    @test length(search(reflg,consts)) == expected
end

###########################
# CONSTANT (arrow from 1) #
###########################

srch_const2 = cs -> [crels_to_cset(const2, r) for r in search(const2, cs)]

# Only one option if codomain is singleton
@test length(search(const2, [1,1])) == 1

# with two constants, they can either be the same or different
for i in 2:5
    @test length(search(const2, [1,i])) == 2
end

# Nothing for the terminal object to map to
@test length(search(const2, [1,0])) == 0

# Terminal object must be singleton, no more, no less
for i in 0:5
    @test length(search(const2, [0,i])) == 0
    @test length(search(const2, [2,i])) == 0
end

####################################
# Infinite sequence of a's and b's #
####################################

srch_inflist = cs -> [crels_to_cset(inflist, r) for r in search(inflist, cs)]
srch_inflist([1,0,0])

for combo in cartprod(0:3,0:3,0:3)
    n1, nd, nl = combo
    # terminal object must have 1
    # 1 must have something to send a,b to
    # no results if either of these don't happen
    if n1 != 1 || nd == 0
        expected = 0
    # |L|=|L|x|D| satisfied trivially when |L|=0
    # just have the number of ways for a,b to map to D
    elseif nl == 0 #
        expected = nd == 1 ? 1 : 2 # either to same elem or diff, if |D| > 1
    elseif nd == 1
        expected = iso_perms[nl]
    else
        expected = 0 # nd > 1, and nl >0, so |L|=|L|x|D| impossible to solve
    end
    length(search(inflist,collect(combo))) == expected
end

srch_trips = cs -> [crels_to_cset(trips, r) for r in search(trips, cs)]

##############
# SEMIGROUPS #
##############

srch_sg = cs -> [crels_to_cset(semig, r) for r in search(semig, cs)]
# too big to solve [2,4,8].

##############
# Categories #
##############
"""
this was defined using an old implementation that used graph homomorphisms,
should be redone using LabeledGraphs
#------------------------------------------------
# 10.1.5: Categories

# note composition is defined "backwards", i.e. g∘f
CatG = Graph(4) # ob arr composable_arrows composable_triples
o,a,a²,a³ = 1:4 # shorthand
add_edges!(CatG, [o,a,a,a²,a²,a²,a³,a³,a³,a³, a³, a³, a³, a,    a,   a],
                 [a,o,o,a, a, a, a, a, a, a², a², a², a², a²,   a²,  a])
                # u s t π₁ π₂ c  p₁ p₂ p₃ p₁₂ p₂₃ p₁c cp₃ utid idus, ida
u,s,t,π₁,π₂,c,p₁,p₂,p₃,p₁₂,p₂₃,p₁c,cp₃,utid,idus,ida,v = 1:17
# There's a clear typo in the second diagram in CTCS
CD = T[
    T(OutwardTri,     CatG; V=[a³,a²,a,a],    E=[p₁₂,p₁,p₂,π₁,π₂]),
    T(OutwardTri,     CatG; V=[a³,a²,a,a],    E=[p₂₃,p₂,p₃,π₁,π₂]),
    T(SquareTriangle, CatG; V=[a³,a²,a,a²,a], E=[cp₃,p₃,p₁₂,π₂,π₁,c]),
    T(SquareTriangle, CatG; V=[a³,a²,a,a²,a], E=[p₁c,p₁,p₂₃,π₁,π₂,c]),
    T(SquareTriangle, CatG; V=[a,a²,a,o,a],   E=[utid,ida,t,π₂,π₁,u]),
    T(SquareTriangle, CatG; V=[a,a²,a,o,a],   E=[idus,ida,s,π₁,π₂,u]),
    T(InvertedTri,    CatG; V=[a,a²,a,a],     E=[idus,ida,utid,ida,c]), # unit
    T(Square,         CatG; V=[a³,a²,a²,a],   E=[p₁c,cp₃,c,c]) # associativity
]

CLim1 = T(Square, CatG; V=[a²,a,a,o], E=[π₁,π₂,s,t])
CLim2G = Graph(6)
add_edges!(CLim2G, [1,1,1,2,3,3,4],[2,3,4,5,5,6,6])
CLim2 = T(CLim2G, CatG; V=[a³,a,a,a,o,o], E=[p₁,p₂,p₃,s,t,s,t])

CatSketch= FLSketch(CatG, CD, Pair{Int,T}[1 => CLim1, 1 => CLim2])
"""


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.
#*-----------------------------------------------------------------------
#* Author: Tiago Os򱨯rio
#* $Id$
#*---------------------------------------------------------------------
type IAE
	IAE()=begin
		[
			:(output=0*"s"),
		],
		[
			"",
		],
		new(
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"input signal"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"output signal",
				:Unit=>"s"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"setpoint"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"absolute error"
			)),
			[
				:(absError = abs(setPoint-input)),
				:(diff(output)=absError),
			],
			[
				"Absolute error definition","Calculate output",
			],
			[:input,:output,:setPoint,:absError,]
		)
	end
	input::DanaReal
	output::DanaReal
	setPoint::DanaReal
	absError::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export IAE
function setEquationFlow(in::IAE)
	addEquation(1)
	addEquation(2)
end
function initial(in::IAE)
	addEquation(1)
end
function atributes(in::IAE,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/IAE"
	fields[:Brief]="Model IAE."
	fields[:Info]="== Inputs ==
* One input signal.
* One setpoint signal.

== Outputs ==
* One output signal.
"
	drive!(fields,_)
	return fields
end
IAE(_::Dict{Symbol,Any})=begin
	newModel=IAE()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using VariationalInequality
using JuMP

using Base.Test

info("Nagurney Model")
include(joinpath(dirname(dirname(@__FILE__)), "example", "nagurney.jl"))
info("Traffic Equilibrium")
include(joinpath(dirname(dirname(@__FILE__)), "example", "traffic.jl"))
info("Fukushima1")
include(joinpath(dirname(dirname(@__FILE__)), "example", "fukushima1.jl"))


text:
mutable struct Episode{E<:AbstractEnvironment,P<:AbstractPolicy,F<:AbstractFloat}
  env::E
  policy::P
  total_reward::F # total reward of the episode
  last_reward::F
  niter::Int      # current step in this episode
  freq::Int       # number of steps between choosing actions
  maxn::Int       # max steps in an episode - should be constant during an episode
end

Episode(env::AbstractEnvironment, π::AbstractPolicy; freq = 1, maxn = maxsteps(env)) =
  Episode(env, π, 0.0, 0.0, 1, freq, maxn)

function _start!(ep::Episode{E,P,F}) where {E,P,F}
  reset!(ep.env)
  reset!(ep.policy)
  ep.total_reward = zero(F)
  ep.niter = 1
end

_done(ep::Episode, i) =
  (ep.maxn != 0 && ep.niter >= ep.maxn) || finished(ep.env, state(ep.env))

# take one step in the enviroment after querying the policy for an action
function Base.iterate(ep::Episode, i = _start!(ep))
  _done(ep::Episode, i) && return nothing

  env = ep.env
  π = ep.policy
  s = state(env)
  A = actions(env, s)  # action space
  r = reward(env)
  a = action(π, r, s, A)

  @assert(a ∈ A, "action $a is not in $A")

  # take freq steps using action a
  last_reward = 0.0
  s′ = s
  for _ ∈ 1:ep.freq
    r, s′ = step!(env, s′, a)
    last_reward += r
    _done(ep, ep.niter) && break
  end

  ep.total_reward += last_reward
  ep.last_reward = last_reward
  ep.niter = i

  (s, a, r, s′), i + 1
end

"""
    run_episode(f, env, policy)
    run_episode(env, policy) do (s, a, r, s′)
      # render or something else
    end

Helper function for running an episode,
and return the total reward gained in this episode.
"""
function run_episode(f::Base.Callable, env::AbstractEnvironment, π::AbstractPolicy)
  ep = Episode(env, π)
  for sars in ep
    f(sars)
  end
  ep.total_reward
end

# ---------------------------------------------------------------------
# iterate through many episodes

mutable struct Episodes
  env
  kw

  # note: we have different groups of strategies depending on when they should be applied
  episode_strats    # learning strategies for each episode
  epoch_strats      # learning strategies for each complete episode
  iter_strats       # learning strategies applied at every iteration
end

function Episodes(env;
                  episode_strats = [],
                  epoch_strats = [],
                  iter_strats = [],
                  kw...)
  Episodes(
    env,
    kw,
    MetaLearner(episode_strats...),
    MetaLearner(epoch_strats...),
    MetaLearner(iter_strats...)
  )
end

length_state(eps::Episodes) = length(state(eps.env)) + length(eps.last_action)

# the main function... run episodes until stopped by one of the epoch/iter strats
function learn!(policy, eps::Episodes)
  # setup
  setup!(eps.epoch_strats, policy)
  setup!(eps.iter_strats, policy)

  # loop over epochs until done
  done = false
  epoch = 1
  iter = 1

  while !done
    # one episode
    setup!(eps.episode_strats, policy)
    ep = Episode(eps.env, policy; eps.kw...)
    for sars′ in ep
      learn!(policy, sars′...)

      # learn steps
      for metalearner in (eps.episode_strats, eps.epoch_strats, eps.iter_strats)
        for strat in metalearner.managers
          learn!(policy, strat, sars′)
        end
      end

      # iter steps
      timestep = ep.niter
      hook(eps.episode_strats, ep, timestep)
      hook(eps.epoch_strats, ep, epoch)
      hook(eps.iter_strats, ep, iter)

      # finish the timestep with checks
      if finished(eps.episode_strats, policy, timestep)
        break
      end
      if finished(eps.epoch_strats, policy, epoch) || finished(eps.iter_strats, policy, iter)
        done = true
        break
      end
      iter += 1
    end
    info("Finished episode $epoch after $(ep.niter) steps. Reward: $(ep.total_reward) mean(Reward): $(ep.total_reward/max(ep.niter,1))")
    cleanup!(eps.episode_strats, policy)
    epoch += 1

  end

  # tear down
  cleanup!(eps.epoch_strats, policy)
  cleanup!(eps.iter_strats, policy)
  return
end

# function hook(policy, ep::Episodes, i)
#     if ep.should_reset
#         reset!(ep.env)
#         reset!(policy)
#         ep.should_reset = false
#         ep.total_reward = 0.0
#         ep.nsteps = 0
#         for learner in ep.learners
#             setup!(learner, policy)
#         end
#     end
#
#     # take one step in the enviroment after querying the policy for an action
# 	env = ep.env
# 	s = state(env)
#     A = actions(env, s)
#     r = reward(env)
# 	a = action(policy, r, s, A)
#     if !(a in A)
#         warn("action $a is not in $A")
#         # a = rand(A)
#     end
#     @assert a in A
# 	r, s′ = step!(env, s, a)
# 	ep.total_reward += r
#
#     # "sars" learn step for the policy...
#     #   note: ensures that the final reward is included in the learning
#     learn!(policy, s, a, r, s′)
#
#     ep.nsteps += 1
#     for learner in ep.learners
#         learn!(policy, learner, ep.nsteps)
#         hook(learner, policy, ep.nsteps)
#     end
#
#     # if this episode is done, just flag it so we reset next time
#     if finished(env, s′) || any(learner -> finished(learner, policy, ep.nsteps), ep.learners)
#         ep.should_reset = true
#         ep.nepisode += 1
#         for learner in ep.learners
#             cleanup!(learner, policy)
#         end
#         info("Finished episode $(ep.nepisode) after $(ep.nsteps) steps. Reward: $(ep.total_reward)")
#     end
#     return
# end


text:
using ITensors,
      Test

@testset "MPS Basics" begin

  N = 10
  sites = SiteSet(N,2)
  psi = MPS(sites)
  @test length(psi) == N

  @test siteindex(psi,2) == sites[2]
  @test hasindex(psi[3],linkindex(psi,2))
  @test hasindex(psi[3],linkindex(psi,3))

  psi[1] = ITensor(sites[1])
  @test hasindex(psi[1],sites[1])

  @testset "RandomMPS" begin
    phi = randomMPS(sites)
    @test hasindex(phi[1],sites[1])
    @test norm(phi[1])≈1.0
    @test hasindex(phi[4],sites[4])
    @test norm(phi[4])≈1.0
  end

  @testset "inner different MPS" begin
    phi = randomMPS(sites)
    psi = randomMPS(sites)
    phipsi = dag(phi[1])*psi[1]
    for j = 2:N
      phipsi *= dag(phi[j])*psi[j]
    end
    @test phipsi[] ≈ inner(phi,psi)
  end

  @testset "inner same MPS" begin
    psi = randomMPS(sites)
    psidag = dag(psi)
    primelinks!(psidag)
    psipsi = psidag[1]*psi[1]
    for j = 2:N
      psipsi *= psidag[j]*psi[j]
    end
    @test psipsi[] ≈ inner(psi,psi)
  end

end


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

mutable struct PolicyAssignmentProperties <: SwaggerModel
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    policyDefinitionId::Any # spec type: Union{ Nothing, String } # spec name: policyDefinitionId
    scope::Any # spec type: Union{ Nothing, String } # spec name: scope
    notScopes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: notScopes
    parameters::Any # spec type: Union{ Nothing, ParameterValues } # spec name: parameters
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    metadata::Any # spec type: Union{ Nothing, Any } # spec name: metadata
    enforcementMode::Any # spec type: Union{ Nothing, String } # spec name: enforcementMode

    function PolicyAssignmentProperties(;displayName=nothing, policyDefinitionId=nothing, scope=nothing, notScopes=nothing, parameters=nothing, description=nothing, metadata=nothing, enforcementMode=nothing)
        o = new()
        validate_property(PolicyAssignmentProperties, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(PolicyAssignmentProperties, Symbol("policyDefinitionId"), policyDefinitionId)
        setfield!(o, Symbol("policyDefinitionId"), policyDefinitionId)
        validate_property(PolicyAssignmentProperties, Symbol("scope"), scope)
        setfield!(o, Symbol("scope"), scope)
        validate_property(PolicyAssignmentProperties, Symbol("notScopes"), notScopes)
        setfield!(o, Symbol("notScopes"), notScopes)
        validate_property(PolicyAssignmentProperties, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        validate_property(PolicyAssignmentProperties, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(PolicyAssignmentProperties, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(PolicyAssignmentProperties, Symbol("enforcementMode"), enforcementMode)
        setfield!(o, Symbol("enforcementMode"), enforcementMode)
        o
    end
end # type PolicyAssignmentProperties

const _property_map_PolicyAssignmentProperties = Dict{Symbol,Symbol}(Symbol("displayName")=>Symbol("displayName"), Symbol("policyDefinitionId")=>Symbol("policyDefinitionId"), Symbol("scope")=>Symbol("scope"), Symbol("notScopes")=>Symbol("notScopes"), Symbol("parameters")=>Symbol("parameters"), Symbol("description")=>Symbol("description"), Symbol("metadata")=>Symbol("metadata"), Symbol("enforcementMode")=>Symbol("enforcementMode"))
const _property_types_PolicyAssignmentProperties = Dict{Symbol,String}(Symbol("displayName")=>"String", Symbol("policyDefinitionId")=>"String", Symbol("scope")=>"String", Symbol("notScopes")=>"Vector{String}", Symbol("parameters")=>"ParameterValues", Symbol("description")=>"String", Symbol("metadata")=>"Any", Symbol("enforcementMode")=>"String")
Base.propertynames(::Type{ PolicyAssignmentProperties }) = collect(keys(_property_map_PolicyAssignmentProperties))
Swagger.property_type(::Type{ PolicyAssignmentProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PolicyAssignmentProperties[name]))}
Swagger.field_name(::Type{ PolicyAssignmentProperties }, property_name::Symbol) =  _property_map_PolicyAssignmentProperties[property_name]

const _allowed_PolicyAssignmentProperties_enforcementMode = ["Default", "DoNotEnforce"]

function check_required(o::PolicyAssignmentProperties)
    true
end

function validate_property(::Type{ PolicyAssignmentProperties }, name::Symbol, val)
    if name === Symbol("enforcementMode")
        Swagger.validate_param(name, "PolicyAssignmentProperties", :enum, val, _allowed_PolicyAssignmentProperties_enforcementMode)
    end
end


text:
using MarketData

facts("find methods") do

  context("findall returns correct row numbers array") do
    @fact cl[findall(cl .> op)].timestamp[1] => Date(2000,1,3)
    @fact length(findall(cl .> op))          => 244
  end

  context("findwhen returns correct Dates array") do
    @fact findwhen(cl .> op)[1]      => Date(2000,1,3)
    @fact length(findwhen(cl .> op)) => 244
  end
end

facts("split date operations") do

  context("from and to correctly subset") do
    @fact length(from(cl, 2001,12,28)) => 2
    @fact length(to(cl, 2000,1,4))     => 2
  end

  context("bydate methods correctly subset") do
     @fact by(cl,2001, period=year).timestamp[1]   => Date(2001,1,2)
     @fact by(cl,2, period=month).timestamp[1]     => Date(2000,2,1)
     @fact by(cl,4, period=day).timestamp[1]       => Date(2000,1,4)
     @fact by(cl,5, period=dayofweek).timestamp[1] => Date(2000,1,7)
     @fact by(cl,4, period=dayofyear).timestamp[1] => Date(2000,1,4)
  end
end

facts("element wrappers") do

  context("type element wrappers isolate elements") do
     @fact isa(timestamp(cl), Array{Date,1}) => true
     @fact isa(values(cl), Array{Float64,1})              => true
     @fact isa(values(ohlc), Array{Float64,2})            => true
     @fact isa(colnames(cl), Array{UTF8String, 1})       => true
  end
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct ArrayFeatureExtractor <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ArrayFeatureExtractor(; kwargs...)
        obj = new(meta(ArrayFeatureExtractor), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ArrayFeatureExtractor
const __meta_ArrayFeatureExtractor = Ref{ProtoMeta}()
function meta(::Type{ArrayFeatureExtractor})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ArrayFeatureExtractor)
            __meta_ArrayFeatureExtractor[] = target = ProtoMeta(ArrayFeatureExtractor)
            pack = Symbol[:extractIndex]
            allflds = Pair{Symbol,Union{Type,String}}[:extractIndex => Base.Vector{UInt64}]
            meta(target, ArrayFeatureExtractor, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ArrayFeatureExtractor[]
    end
end
function Base.getproperty(obj::ArrayFeatureExtractor, name::Symbol)
    if name === :extractIndex
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{UInt64}
    else
        getfield(obj, name)
    end
end

export ArrayFeatureExtractor


text:
module CryptoGroups

include("utils.jl")

function spec end
function specialize end
function order end
function bitlength end
function <| end
function tobits end

function validate end # isvalid
function value end
function modulus end

function cofactor end

function a end
function b end

# In most cases probably using sizeof directly is more appropriate
bitlength(::Type{T}) where T <: Integer = sizeof(T) * 8

function bitlength(p::Integer)

    bits = bitstring(p)
    start = findfirst(x -> x == '1', bits)
    N = length(bits) - start + 1

    return N
end

function bitlength(p::BigInt)

    bytes = int2bytes(p)
    bits = bitstring(bytes[end])
    start = findfirst(x -> x == '1', bits)
    N = length(bytes) * 8  - (start - 1)

    return N
end

include("Fields/Fields.jl")
include("Curves/Curves.jl")
include("Specs/Specs.jl")

include("groups.jl")
include("elgamal.jl")
include("spec.jl")

export @bin_str

export spec, specialize, order, bitlength

import .Fields: Field, BinaryField, FP, F2GNB, F2PB, PrimeField
export Field, FP, F2GNB, F2PB # frombits and tobits method shall be dealt with convert

import .Curves: AbstractPoint, ECPoint, AffinePoint, BinaryCurve, validate, gx, gy, a, b, oncurve
export AbstractPoint, ECPoint, AffinePoint, BinaryCurve, validate #, gx, gy, a, b, oncurve

import .Specs: MODP, Koblitz, ECP, EC2N, PB, GNB, generator, Hash, PRG, RO, ROPRG
export MODP, Koblitz, ECP, EC2N, PB, GNB, generator, Hash, PRG, RO, ROPRG

export ElGamal

end


text:
@testset "SimpleCodec" begin
    sc = SimpleCodec(x -> x .^ 2, x -> sqrt.(x))
    @test encode(sc, 3) == 9
    @test encode(sc, [1, 2]) == [1, 4]
    @test decode(sc, 100) == 10
    @test decode(sc, [1, 9]) == [1, 3]
end


text:
@testset "v0.4 deprecations" begin
    x_inv_microns = [10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.6, 4.0, 3.0, 2.78, 2.27, 1.82, 1.43, 1.11, 0.8, 0.63, 0.46]
    wave = 1e4 ./ x_inv_microns
    for law in [ccm89, od94]
        @test_deprecated law(wave)
    end
end

@testset "v0.7 deprecations" begin
    wave = 3000
    flux = 1

    @test_deprecated redden(wave, flux, 0.3)
    @test_deprecated deredden(wave, flux, 0.3)

    for law in [ccm89, od94, cal00, vcg04, gcc09]
        @test_deprecated law(wave)
    end
end


text:
#!/home/s1/bin/julia
#filesearch.jl
using JSON, DocOpt, Random, DataFrames, CSV, Match

export main
doc="""fileSearch.

Usage:
  fileSearch.jl [--path <path>] [--filter <filter>] [--random <random>] [--save <savefile>] [--pics <picsdir>] [--last <last>]
  fileSearch.jl --help
  fileSearch.jl --version

Options:
  -H --help                Show this screen
  -V --version             Show version
  -P --path=<path>         Path for search [default: ./ ]
  -F --filter=<filter>     Regex to filter
  -R --random=<random>     Randomly pick files from search
  -S --save=<savefile>     File to save search results
  --pics=<picsdir>      Create a pictures dir
  -L --last=<last>        Return last <count> dated files
  Note: If --last and --random are both selected, <random> files are selected from the <last> files
"""

#kwargs
Base.@kwdef mutable struct Options
  path::String=""
  filter=r".*\.mp4"
  random::Union{Int,Nothing}=nothing
  save=false
  pics::Union{String,Nothing}=nothing
  last::Union{Int,Nothing}=nothing
end

"Set values of Option to match command line flags"
function setup!(flag::String, value::Any,  opt::Options)
  @match (flag, value)  begin
  ("--path"   || "-P",  p::String) => (opt.path=p;println("Path set to $(opt.path)"))
  ("--filter" || "-F",  f::String) => (opt.filter=Regex(f))
  ("--random" || "-R",  r::String) => (opt.random=parse(Int64, r);println("Picking $r random files"))
  ("--save"   || "-S",  s::String) => (opt.save=s)
  ("--pics" ,           p::String) => (opt.pics=p)
  ("--last"   || "-L",  l::String) => (opt.last=parse(Int64, l);println("Picking last $l files"))
#  bad                     => println("Unknown argument: $bad")
  end
end

"Gets all the file names in a directory tree"
function getAllFiles(path;ignore_hidden=true)
  root, dirs, files = first(walkdir(path;follow_symlinks=true))
  allFiles = map(f -> joinpath(path, f), files)
  for dir in dirs
    if ignore_hidden && dir[1] == '.' continue end
    dirFiles=getAllFiles(joinpath(root, dir))
    append!(allFiles, dirFiles)
  end
  allFiles
end

function write(files, options)
  shFile=joinpath(options.path, "vlc.sh")
  count=length(files)
  println("Writing $count files to $shFile")
  open(shFile, "w") do io
    println(io, "#!/bin/bash")
    for fname in files
      pFile="file:///" * "'" * fname * "'"
      rs="vlc $pFile &"
      println(io, rs)
    end
  end
  chmod(shFile, 0o777)
end

function getRandom(files, options)
  randomFiles=[]
  for r in 1:options.random
    r=rand(1:size(files)[1])
#    println("Random file: $r")
    push!(randomFiles, files[r])
  end
  randomFiles
end

getFileCTime(file)=(ctime(stat(file)), file)
getFileTimes(allFiles)=map(getFileCTime, allFiles)

"Get the last created files by timestamp"
function getLast(files, options)
  fileTimes=map(getFileCTime, files)
  sortFiles=sort(fileTimes, rev=true)
  lastFiles=[sortFiles[r][2] for r in 1:options.last]
  #println(lastFiles)
  lastFiles
end

is(something)=!isnothing(something)

function main()
  options=Options()
  args=docopt(doc)
  [setup!(arg.first, arg.second, options) for arg in args]
  println("Options: $options")
  allFiles=getAllFiles(options.path)
  isMatch(str)=!isnothing(match(options.filter, str))
  filterFiles=filter(isMatch, allFiles)
  fileLength=size(filterFiles)[1]
  println("There are $(fileLength) files in the directory tree. ")
  if fileLength == 0
    println("No files found in path. Exiting")
    exit()
  end

  writeFiles=filterFiles
  if is(options.last)
    writeFiles=getLast(filterFiles, options)
    println("Selected last $(length(writeFiles)) files")
    if is(options.random)
      writeFiles=getRandom(writeFiles, options)
      println("Selected last random $(length(writeFiles)) file")
    end
  elseif is(options.random)
    writeFiles=getRandom(filterFiles, options)
  end
  if (is(options.random) || is(options.last)) write(writeFiles, options)  end
#  if !isnothing(options.save) CSV.write(options.save, allFiles) end
end

main()
#args=docopt(doc)
#println(args)


text:
const TEXTURE_PARAMS = @map Symbol, GLenum {
    depth_mode   = GL_DEPTH_STENCIL_TEXTURE_MODE,
    base_level   = GL_TEXTURE_BASE_LEVEL,
    compare_func = GL_TEXTURE_COMPARE_FUNC,
    compare_mode = GL_TEXTURE_COMPARE_MODE,
    lod_bias     = GL_TEXTURE_LOD_BIAS,
    min_filter   = GL_TEXTURE_MIN_FILTER,
    mag_filter   = GL_TEXTURE_MAG_FILTER,
    min_lod      = GL_TEXTURE_MIN_LOD,
    max_lod      = GL_TEXTURE_MAX_LOD,
    max_level    = GL_TEXTURE_MAX_LEVEL,
    swizzle_r    = GL_TEXTURE_SWIZZLE_R,
    swizzle_g    = GL_TEXTURE_SWIZZLE_G,
    swizzle_b    = GL_TEXTURE_SWIZZLE_B,
    swizzle_a    = GL_TEXTURE_SWIZZLE_A,
    wrap_s       = GL_TEXTURE_WRAP_S,
    wrap_t       = GL_TEXTURE_WRAP_T,
    wrap_r       = GL_TEXTURE_WRAP_R,
}

function texture2D(pixels::Union{AbstractArray, Ptr}, width::Integer, height::Integer)
    texture, type = GLuint'[], g"TEXTURE 2D"
    glGenTextures(1, texture)
    glBindTexture(type, texture')
    glTexImage2D(type, 0, g"RGB", width, height, 0, g"RGB", g"UNSIGNED BYTE", pixels)
    glGenerateMipmap(g"texture 2d")
    return texture'
end

function texture2D(pixels::Union{AbstractArray, Ptr}, width::Integer, height::Integer, params::Dict{Symbol, GLenum})
    texture = texture2D(pixels, width, height)
    for (param, value) in params
        param ∉ keys(TEXTURE_PARAMS) && error("Unkown texture parameter '$param'")
        glTexParameteri(g"TEXTURE 2D", TEXTURE_PARAMS[param], value)
    end
    return texture
end


text:
using Base: @propagate_inbounds
using CUDA
using Adapt
using OffsetArrays
using Statistics

using Oceananigans.Architectures
using Oceananigans.Utils
using Oceananigans.Grids: interior_indices, interior_parent_indices

import Base: minimum, maximum, extrema
import Statistics: mean
import Oceananigans: location, instantiated_location
import Oceananigans.Architectures: architecture
import Oceananigans.Grids: interior_x_indices, interior_y_indices, interior_z_indices
import Oceananigans.Grids: total_size, topology, nodes, xnodes, ynodes, znodes, xnode, ynode, znode
import Oceananigans.Utils: datatuple

const ArchOrNothing = Union{AbstractArchitecture, Nothing}
const GridOrNothing = Union{AbstractGrid, Nothing}

"""
    AbstractField{LX, LY, LZ, G, T, N}

Abstract supertype for fields located at `(LX, LY, LZ)`
and defined on a grid `G` with eltype `T` and `N` dimensions.

Note: we need the parameter `T` to subtype AbstractArray.
"""
abstract type AbstractField{LX, LY, LZ, G <: GridOrNothing, T, N} <: AbstractArray{T, N} end

Base.IndexStyle(::AbstractField) = IndexCartesian()

# Endpoint for recursive `datatuple` function:
data(a) = a
@inline datatuple(obj::AbstractField) = data(obj)

#####
##### AbstractField functionality
#####

"Returns the location `(LX, LY, LZ)` of an `AbstractField{LX, LY, LZ}`."
@inline location(a) = (Nothing, Nothing, Nothing) # used in AbstractOperations for location inference
@inline location(::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = (LX, LY, LZ) # note no instantiation
@inline location(f::AbstractField, i) = location(f)[i]
@inline instantiated_location(::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = (LX(), LY(), LZ())

"Returns the architecture of on which `f` is defined."
architecture(f::AbstractField) = architecture(f.grid)

"Returns the topology of a fields' `grid`."
@inline topology(f::AbstractField, args...) = topology(f.grid, args...)

"""
    size(f::AbstractField)

Returns the size of an `AbstractField{LX, LY, LZ}` located at `LX, LY, LZ`.
This is a 3-tuple of integers corresponding to the number of interior nodes
of `f` along `x, y, z`.
"""
Base.size(f::AbstractField) = size(location(f), f.grid)
Base.length(f::AbstractField) = prod(size(f))
Base.parent(f::AbstractField) = f

"""
    total_size(field::AbstractField)

Returns a 3-tuple that gives the "total" size of a field including
both interior points and halo points.
"""
total_size(f::AbstractField) = total_size(location(f), f.grid)

interior(f::AbstractField) = f

#####
##### Coordinates of fields
#####

@propagate_inbounds xnode(i, ψ::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = xnode(LX(), i, ψ.grid)
@propagate_inbounds ynode(j, ψ::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = ynode(LY(), j, ψ.grid)
@propagate_inbounds znode(k, ψ::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = znode(LZ(), k, ψ.grid)

@propagate_inbounds xnode(i, j, k, ψ::AbstractField) = xnode(instantiated_location(ψ)..., i, j, k, ψ.grid)
@propagate_inbounds ynode(i, j, k, ψ::AbstractField) = ynode(instantiated_location(ψ)..., i, j, k, ψ.grid)
@propagate_inbounds znode(i, j, k, ψ::AbstractField) = znode(instantiated_location(ψ)..., i, j, k, ψ.grid)

xnodes(ψ::AbstractField) = xnodes(location(ψ, 1), ψ.grid)
ynodes(ψ::AbstractField) = ynodes(location(ψ, 2), ψ.grid)
znodes(ψ::AbstractField) = znodes(location(ψ, 3), ψ.grid)

nodes(ψ::AbstractField; kwargs...) = nodes(location(ψ), ψ.grid; kwargs...)

#####
##### Some conveniences
#####

for f in (:+, :-)
    @eval Base.$f(ϕ::AbstractArray, ψ::AbstractField) = $f(ϕ, interior(ψ))
    @eval Base.$f(ϕ::AbstractField, ψ::AbstractArray) = $f(interior(ϕ), ψ)
end

function Statistics.norm(a::AbstractField)
    r = zeros(a.grid, 1)
    Base.mapreducedim!(x -> x * x, +, r, a)
    return CUDA.@allowscalar sqrt(r[1])
end


text:
# Program to Swap two numbers with temp

print("Enter first value a : ")
a = parse(Int64,readline())
print("Enter second value b: ")
b = parse(Int64,readline())

print("\n")
print("Before swapping : a = $a and b = $b")

temp = a
a = b
b = temp

print("\n")
print("After swapping : a = $a and b = $b")

max_stars_repo_name: cashewcodes/Basic-Julia-Projects

max_issues_repo_name: cashewcodes/Basic-Julia-Projects

max_forks_repo_name: cashewcodes/Basic-Julia-Projects


text:
using SparseArrays
using BlockArrays
using IterTools

struct DampedWaveNet
	incidence_matrix    ::SparseMatrixCSC{Int8,Int64}
	edge_parameters     ::NamedTuple{(:a,:b,:d,:l,:n),
	                           Tuple{Vector{Float64},
	                                 Vector{Float64},
	                                 Vector{Float64},
	                                 Vector{Float64},
	                                 Vector{Int64  }}}
	boundary_conditions ::Vector{Char}

	function DampedWaveNet(inc_mat,ep,bc)
		inc_mat, ep, bc = convert(Tuple{fieldtypes(DampedWaveNet)...},(inc_mat,ep,bc))

		@assert all(in.(nonzeros(inc_mat),Ref([-1,1]))) "" *
		"Invalid incidence matrix: found value(s) other than {-1,0,1}"

		@assert all(nnz.(eachcol(inc_mat)) .== 2) "" *
		"Invalid incidence matrix: found column(s) with other than 2 entries"

		@assert all(sum(abs.(inc_mat),dims=2) .> 0) "" *
		"Invalid incidence matrix: found disconnected vertices"

		@assert sum(sum(abs.(inc_mat),dims=2) .== 1) > 0 "" *
		"Invalid incidence matrix: found no boundary vertices"

		@assert all(length.(values(ep)) .== size(inc_mat)[2]) "" *
		"Invalid edge parameters: need same number as edges"

		@assert all(ep.d .>= 0) "" *
		"Invalid edge parameters: found damping coefficient(s) < 0"

		@assert all(ep.l .> 0) "" *
		"Invalid edge parameters: found length(s) <= 0"

		@assert all(ep.n .> 0) "" *
		"Invalid edge parameters: found cell number(s) < 1"

		@assert length(bc) == sum(sum(abs.(inc_mat),dims=2) .== 1) "" *
		"Invalid boundary conditions: need same number as boundary nodes"

		@assert all(in.(bc,Ref(['p','m']))) "" *
		"Invalid boundary conditions: found identifier other than {p,m}"

		return new(inc_mat,ep,bc)
	end
end

function DampedWaveNet(id::String)
	if      id == "pipe"
		inc_mat = reshape([ 1;-1],:,1)
		ep = (a=[1],b=[1],d=[1],l=[1],n=[10])
		bc = ['p','p']

	elseif id == "fork"
		inc_mat = [-1  0  0;
		       0  1  0;
		       0  0  1;
		       1 -1 -1]
		ep = (a=[1  ,1  ,1 ],
		      b=[1  ,1  ,1 ],
		      d=[1  ,1  ,1 ],
		      l=[2  ,1  ,10],
		      n=[40 ,30 ,90])
		bc = ['p','p','p']

	elseif id == "diamond"
		inc_mat = [-1  0  0  0  0  0  0;
		       0  0  0  0  0  0  1;
		       1 -1 -1  0  0  0  0;
		       0  1  0 -1 -1  0  0;
		       0  0  1  1  0 -1  0;
		       0  0  0  0  1  1 -1]
		ep = (a=[4  ,4  ,1  ,1  ,1  ,4  ,4  ]     ,
		      b=[1/4,1/4,1  ,1  ,1  ,1/4,1/4]     ,
		      d=[1  ,1  ,8  ,8  ,8  ,1  ,1  ]./80 ,
		      l=[1  ,1  ,1  ,1  ,1  ,1  ,1  ]     ,
		      n=[1  ,1  ,1  ,1  ,1  ,1  ,1  ].*500)
		bc = ['p','p']

	else
		throw("Config id \'"*id*"\' not recognized")
	end

	return DampedWaveNet(inc_mat,ep,bc)
end

function build(problem::DampedWaveNet)
	#Convenience
	inc_mat = sparse(problem.incidence_matrix')
	ep = problem.edge_parameters
	bc = problem.boundary_conditions

	#Index calculations
	n_p  = sum(ep.n)                 #Number of pressure variables
	n_m  = sum(ep.n .+ 1)            #Number of mass flow variables
	n_b  = length(bc)                #Number of boundary conditions
	n_ip = sum(nnz.(eachcol(inc_mat)).-1) #Number of internal conditions for p
	n_im = size(im)[2] - n_b         #Number of internal conditions for m
	n_x  = n_p+n_m+n_ip+n_im+n_b     #Number of state variables
	n_u  = n_b                       #Number of input variables
	n_y  = n_b                       #Number of output variables

	i_ep = collect(eachblock(BlockArray(1:n_p,ep.n   ))) #Edge indices for p
	i_em = collect(eachblock(BlockArray(1:n_m,ep.n.+1))) #Edge indices for m

	#Global FEM system
	E   = spzeros(n_x,n_x)
	A   = spzeros(n_x,n_x)
	B   = spzeros(n_x,n_u)

	Ew  = PseudoBlockArray(E,[n_p,n_m,n_ip+n_im+n_b],[n_p,n_m,n_ip+n_im+n_b])
	Aw  = PseudoBlockArray(A,[n_p,n_m,n_ip,n_im,n_b],[n_p,n_m,n_ip,n_im,n_b])
	Bw  = PseudoBlockArray(B,[n_p+n_m+n_ip+n_im,n_b],[n_b])

	Mp  = view(Ew,Block(1,1)) #Mass matrix for p
	Mm  = view(Ew,Block(2,2)) #Mass matrix for m
	Gp  = view(Aw,Block(2,1)) #Negative gradient matrix for p
	Gm  = view(Aw,Block(1,2)) #Negative gradient matrix for m
	Dm  = view(Aw,Block(2,2)) #Damping matrix for m
	Cip = view(Aw,Block(3,1)) #Internal condition matrix for p
	Cim = view(Aw,Block(4,2)) #Internal condition matrix for m
	Cbp = view(Aw,Block(5,1)) #Boundary condition matrix for p
	Cbm = view(Aw,Block(5,2)) #Boundary condition matrix for m
	Bu  = view(Bw,Block(2,1)) #Input matrix

	#Edge physics
	for (e,h) in enumerate(ep.l./ep.n)
		#Local FEM system on edge e
		Mp_loc = [ 1]*h*ep.a[e]
		Mm_loc = [ 2  1; 1  2]*h/6*ep.b[e]
		Gp_loc = [-1; 1]
		Gm_loc = [ 1 -1]
		Dm_loc = [-2 -1;-1 -2]*h/6*ep.d[e]

		#Contributions for each element
		for (i_p,i_m) in zip(eachrow(i_ep[e]),collect.(partition(i_em[e],2,1)))
			Mp[i_p,i_p] += Mp_loc
			Mm[i_m,i_m] += Mm_loc
			Gp[i_m,i_p] += Gp_loc
			Gm[i_p,i_m] += Gm_loc
			Dm[i_m,i_m] += Dm_loc
		end
	end

	#Algebraic conditions
	i_ip, i_im, i_b = (1,1,1) #Counter for each type of condition

	for v in eachcol(inc_mat)
		#Indices, directions of edges connected to v, corresponding variable indices
		es, ds  = (rowvals(v),nonzeros(v))
		js(i_e) = ifelse.(ds .< 0,first.(i_e[es]),last.(i_e[es]))

		if length(es) > 1
			#Internal conditions for p
			for j_pair in collect.(partition(js(i_ep),2,1))
				Cip[i_ip,j_pair] = [1,-1]
				i_ip += 1
			end

			#Internal conditions for m
			Cim[i_im,js(i_em)] .= ds
			i_im += 1
		else
			#Boundary conditions
			Cb, j       = (bc[i_b] == 'p') ? (Cbp,js(i_ep)) : (Cbm,js(i_em))
			Cb[i_b,j  ] = ds
			Bu[i_b,i_b] = 1
			i_b += 1
		end
	end

	#Apply Lagrange multiplier: L .= C'
	view(Aw,Block(1,3)) .= Cip'
	view(Aw,Block(2,4)) .= Cim'
	view(Aw,Block(1,5)) .= Cbp'
	view(Aw,Block(2,5)) .= Cbm'

	#Packaging
	return (shape=(n_x,n_u,n_y),E=E,A=A,B=B)
end


text:
#tools.jl - a few tools which make unum work much smoother...  To be reassigned
#in the directory structure

doc"""
the `@gen_code` macro rejigs the standard julia `@generate` macro so that at the
end the function expects a `code` expression variable that can be created and
automatically extended using the `@code` macro.
"""
macro gen_code(f)
  #make sure this macro precedes a function definition.
  isa(f, Expr) || error("gen_code macro must precede a function definition")
  (f.head == :function) || error("gen_code macro must precede a function definition")

  #automatically generate a 'code-creation' statement at the head of the function.
  unshift!(f.args[2].args, :(code = :(nothing)))
  #insert the code release statement at the tail of the function.
  push!(f.args[2].args, :(code))

  #return the escaped function to the parser so that it generates the new function.
  ##next, wrap the function f inside of the @generated macro and escape it
  esc(:(@generated $f))
end

doc"""
  `@unumbers` imports a series of useful int64 and int16 constants:
 * o64:  one(UInt64)
 * z64:  zero(UInt64)
 * t64:  0x8000_0000_0000_0000
 * f64:  0xFFFF_FFFF_FFFF_FFFF
 * o16:  one(UInt16)
 * z16:  zero(UInt16)
"""
macro unumbers()
  :(import Unums: o64, t64, z64, f64, o16, z16)
end
export @unumbers

#a helper function that makes the parsing process much cleaner.
function __append_code(a, b)
  return :($a; $b)
end

doc"""
the `@code` macro is used in conjunction to put a line of code (or a quote) onto
the growing `code` variable that is expected by the `@gen_code` macro.
"""
macro code(e)
  isa(e, Expr) || error("can't generate code from non-expressions")
  #since code variable has already been created, instead append it to the previous code statement.
  return esc(:(code = __append_code(code, $e)))
end

#fname extracts the function name from the expression
function fname(ex::Expr)
  ex.args[1].args[1]
end
#vfunc generates a type-parameter variadic function head.
function vfunc(fn)
  :($fn{ESS,FSS})
end

if options[:devmode]
  to16(n::Int64) = UInt16(n)
else
  to16(n::Int64) = UInt16(reinterpret(UInt64, n) & 0x0000_0000_0000_FFFF)
end
doc"""
  `Unums.to16(::Int64)` converts a signed Int64 value to an unsigned Int16.
  Normally one would do this by doing a direct conversion using call(Int16, ...)
  however, this potentially throws inexact errors, so we will craftily hint to
  LLVM that we are not worried about this possibility, except in dev mode.
"""
to16

doc"""
  the `@universal` macro is prepended to a function defined with parameters that
  are generic Unum types (Unum, Ubound, Gnum), and generates two functions,
  with each of the \*Small{ESS,FSS} and \*Large{ESS,FSS} suffixes.
"""
#creates a universal function f that operates across all types of unums
macro universal(f)
  if (f.head == :(=))
    (f.args[1].head == :call) || throw(ArgumentError("@universal macro must operate on a function"))
  elseif (f.head == :function)
    nothing  #we're good.
  else
    throw(ArgumentError("@universal macro must operate on a function"))
  end

  #copy f into two function definitions, fsmall and flarge
  fsmall = copy(f)
  flarge = copy(f)

  #extract the functionname and append the {ESS,FSS} signature onto the functionname
  functionname = fname(f)
  functioncall = vfunc(functionname)

  fsmall.args[1].args[1] = functioncall
  flarge.args[1].args[1] = functioncall

  #next work with the parameters
  parameters = f.args[1].args

  stypedefs = quote
    G = GnumSmall{ESS,FSS}
    B = UboundSmall{ESS,FSS}
    U = UnumSmall{ESS,FSS}
  end

  ltypedefs = quote
    G = GnumLarge{ESS,FSS}
    B = UboundLarge{ESS,FSS}
    U = UnumLarge{ESS,FSS}
  end

  #append these type definitions onto fsmall and flarge.
  unshift!(fsmall.args[2].args, stypedefs)
  unshift!(flarge.args[2].args, ltypedefs)

  for idx = 2:length(parameters)
    if (isa(parameters[idx], Expr)
         && (parameters[idx].head == :(::)))
      utypeparam = parameters[idx].args[2]
      if (utypeparam in [:Unum, :Ubound, :Gnum])
        stype = Symbol(utypeparam, :Small)
        ltype = Symbol(utypeparam, :Large)
        fsmall.args[1].args[idx].args[2] = :($stype{ESS,FSS})
        flarge.args[1].args[idx].args[2] = :($ltype{ESS,FSS})
      elseif isa(utypeparam, Expr)
        if (utypeparam.head == :curly) && (utypeparam.args[1] == :Type) && (utypeparam.args[2] in [:Unum, :Ubound, :Gnum])
          stype = Symbol(utypeparam.args[2], :Small)
          ltype = Symbol(utypeparam.args[2], :Large)
          fsmall.args[1].args[idx].args[2].args[2] = :($stype{ESS,FSS})
          flarge.args[1].args[idx].args[2].args[2] = :($ltype{ESS,FSS})
        end
      end
    end
  end

  return esc(quote
    Base.@__doc__ $fsmall
    $flarge
  end)
end

################################################################################
# unsigned integer shorthands

#sixteen bit numbers
const z16 = zero(UInt16)
const o16 = one(UInt16)
const f16 = UInt16(0xFFFF)

#64 bit numbers
const z64 = zero(UInt64)
const o64 = one(UInt64)
const t64 = 0x8000_0000_0000_0000               #top bit
const f64 = 0xFFFF_FFFF_FFFF_FFFF               #full bits

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
module TTCrystal

using FileIO, Distributed, DistributedArrays, LinearAlgebra, SharedArrays, StaticArrays, Printf, OffsetArrays, Elliptic, Roots, NLSolversBase, Optim, Plots, Quadmath, GenericSchur, QuadGK, JLD2, FFTW, RecursiveArrayTools
using LoopVectorization

using CustomQuad

import InteractiveUtils: subtypes
using LineSearches
#import LineSearches: MoreThuente
import Base.Iterators: product
import Base: ==, <, >, <=, >=

include("parameters.jl")
include("repulsion.jl")
include("matrices.jl")
include("structure_generation.jl")
include("obj_functions.jl")
include("optimizers.jl")
include("analytic_continuation.jl")

export ParamsB1, ParamsB1_pinned, ParamsNoB1, ParamsNoB1_pinned, ReducedDispersion, Dispersion, Saver
export get_u₀, get_psamples, widen
export generate_shared_cash, generate_psamples, get_psamples, get_psamples_old, get_psamples_new, separate_psamples
export construct_objective, process_bs
export G_Cash, G_Cash2
export compute_single_period!, compute_full_span!
export precompute_step!, precompute_step, compute_grad_components!
export isB1, output, save_data, log_data, finalize!, finalize
export B, B′₀, B′₁, B′′₀₀, B′′₁₁, B′′₀₁
export bfgsObjFunc
export @unpack

end # module


text:
# ------------------------------------------------
# import
# ------------------------------------------------
using Random
using Printf
#using Distributions

# ------------------------------------------------
# main
# ------------------------------------------------
function main()
    # -------------------------------------
    # file
    # -------------------------------------
    f_train = "train.txt"
    f_test  = "test.txt"

    f_re_train = "corr_history_train"
    f_re_test  = "corr_history_test"

    # -------------------------------------
    # input & output
    # -------------------------------------
    in_n       = 10           # +1 for b
    fout_n     = 2            # final output
    out_n_at_m = [in_n, 6, fout_n]          # number of cell at middle layers

    # -------------------------------------
    # nn para
    # -------------------------------------
    num_mlayers = 1   # middle layer
    mu = 0.9
    ep = 0.5

    # -------------------------------------
    # init weight
    # -------------------------------------
    mean  = 0.0    # average
    stdev = 1.0    # variance value

    # -------------------------------------
    # nn loop
    # -------------------------------------
    main_ln  = 1000
    online_ln = 10

    # -------------------------------------
    # cal
    # -------------------------------------
    for i in 1:length(out_n_at_m)
        out_n_at_m[i] += 1
    end
    out_n = maximum(out_n_at_m)
    num_layers = num_mlayers + 2

    # ---------------------------------------------------------------------------------------
    # -------------------------------------
    # setup w
    # -------------------------------------
    w_old  = ones(num_layers-1, out_n, out_n)
    w_new  = ones(num_layers-1, out_n, out_n)
    w_temp = ones(num_layers-1, out_n, out_n)

    # dis of gauss
    for i in 1:num_layers-1
        for j in 2:out_n
            for k in 1:out_n
                #rng = MersenneTwister(1234)
                w_old[i,j,k] = rand_normal(mean, stdev)
                w_new[i,j,k] = rand_normal(mean, stdev)
            end
        end
    end

    # -------------------------------------
    # setup for nn
    # -------------------------------------
    u     = zeros(num_layers  , out_n)   # sum result
    z     = zeros(num_layers  , out_n)   # output
    temp  = zeros(out_n)                 # output
    delta = zeros(num_layers  , out_n)
    dEdw  = zeros(num_layers-1, out_n, out_n)

    roca_test  = zeros(main_ln)
    roca_train = zeros(main_ln)

    # -------------------------------------
    # read file
    # -------------------------------------
    train_re, train_data, train_ans = read_file(f_train)
    test_re,  test_data,  test_ans  = read_file(f_test)
    s = size(train_re)[1]

    # -------------------------------------
    # machine learning by NN
    # -------------------------------------
    for t in 1:main_ln
        for l in 1:online_ln
            n = rand(1:s)

            # store init input para
            z[1,1] = 0.0
            for j in 1:in_n-1
                z[1,j+1] = train_data[n,j]
            end

            # Forward propagation
            for i in 2:num_layers

                for k in 1:out_n_at_m[i]
                    u[i,k] = 1.0
                    for j in 1:out_n_at_m[i-1]
                        u[i,k] += w_new[i-1,j,k] * z[i-1,j]
                    end
                end
                """
                u[0] = 1.0 + train_data[i][0] * w_new[0,1] + train_data[i][1] * w_new[0,2]
                u[1] = 1.0 + train_data[i][0] * w_new[1,1] + train_data[i][1] * w_new[1,2]
                """

                # activation function
                if i != num_layers
                    # middle layers
                    #rectifier_func(z, u, i)
                    z = sigm_func(z, u, i)
                else
                    # output layer
                    z = softmax_func(z, u, i)
                end

                # set z[i,0]
                z[i,1] = 0.0
            end

            # Back propagation
            # delta at output layer
            i = num_layers - 1
            for k in 2:out_n_at_m[num_layers]
                delta[i,k] = z[i,k] - train_re[n,k-1]
            end

            # delta
            for ii in 1:num_layers-1
                i = num_layers - ii+1
                for k in 1:out_n_at_m[i]
                    #df = rectifier_func_derivation(u[i,k])
                    df = sigm_func_derivation(u[i,k])
                    for j in 1:out_n_at_m[i-1]
                        delta[i-1,k] += delta[i,j] *(w_new[i-1,j,k] * df )
                    end
                end
            end
            """
            delta[0] = 0.0
            delta[1] = z[1] - train_re[i,0]
            delta[2] = z[2] - train_re[i,1]
            """

            for ii in 1:num_layers-1
                i = num_layers - ii+1

                # slope
                for k in 1:out_n_at_m[i]
                    for j in 1:out_n_at_m[i-1]
                        dEdw[i-1,j,k] += delta[i,j] * z[i,k]
                    end
                end
                """
                dEdw[0,0] += delta[0] * 0.0
                dEdw[0,1] += delta[0] * train_data[i,0]
                dEdw[0,2] += delta[0] * train_data[i,1]
                dEdw[1,0] += delta[1] * 0.0
                dEdw[1,1] += delta[1] * train_data[i,0]
                dEdw[1,2] += delta[1] * train_data[i,1]
                """
            end

            # update w
            for i in 1:num_layers-1
                for j in 2:out_n
                    for k in 1:out_n
                        w_temp[i,j,k] = w_old[i,j,k] + mu*(w_new[i,j,k] - w_old[i,j,k]) - ep*dEdw[i,j,k] / online_ln
                    end
                end
            end
            #w_temp = w_old + mu*(w_new - w_old) - ep*dEdw / online_ln
            w_old  = copy(w_new)
            w_new  = copy(w_temp)

            # reset
            for i in 1:num_layers-1
                for j in 1:out_n
                    for k in 1:out_n
                        dEdw[i,j,k] = 0.0
                    end
                end
            end
            for i in 1:num_layers
                for k in 1:out_n
                    delta[i,k] = 0.0
                end
            end
        end

        # check rate of correct answer
        roca_test[t]  = test_nn( test_data,  test_ans, w_new, in_n, out_n, out_n_at_m, num_layers,  f_re_test)
        roca_train[t] = test_nn(train_data, train_ans, w_new, in_n, out_n, out_n_at_m, num_layers, f_re_train)
    end

    # output
    output_result(f_re_train, main_ln, roca_train)
    output_result(f_re_test, main_ln, roca_test)
end

# ------------------------------------
# read file
# ------------------------------------
function read_file(infile)
    skipnum = 1

    fff = []
    open(infile, "r") do f
        fff = read(f,String)
    end
    fff = split(fff, "\n", keepempty=false)
    num_data = length(fff) - skipnum

    temp   = split(fff[2],",")
    data   = zeros(num_data, length(temp)-2)
    result = zeros(num_data, 2)
    ans    = zeros(num_data)

    for i in 1+skipnum:length(fff)
        ii = i - 1
        fff[i] = replace(fff[i]," \r" => "")

        temp = split(fff[i],",")

        result[ii, (parse(Int,temp[2]) + 1)] = 1
        for j in 1:length(temp)-2
            data[ii,j] = parse(Float64, temp[j+2])
        end

        ans[ii] = parse(Int, temp[2])
    end

    #println(size(data))
    #println(size(result))
    #println(size(ans))

    return result, data, ans
end

# ------------------------------------
# output file
# ------------------------------------
function output_result(outf, main_ln, roca)
    open(outf,"w") do f
        for t in 1:main_ln
            a = @sprintf("%8.8e", roca[t])
            write(f, a)
            write(f, "\n")
        end
    end
end

# ------------------------------------------------
# whitening
# ------------------------------------------------
function whitening(data)
    nl = length(data)
    nd = length(data[1])

    # ave = 0
    for i in 1:nl
        ave = 0.0
        for j in 1:nd
            ave += data[i,j]
        end
        ave /= nd
        for j in 1:nd
            data[i,j] -= ave
        end
    end

    # variance value = 0
    for i in 1:nl
        sigma = 0.0
        for j in 1:nd
            sigma += data[i,j]^2
        end
        sigma /= nd
        sigma = (sigma)^0.5

        if sigma == 0.0
            sigma = 1.0
        end
        for j in 1:nd
            data[i,j] /= sigma
        end
    end

    return data
end

# ------------------------------------------------
# random on gauss distribution
# ------------------------------------------------
function rand_normal(mean, stdev)
    if stdev <= 0.0
        error("standard deviation must be positive")
    end
    u1 = rand()
    u2 = rand()
    r  = sqrt( -2.0*log(u1) )
    theta = 2.0*pi*u2
    rg = mean + stdev*r*sin(theta)
    return rg
end
# ------------------------------------------------
# logistic sigmoid function
# ------------------------------------------------
function sigm_func(z, u, i)
    ld = length(u[i])
    for j in 1:ld
        z[i,j] = 1.0 / ( 1.0 + exp(-u[i,j]) )
    end
    return z
end

# ------------------------------------------------
# Derivation of logistic sigmoid function
# ------------------------------------------------
function sigm_func_derivation(ui)
    y = 1 / ( 1 + exp(-ui) )
    y = (1-y)*y
    return y
end

# ------------------------------------------------
# hyperbolic tangent function
# ------------------------------------------------
function hyperbolic_tangent_func(z, u, i)
    ld = length(u[i])
    for j in 1:ld
        z[i,j] = ( exp(u[i,j]) - exp(-u[i,j]) ) / ( exp(u[i,j]) + exp(-u[i,j]) )
    end
    return z
end

# ------------------------------------------------
# Rectified linear function
# ------------------------------------------------
function rectifier_func(z, u, i)
    ld = length(u[i])
    for j in 1:ld
        z[i,j] = max( u[i,j], 0.0 )
    end
    return z
end

# ------------------------------------------------
# Derivation of Rectified linear function
# ------------------------------------------------
function rectifier_func_derivation(ui)
    fd = 1.0

    if ui <= 0.0
        fd = 0.0
    end
    return fd
end

# ------------------------------------------------
# linear mapping ( function )
# ------------------------------------------------
function linear_func(z, u, i)
    ld = length(u[i])
    for j in 1:ld
        z[i,j] = u[i,j]
    end
    return z
end

# ------------------------------------------------
# softmax function
# ------------------------------------------------
function softmax_func(z, u, i)
    ld = length(u[i])

    s = 0.0
    for j in 1:ld
        s += exp(u[i,j])
    end
    for j in 1:ld
        z[i,j] = exp(u[i,j]) / s
    end
    return z
end

# ------------------------------------------------
# test
# ------------------------------------------------
function test_nn(test_data, corr_data, w, in_n, out_n, out_n_at_m, num_layers, of)
    size = Int(length(test_data) / in_n)
    ans  = zeros(size)
    u    = zeros(num_layers, out_n) # sum result
    z    = zeros(num_layers, out_n) # output

    point = 0
    for n in 1:size
        # store init input para
        z[1,1] = 0.0
        for j in 1:in_n
            z[1,j+1] = test_data[n,j]
        end

        # Forward propagation
        for i in 2:num_layers

            for k in 1:out_n_at_m[i]
                u[i,k] = 1.0
                for j in 1:out_n_at_m[i-1]
                    u[i,k] += w[i-1,j,k] * z[i-1,j]
                end
            end
            """
            u[0] = 1.0 + train_data[i][0] * w_new[0,1] + train_data[i][1] * w_new[0,2]
            u[1] = 1.0 + train_data[i][0] * w_new[1,1] + train_data[i][1] * w_new[1,2]
            """

            # activation function
            if i != num_layers
                # middle layers
                z = rectifier_func(z, u, i)
            else
                # output layer
                z = softmax_func(z, u, i)
            end

            # set z[i,0]
            z[i,1] = 0.0
        end

        # if corr or not
        max_index = argmax(z[num_layers-1])
        ans[n] = max_index - 1

        if ans[n] == corr_data[n]
            point += 1
        end
        """
        print(z[num_layers-1])
        print(ans)
        raise ValueError("error!")
        """
    end

    # rate of correct answer
    roca = point / size

    """
    print("------------------------------")
    print(ans[:50])
    print(corr_data[:50])
    print(point)
    print(size)
    print(roca)
    """

    return roca
end

# ---------------------------------------------
# ---------------------------------------------
main()


text:
using Test
using Market

x = LogitDemand("Hej", 3)
y = NestedLogitDemand("Hej", 3)

estimate(x)

estimate(y)


text:
using Base.Test
using JLab
import PyPlot.plt

#Sample data
y = [0,0,4,4,4,4,6,8,10,11,10,8,6,4,2,0,0,0,-2,-4,-6,-4,-6,-4,-2,0,2,4,6,8,10,4,6,4,2,0,0,-2,-4,-4]
peakloc = [10,22,31,33]
lbnds_h = [6,21,23,32]
lbnds_m = [2,21,23,32]
rbnds_h = [16,23,32,36]
rbnds_m = [21,23,32,39]

#  plt.plot(y)

# #############################################################################
#Testing naive peak finding
println("########## Testing naive peak finding ##########")

#Testing hills
println("Testing hills ...")
pks_h = gethills(y)
@test pks_h.summits == peakloc
@test pks_h.lbnds == lbnds_h
@test pks_h.rbnds == rbnds_h

#Testing mountains
println("Testing mountains ...")
pks_m = getmountains(y)
@test pks_m.summits == peakloc
@test pks_m.lbnds == lbnds_m
@test pks_m.rbnds == rbnds_m

#Testing threshhold
println("Testing threshhold ...")
pks_h = gethills(y,thresh=2)
@test pks_h.summits == [10]
@test pks_h.lbnds == [2]
@test pks_h.rbnds == [32]
pks_h = gethills(y,thresh=4)
@test pks_h.summits == []

pks_m = getmountains(y,thresh=2)
@test pks_m.summits == [10]
@test pks_m.lbnds == [2]
@test pks_m.rbnds == [32]
pks_m = getmountains(y,thresh=4)
@test pks_m.summits == []

println("Functions \"getmountains\" and \"gethills\" passed tests.")


text:
function gencost(m::JuMP.AbstractModel, variable::JuMPArray{JuMP.VariableRef}, cost_component::Function)

    store = Array{JuMP.AbstractJuMPScalar,1}(undef,length(variable.axes[1]))

    for (ix, element) in enumerate(variable)
        store[ix] = cost_component(element)
    end

    gen_cost = sum(store)

    return gen_cost

end

function gencost(m::JuMP.AbstractModel, variable::JuMPArray{JuMP.VariableRef}, cost_component::Float64)

    gen_cost = sum(variable)*cost_component

    return gen_cost

end

function pwlgencost(m::JuMP.AbstractModel, variable::VariableRef, cost_component::Array{Tuple{Float64, Float64}})

    pwlvars = @variable(m, [i = 1:(length(cost_component)-1)], base_name = "pwl_{$(variable)}", start = 0.0, lower_bound = 0.0, upper_bound = (cost_component[i+1][1] - cost_component[i][1]))
     for (ix, pwlvar) in enumerate(pwlvars)
        c = @constraint(m, pwlvar <= cost_component[ix + 1][2])
        c = @constraint(m, pwlvar >= 0)
    end
    c = @constraint(m, variable == sum(pwlvars[ix] for (ix, pwlvar) in enumerate(pwlvars)))

    # TODO: Check for performance this syntax, the changes in GenericAffExpr might require refactoring

    gen_cost = AffExpr()

    for (ix, pwlvar) in enumerate(pwlvars)

        gen_cost = JuMP.add_to_expression!(gen_cost, (
                cost_component[ix + 1][1] * cost_component[ix + 1][2] - cost_component[ix][1] * cost_component[ix][2]
            ) / ( cost_component[ix + 1][2] - cost_component[ix][2] ) * pwlvar
        )

    end

    return gen_cost
end

function gencost(m::JuMP.AbstractModel, variable::JuMPArray{JuMP.VariableRef}, cost_component::Array{Tuple{Float64, Float64}})

    gen_cost = JuMP.AffExpr()

    for var in variable
        c = pwlgencost(m, var, cost_component)
        JuMP.add_to_expression!(gen_cost,c)
    end

    return gen_cost

end


text:
#
# Attributes
#

# TODO: fix inference when there are attributes
immutable WrittenAttribute{DS<:WriteDataspace,H5T<:H5Datatype,T}
    name::Symbol
    dataspace::DS
    datatype::H5T
    data::T
end

function WrittenAttribute{T}(f::JLDFile, name::Symbol, data::T)
    WrittenAttribute(name, WriteDataspace(f, data, objodr(data)), h5type(f, data), data)
end

immutable ReadAttribute
    name::Symbol
    dataspace::ReadDataspace
    datatype_class::UInt8
    datatype_offset::Int64
    data_offset::Int64
end

const EMPTY_READ_ATTRIBUTES = ReadAttribute[]

immutable AttributeHeader
    version::UInt8
    flags::UInt8
    name_size::UInt16
    datatype_size::UInt16
    dataspace_size::UInt16
end
define_packed(AttributeHeader)

Base.sizeof(attr::WrittenAttribute) = 8 + symbol_length(attr.name) + 1 + sizeof(attr.datatype) + sizeof(attr.dataspace) +
                                     numel(attr.dataspace) * sizeof(objodr(attr.data))

function write_attribute(io::IO, f::JLDFile, attr::WrittenAttribute, wsession::JLDWriteSession)
    namelen = symbol_length(attr.name)
    write(io, AttributeHeader(0x02, isa(attr.datatype, CommittedDatatype), namelen+1,
                              sizeof(attr.datatype), sizeof(attr.dataspace)))
    unsafe_write(io, Base.unsafe_convert(Ptr{Cchar}, attr.name), namelen)
    write(io, UInt8(0))
    write(io, attr.datatype)
    write(io, attr.dataspace)
    write_data(f, attr.data, objodr(attr.data), wsession)
end

function read_attribute(io::IO, f::JLDFile)
    ah = read(io, AttributeHeader)
    ah.version == 0x02 || throw(UnsupportedVersionException())
    committed = ah.flags == 1
    !committed && ah.flags != 0 && throw(UnsupportedFeatureException())

    name = Symbol(read(io, UInt8, ah.name_size-1))
    read(io, UInt8) == 0 || throw(InvalidDataException())

    datatype_end = position(io) + ah.datatype_size
    datatype_class, datatype_offset = read_datatype_message(io, f, committed)
    seek(io, datatype_end)

    dataspace_end = position(io) + ah.dataspace_size
    dataspace = read_dataspace_message(io)
    seek(io, dataspace_end)

    ReadAttribute(name, dataspace, datatype_class, datatype_offset, position(io))
end


text:
"""
making it possible to invoke metrics in simple way
"""
module MainAbstractions
using CUDA
using Revise, Parameters, Logging, Test
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..BasicStructs
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates
using ..MainOverlap, ..RandIndex , ..ProbabilisticMetrics , ..VolumeMetric ,..InformationTheorhetic
using ..CUDAAtomicUtils, ..TpfpfnKernel, ..InterClassCorrKernel,..MeansMahalinobis
export prepareMetrics, calcMetricGlobal
"""
function created to prepare  kernel launch - required data is summarized in the ConfigurtationStruct
    confStruct- specified by the user - where is the summary of what metrics should be measured
    numberToLookFor - the number representing label which will be checked in the target arrays
"""
function prepareMetrics(confStruct::ConfigurtationStruct)

    dict = Dict()

    isAnyConfusionMatrixMetric = (confStruct.dice || confStruct.jaccard|| confStruct.gce|| confStruct.dice|| confStruct.vol|| confStruct.randInd|| confStruct.kc|| confStruct.mi || confStruct.vi)

    if(!confStruct.sliceWiseMetrics &&   isAnyConfusionMatrixMetric)
        dict[1]= prepareForconfusionTableMetricsNoSliceWise(confStruct)
    end
    if(confStruct.ic)
        dict[2]=InterClassCorrKernel.prepareInterClassCorrKernel()
    end
    if(confStruct.md)
        dict[3]= MeansMahalinobis.prepareMahalinobisKernel()
    end

    return dict
end

"""
executing metrics calculations
    preparedDict - dictionary created in the preparatory step
    confStruct- specified by the user - where is the summary of what metrics should be measured
    goldGPU,segmGPU - 3 dimensional arrays representing gold standard and output of our algorithm
    numberToLookFor - number representing the label of intrest in the arrays
"""
function calcMetricGlobal(preparedDict, confStruct,goldGPU,segmGPU, numberToLookFor)::ResultMetrics
    res = ResultMetrics()

    numberToLooFor=numberToLookFor
    isAnyConfusionMatrixMetric = (confStruct.dice || confStruct.jaccard|| confStruct.gce|| confStruct.dice|| confStruct.vol|| confStruct.randInd|| confStruct.kc|| confStruct.mi || confStruct.vi)

    if(isAnyConfusionMatrixMetric)
        argss,threads,blocks,metricsTuplGlobal=preparedDict[1]
        TpfpfnKernel.getTpfpfnData!(goldGPU,segmGPU   ,argss,threads,blocks,metricsTuplGlobal,numberToLooFor,confStruct)
        res.dice = metricsTuplGlobal[4][1]
        res.jaccard = metricsTuplGlobal[5][1]
        res.gce = metricsTuplGlobal[6][1]
        res.randInd = metricsTuplGlobal[7][1]
        res.kc = metricsTuplGlobal[8][1]
        res.vol = metricsTuplGlobal[9][1]
        res.mi = metricsTuplGlobal[10][1]
        res.vi = metricsTuplGlobal[11][1]
    end
    if(confStruct.ic)
        argsMain, threads,blocks, totalNumbOfVoxels=preparedDict[2]
        res.ic = InterClassCorrKernel.calculateInterclassCorr(goldGPU,segmGPU,threads,blocks,argsMain,numberToLooFor)
    end
    if(confStruct.md)
        argsMahal,threadsMahal ,blocksMahal= preparedDict[3]
        res.md = MeansMahalinobis.calculateMalahlinobisDistance(goldGPU,segmGPU,argsMahal,threadsMahal ,blocksMahal,numberToLooFor)
    end

    return res
    end#calcMetricGlobal

end


text:
module Foo

using Comonicon

"""
ArgParse example implemented in Comonicon.

# Arguments

- `x`: an argument, an argument

# Options

- `--opt1 <arg>`: an option
- `--opt2, -o <arg>`: another option

# Flags

- `--flag, -f`: a flag
"""
@main function main(x; opt1 = 1, opt2::Int = 2, flag = false)
    println("Parsed args:")
    println("flag=>", flag)
    println("arg=>", x)
    println("opt1=>", opt1)
    println("opt2=>", opt2)
end

end


text:
@testset "exgraph" begin
    let
        g = ExGraph()
        g = ExGraph(;ctx=[:foo => 42], x=1)
        @test g.ctx[:foo] == 42
        @test getvalue(g[1]) == 1
    end

    let
        g = ExGraph(:(z = (x + y)^2); x=1, y=1)
        @test length(g) == 5     # check that fuse_assigned() removed extra var
        @test varname(g[5]) == :z
        @test evaluate!(g) == 4

        g = ExGraph(:(z = (x .+ y).^2); x=ones(3), y=ones(3))
        @test evaluate!(g) == [4.0, 4.0, 4.0]
    end

    let
        ex = :(M = u'v)
        g = ExGraph(ex; u=rand(3), v=rand(3))
        @test getcategory(g[3]) == :call
        @test getexpr(g[3]) == :(transpose(u))
    end

    let
        g = ExGraph(:(x = u + v; z = 2x))
        nd = ExNode{:call}(:y, :(x - 1))
        insert!(g, 2, nd)
        @test varname(g[2]) == :y

        nds = ExGraph(:(t1 = u + x; t2 = t1 - x)).tape
        insert!(g, 3, nds)
        @test varname(g[3]) == :t1
        @test varname(g[4]) == :t2

        delete!(g, 3)
        @test length(g) == 5
        delete!(g, :t2)
        @test length(g) == 4

    end

    let
        # test graph with :opaque nodes
        g = ExGraph(:(y = 2x); x=rand(2))
        push!(g, ExNode{:opaque}(:z, :(x .* 3y)))
        rep_g = reparse(g)
        @test getvar(g[3]) == getvar(rep_g[3])
        @test evaluate!(g) == evaluate!(rep_g)

        rw_nd = rewrite(g[3], :(Z = X * Y), :(Z = X .* Y'); phs=[:X, :Y, :Z])
        @test getcategory(rw_nd) == :opaque
        @test getvalue(rw_nd) == nothing        # value should be cleared
    end
end

@testset "var renaming" begin

    ex = quote
        x = 0
        x = x + 1
        y = 4
        x = 2x + y
        y = 5
        z = x + y
    end

    g = ExGraph(ex)
    @test to_expr(g[end]) == :(z = x3 + y2)
end


text:
@doc """
SPLS Univariate Sparse Partial Least squares Algorithms

Algorithm first outlined in:
    Sparse and robust PLS for binary classification,
    I. Hoffmann, P. Filzmoser, S. Serneels, K. Varmuza,
    Journal of Chemometrics, 30 (2016), 153-162.

Parameters:
-----------
`eta`: float. Sparsity parameter in [0,1)

`n_components`: int, min 1. Note that if applied on data, n_components shall
    take a value <= minimum(size(x_data))

`verbose`: Bool: to print intermediate set of columns retained

`centre`: How to internally centre data. Accepts strings ("mean","median"),
    functions (e.g. Statistics.mean) or a location vector.

`scale`: How to internally scale data. Accepts strings ("std"),
        functions (e.g. Statistics.std) or a scale vector. Enter "none" both for
        centre and scale to work with the raw data.

`fit_algorithm`: Either "snipls", "sim_dense" or "sim_sparse".

`copy`: Bool, whether to copy data

Values:
-------
The mutable struct called by SNIPLS() will be populated with model results, such
    as coef_ (regression coefficients), x_scores_, etc., as well as estimated
    locations and scales for both X and y.

The module is consistent with the ScikitLearn API, e.g.

    import ScikitLearn.GridSearch:GridSearchCV
    gridsearch = GridSearchCV(spls.SPLS(), Dict(:eta => collect(0.9:-0.01:0.1),
                    :n_components => collect(1:4), :verbose => false))
    ScikitLearn.fit!(gridsearch,X,y)
    gridsearch.best_params_
    ScikitLearn.predict(gridsearch.best_estimator_,Xt)

Written by Sven Serneels.

""" ->
module spls

    using Statistics, DataFrames, Parameters, Random
    import ScikitLearnBase: BaseRegressor, predict, fit!, @declare_hyperparameters
    include("._sreg_utils.jl")
    export spls, SPLS, autoscale

    @with_kw mutable struct SPLS <: BaseRegressor
        eta::Number = .5
        n_components::Int = 1
        verbose::Bool = true
        centre = "mean"
        scale = "none"
        copy::Bool = true
        fit_algorithm = "snipls"
        X_Types_Accept = [DataFrame,Array{Number,2},Array{Float64,2},
                        Array{Float32,2},Array{Int32,2},Array{Int64,2},
                        Array{Union{Missing, Float64},2}]
        y_Types_Accept = [DataFrame,Array{Number,1},Array{Float64,1},
                        Array{Float32,1},Array{Int32,1},Array{Int64,1},
                        Array{Union{Missing, Float64},1}]
        X0 = nothing
        y0 = nothing
        x_weights_ = nothing
        x_loadings_ = nothing
        C_ = nothing
        x_scores_ = nothing
        coef_ = nothing
        coef_scaled_ = nothing
        intercept_ = nothing
        x_ev_ = nothing
        y_ev_ = nothing
        fitted_ = nothing
        residuals_ = nothing
        x_Rweights_ = nothing
        colret_ = nothing
        x_loc_ = nothing
        y_loc_ = nothing
        x_sca_ = nothing
        y_sca_ = nothing
        x_names = nothing
        y_name = nothing
    end

    @declare_hyperparameters(SPLS, [:n_components, :eta, :centre, :scale, :verbose, :fit_algorithm])

    @doc """

        Dummy function equivalent to directly creating a SNIPLS struct

        """ ->
    function call(;kwargs...)

        self = SSIMPLS()
        if length(kwargs) > 0
            allkeys = keys(kwargs)
            for j in allkeys
                setfield!(self,j,kwargs[j])
            end #for
        end #if
        return(self)
    end #snipls

    @doc """

        Fit SNIPLS model to data X and y.

        """ ->
    function fit!(self,X,y)

        @assert self.fit_algorithm in ["snipls"] "Algorithm has to be `snipls`; `sim_dense`, `sim_sparse` still under construction]"

        if typeof(self.centre)==String
            if self.centre in ["mean","median"]
                self.centre = getfield(Statistics,Symbol(self.centre))
            else
                @assert self.centre=="none" "Only supported strings for median:" * "\n" *
                     "'mean', 'median', 'none'" * "\n" *
                     "Alternatively pass a function"
                # other location estimators can be included
            end
        end

        if typeof(self.scale)==String
            if self.scale in ["std"]
                self.scale = getfield(Statistics,Symbol(self.scale))
            else
                @assert self.scale=="none" "Only supported strings for scale:" * "\n" *
                     "'std','none'" * "\n" *
                     "Alternatively pass a function"
            end
        end

        X_Type = typeof(X)
        y_Type = typeof(y)
        @assert X_Type in self.X_Types_Accept "Supply X data as DataFrame or Matrix"
        @assert y_Type in self.y_Types_Accept "Supply y data as DataFrame or Vector"

        if self.copy
            setfield!(self,:X0,deepcopy(X))
            setfield!(self,:y0,deepcopy(y))
        end

        X0 = X
        y0 = y

        if y_Type == DataFrame
            ynames = true
            y_name = names(y0)
            y0 = y[:,1]
        else
            ynames = false
        end

        if X_Type == DataFrame
            Xnames = true
            X_name = names(X0)
            X0 = Matrix(X0)
        else
            Xnames = false
            X_name = nothing
        end

        (n,p) = size(X0)
        ny = size(y0,1)
        @assert ny == n "Number of cases in X and y needs to agree"

        centring_X = autoscale(X0,self.centre,self.scale)
        Xs= centring_X.X_as_
        mX = centring_X.col_loc_
        sX = centring_X.col_sca_

        centring_y = autoscale(y0,self.centre,self.scale)
        ys= centring_y.X_as_
        my = centring_y.col_loc_
        sy = centring_y.col_sca_
        nys = sum(ys.^2)
        s = Xs'*ys

        if self.fit_algorithm == "sim_dense"
            # (W,P,R,C,T,B,Xev,yev,colret) = _fit_dense(self,X,y,n,p,Xs,mX,sX,ys,my,sy,nys,Xnames,X_name,s)
            # experimental
        elseif self.fit_algorithm == "sim_sparse"
            # (W,P,R,C,T,B,Xev,yev,colret) = _fit_sparse(self,X,y,n,p,Xs,mX,sX,ys,my,sy,nys,Xnames,X_name,s)
            # experimental
        elseif self.fit_algorithm == "snipls"
            (W,P,R,C,T,B,Xev,yev,colret) = _fit_snipls(self,X,y,n,p,Xs,mX,sX,ys,my,sy,nys,Xnames,X_name,s)
        end

        B_rescaled = (sy./sX)' .* B
        yp_rescaled = X0*B_rescaled

        if self.centre != "none"
            intercept = self.centre(y .- yp_rescaled)
        else
            intercept = mean(y .- yp_rescaled)
        end

        yfit = yp_rescaled .+ intercept
        r = y .- yfit
        setfield!(self,:x_weights_,W)
        setfield!(self,:x_Rweights_,R)
        setfield!(self,:x_loadings_,P)
        setfield!(self,:C_,C)
        setfield!(self,:x_scores_,T)
        setfield!(self,:coef_,B_rescaled)
        setfield!(self,:coef_scaled_,B)
        setfield!(self,:intercept_,intercept)
        setfield!(self,:x_ev_,Xev)
        setfield!(self,:y_ev_,yev)
        setfield!(self,:fitted_,yfit)
        setfield!(self,:residuals_,r)
        setfield!(self,:x_Rweights_,R)
        setfield!(self,:colret_,colret)
        setfield!(self,:x_loc_,mX)
        setfield!(self,:y_loc_,my)
        setfield!(self,:x_sca_,sX)
        setfield!(self,:y_sca_,sy)
        if Xnames
            setfield!(self,:x_names,X_name)
        end
        if ynames
            setfield!(self,:y_name,y_name)
        end

        return(self)

    end

    function _fit_snipls(self,X0,y0,n,p,Xs,mX,sX,ys,my,sy,nys,Xnames,X_name,wh=[])

        T = zeros(n,self.n_components)
        W = zeros(p,self.n_components)
        P = zeros(p,self.n_components)
        C = zeros(self.n_components,1)
        Xev = zeros(self.n_components,1)
        yev = zeros(self.n_components,1)
        B = zeros(p,1)
        oldgoodies = []
        Xi = Xs #deepcopy ?
        yi = ys

        if length(wh) == 0
            wh = Xs'*ys
        end

        for i in 1:self.n_components
            if i > 1
                wh =  Xi' * yi
            end
            (wh, goodies) =  _find_sparse(wh,self.eta)
            global goodies = union(oldgoodies,goodies)
            oldgoodies = goodies
            if length(goodies)==0
                print("No variables retained at" * string(i) * "latent variables" *
                      "and lambda = " * string(self.eta) * ", try lower lambda")
                break
            end
            elimvars = setdiff(1:p,goodies)
            wh[elimvars] .= 0
            th = Xi * wh
            nth2 = sum(th.^2)
            ch = (yi'*th) / nth2
            ph = (Xi' * Xi * wh) / nth2
            ph[elimvars] .= 0
            yi -= th*ch
            W[:,i] = wh
            P[:,i] = ph
            C[i] = ch
            T[:,i] = th
            Xi -= th * ph'
            Xev[i] = 100 - sum((Xs - T[:,1:i]*P[:,1:i]').^2) / sum(Xs.^2)*100
            yev[i] = 100 - sum((ys - T[:,1:i]*C[1:i]).^2) / nys *100
            if Xnames
                global colret = X_name[setdiff((1:p),elimvars)]
            else
                global colret = goodies
            end
            if self.verbose
                print("Variables retained for " * string(i) *
                        " latent variable(s):" *  "\n" * string(colret) * ".\n")
            end
        end

        if length(goodies) > 0
            R = W * inv(P'*W)
            all0 = findall(P.==0)
            R[all0] .= 0
            B = R * C
        else
            B = zeros(p,1)
            R = B
            T = zeros(n,self.n_components)
        end

        return((W,P,R,C,T,B,Xev,yev,colret))

    end

    @doc """

        Fit SNIPLS model to data X and y and only return the regression
        coefficients.

        """ ->
    function fit(self,X,y)

        if self.coef_ == nothing
            fit!(self,X,y)
        end

        return(self.coef_)

    end

    @doc """

        Predict responses for new predictor data.

        """ ->
    function predict(self,Xn)

        Xn_type = typeof(Xn)
        @assert Xn_type in self.X_Types_Accept "Supply new X data as DataFrame or Matrix"
        if Xn_type == DataFrame
            Xn = Matrix(Xn)
        end
        (n,p) = size(Xn)
        @assert p==length(self.coef_) "New data must have same number of columns as the ones the model has been trained with"
        return((Xn * self.coef_ .+ self.intercept_)[:,1])

    end

    @doc """

        Transform new predictor data to estimated scores.

        """ ->
    function transform(self,Xn)

        Xn_type = typeof(Xn)
        @assert Xn_type in self.X_Types_Accept "Supply new X data as DataFrame or Matrix"
        if Xn_type == DataFrame
            Xn = Matrix(Xn)
        end
        (n,p) = size(Xn)
        @assert p==length(self.coef_) "New data must have same number of columns as the ones the model has been trained with"
        Xnc = autoscale(Xn,self.x_loc_,self.x_sca_)
        return(Xnc*self.x_Rweights_)

    end

    @doc """

        Get all settings from an existing SNIPLS struct (also the ones not
            declared to ScikitLearn)

        """ ->
    function get_all_params(self::T, param_names=[]::Array{Any}) where{T}

        if length(param_names)==0

            params_dict = type2dict(self)

        else

            params_dict = Dict{Symbol, Any}()

            for name::Symbol in param_names
                params_dict[name] = getfield(self, name)
            end

        end

        params_dict

    end

    @doc """

        ScikitLearn similar function to set parameters in an existing SNIPLS
            struct, yet takes a Dict as an argument.

        Compare:
        ScikitLearn.set_params!(lnrj,eta=.3)
        snipls.set_params_dict!(lnrj,Dict(:eta => .3))

        """ ->
    function set_params_dict!(self::T, params; param_names=nothing) where {T}

        for (k, v) in params

            if param_names !== nothing && !(k in param_names)

                throw(ArgumentError("An estimator of type $T was passed the invalid hyper-parameter $k. Valid hyper-parameters: $param_names"))

            end

            setfield!(self, k, v)

        end

        print(self)

    end

    clone_param(v::Any) = v # fall-back

    clone_param(rng::Random.MersenneTwister) = deepcopy(rng) # issue #15698. Solved in 0.5

    function is_classifier(self) return(false) end

    @doc """

        ScikitLearn compatible function to clone an existing SNIPLS struct.

        """ ->
    function clone(self::T) where {T}

        kw_params = Dict{Symbol, Any}()

        # cloning the values is scikit-learn's default behaviour. It's ok?

        for (k, v) in get_params(self) kw_params[k] = clone_param(v) end

        return T(; kw_params...)

    end

end

max_stars_repo_name: SvenSerneels/spls

max_issues_repo_name: SvenSerneels/spls

max_forks_repo_name: SvenSerneels/spls