text:
#!/usr/bin/env julia

#===================================================
#
#  Section to transfer to new file rheodata.jl
#
===================================================#

#=
-------------------------
Log related functionality
-------------------------
=#
struct RheoLogItem
   action      # Nothing, or NamedTuple with fields:
               #     type::Symbol, funct::Symbol, params::NamedTuple, keywords::NamedTuple
   info        # NamedTuple for process output, comments, etc.
end

const RheoLog = Vector{RheoLogItem}

# simply pass text as comment
function RheoLogItem(s::String)
   return(RheoLogItem(Nothing,(comment=s,)))
end

# store info data in log from arbitrary keyword arguments
function RheoLogItem(;kwargs...)
   return(RheoLogItem(Nothing, kwargs.data))
end

function Base.show(io::IO, ::MIME"text/plain", rli::RheoLogItem)
    print("action = "); println(rli.action)
    print("info   = "); println(rli.info)
end

function Base.show(io::IO, ::MIME"text/plain", rl::RheoLog)
    println("RheoLog structure")
    i=0
    for rli in rl
       i+=1
       println("["*string(i)*"]")
       display(rli)
    end
end

# ----------
# Utilities
# ----------

"""
    rheoconvert(t)

Converts if required the scalar value or vector of real number `t` to the type `RheoFloat`.
"""
rheoconvert(t::Real) = RheoFloat(t)
rheoconvert(t::RheoFloat) = t
rheoconvert(t::Vector{T}) where T<:Real = convert(Vector{RheoFloat},t)
rheoconvert(t::Vector{RheoFloat}) = t

#=
-------------------------------
Time data related functionality
-------------------------------
=#
"""
    RheoTimeData(;σ::Vector{T1}, ϵ::Vector{T2}, t::Vector{T3}) where {T1<:Real, T2<:Real, T3<:Real}

`RheoTimeData` struct contains stress, strain and time data.

If preferred, an instance can be generated manually by just providing the three data
vectors in the right order, sampling type will be checked automatically.

# Fields

- `σ`: stress
- `ϵ`: strain
- `t`: time
- `log`: a log of struct's events, e.g. preprocessing
"""
struct RheoTimeData

    σ::Vector{RheoFloat}
    ϵ::Vector{RheoFloat}
    t::Vector{RheoFloat}

    log::Union{RheoLog,Nothing}

end

function RheoTimeData(;strain = RheoFloat[], ϵ::Vector{T1} = strain, stress = RheoFloat[], σ::Vector{T2} = stress, t::Vector{T3} = RheoFloat[], comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing)  where {T1<:Real, T2<:Real, T3<:Real}
    typecheck = check_time_data_consistency(t,ϵ,σ)
    RheoTimeData(convert(Vector{RheoFloat},σ), convert(Vector{RheoFloat},ϵ), convert(Vector{RheoFloat},t),
                    log === nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
end

function Base.show(io::IO, d::RheoTimeData)
    b = length(d.t) > 10
    n = b ? 10 : length(d.t)

    if d.t != RheoFloat[]
        print("t =\t")
        for i = 1:n
            print("$(d.t[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.ϵ != RheoFloat[]
        print("ϵ =\t")
        for i = 1:n
            print("$(d.ϵ[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.σ != RheoFloat[]
        print("σ =\t")
        for i = 1:n
            print("$(d.σ[i])\t")
        end
        println( b ? "..." : "")
    end
end

@enum TimeDataType invalid_time_data=-1 time_only=0 strain_only=1 stress_only=2 strain_and_stress=3

function check_time_data_consistency(t,e,s)
    @assert (length(t)>0)  "Time data empty"

    sdef=(s != RheoFloat[])
    edef=(e != RheoFloat[])
    if (sdef && edef)
        @assert (length(s)==length(t)) && (length(e)==length(t)) "Time data length inconsistent"
        return strain_and_stress
    end

    if (sdef && (!edef))
        @assert (length(s)==length(t)) "Time data length inconsistent"
        return stress_only
    end

    if (edef && (!sdef))
        @assert (length(e)==length(t)) "Time data length inconsistent"
        return strain_only
    end

    if ((!edef) && (!sdef))
        return time_only
    end

    return invalid_time_data

end

function rheotimedatatype(d::RheoTimeData)
    return check_time_data_consistency(d.t,d.ϵ,d.σ)
end

"""
    hastime(d::RheoTimeData)

returns 'true' if d contains a time array
"""
function hastime(d::RheoTimeData)
    return (d.t != RheoFloat[])
end

"""
    hasstress(d::RheoTimeData)

returns 'true' if d contains stress data
"""
function hasstress(d::RheoTimeData)
    return (d.σ != RheoFloat[])
end

"""
    hasstrain(d::RheoTimeData)

returns 'true' if d contains strain data
"""
function hasstrain(d::RheoTimeData)
    return (d.ϵ != RheoFloat[])
end

"""
    gettime(d::RheoTimeData)

returns the time vector.
"""
function gettime(d::RheoTimeData)
    return (d.t)
end

"""
    getstrain(d::RheoTimeData)

returns the strain vector if d contains strain data.
"""
function getstrain(d::RheoTimeData)
    return (d.ϵ)
end

"""
    getstress(d::RheoTimeData)

returns the stress vector if d contains stress data.
"""
function getstress(d::RheoTimeData)
    return (d.σ)
end

# Utilities for in-place functions
#
# Internal - Not be exposed to the user.

function _settime!(data::RheoTimeData, t::Vector{T}) where {T<:Real}
    if hastime(data)
      empty!(data.t)
    end
    append!(data.t, rheoconvert(t))
end

function _setstrain!(data::RheoTimeData, ϵ::Vector{T}) where {T<:Real}
    if hasstrain(data)
      empty!(data.ϵ)
    end
    append!(data.ϵ, rheoconvert(ϵ))
end

function _setstress!(data::RheoTimeData, σ::Vector{T}) where {T<:Real}
    if hasstress(data)
      empty!(data.σ)
    end
    append!(data.σ, rheoconvert(σ))
end

# Constants that are useful for certain processing function to determine which module to use.
@enum LoadingType strain_imposed=1 stress_imposed=2

function +(d1::RheoTimeData, d2::RheoTimeData)

    type1 = rheotimedatatype(d1)
    type2 = rheotimedatatype(d2)
    @assert (type1!=invalid_time_data) "Addition error: first parameter invalid"
    @assert (type2!=invalid_time_data) "Addition error: second parameter invalid"
    @assert (type1==type2) "Addition error: parameters inconsistent"
    @assert (type1!=time_only) "Addition error: time only data cannot be added"
    @assert (d1.t == d2.t) "Addition error: timelines inconsistent"

    log =   if (d1.log === nothing) || (d2.log === nothing)
                nothing
            else
                [ RheoLogItem( (type = :source, funct = :+, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ()  ) ]
            end

    if (type1==strain_only) && (type2==strain_only)
        return RheoTimeData(RheoFloat[], d1.ϵ+d2.ϵ, d1.t, log)
    end

    if (type1==stress_only) && (type2==stress_only)
        return RheoTimeData(d1.σ+d2.σ, RheoFloat[], d1.t, log)
    end

    if (type1==strain_and_stress) && (type2==strain_and_stress)
        return RheoTimeData(d1.σ+d2.σ, d1.ϵ+d2.ϵ, d1.t, log)
    end

end

function -(d1::RheoTimeData, d2::RheoTimeData)

    type1 = rheotimedatatype(d1)
    type2 = rheotimedatatype(d2)
    @assert (type1!=invalid_time_data) "Subtraction error: first parameter invalid"
    @assert (type2!=invalid_time_data) "Subtraction error: second parameter invalid"
    @assert (type1==type2) "Subtraction error: parameters inconsistent"
    @assert (type1!=time_only) "Subtraction error: time only data cannot be added"
    @assert (d1.t == d2.t) "Subtraction error: timelines inconsistent"

    log =   if (d1.log === nothing) || (d2.log === nothing)
                nothing
            else
                [ RheoLogItem( (type = :source, funct = :-, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ()  ) ]
            end

    if (type1==strain_only) && (type2==strain_only)
        return RheoTimeData(RheoFloat[], d1.ϵ-d2.ϵ, d1.t, log)
    end

    if (type1==stress_only) && (type2==stress_only)
        return RheoTimeData(d1.σ-d2.σ, RheoFloat[], d1.t, log)
    end

    if (type1==strain_and_stress) && (type2==strain_and_stress)
        return RheoTimeData(d1.σ-d2.σ, d1.ϵ-d2.ϵ, d1.t, log)
    end

end

function -(d::RheoTimeData)

    type = rheotimedatatype(d)
    @assert (type!=invalid_time_data) "unary - error: parameter invalid"
    @assert (type!=time_only) "unary - error: time only data cannot be manipulated this way"

    log = d.log === nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:-, params=(), keywords=() ), () ) ]

    return RheoTimeData(-d.σ, -d.ϵ, d.t, log)

end

# Required to execute +/- operators for logs

function +(rl1::RheoLog, rl2::RheoLog)
    return(rheologrun(rl1) + rheologrun(rl2))
end

function -(rl1::RheoLog, rl2::RheoLog)
    return(rheologrun(rl1) - rheologrun(rl2))
end

function *(operand::Real, d::RheoTimeData)

    type = rheotimedatatype(d)
    @assert (type!=invalid_time_data) "* error: parameter invalid"
    @assert (type!=time_only) "* error: time only data cannot be manipulated this way"

    log = d.log === nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:*, params=(operand = operand), keywords=() ), () ) ]

    return( RheoTimeData(operand .* d.σ, operand .* d.ϵ, d.t, log) )
end

function *(d::RheoTimeData, operand::Real)
    return(operand * d)
end

function |(d1::RheoTimeData, d2::RheoTimeData)
    # get types and type-check
    type1 = rheotimedatatype(d1)
    type2 = rheotimedatatype(d2)

    @assert ((type1==strain_only && type2==stress_only) || (type1==stress_only && type2==strain_only)) "One RheoTimeData must be stress_only and the other must be strain_only"

    # call recursively if not expected way round
    (type1==strain_only && type2==stress_only) && return d2|d1

    # expected case, type1==stress_only, type2==strain_only
    @assert d1.t==d2.t "Data to be combined must have same time-series"

    log = ((d1.log === nothing) || (d2.log === nothing)) ? nothing : [RheoLogItem((type = :source, funct = :|, params = (rl1 = d1.log, rl2 = d2.log), keywords = ()), ())]

    return RheoTimeData(d1.σ, d2.ϵ, d1.t, log)
end

#=
------------------------------------
Frequency data related functionality
------------------------------------
=#
"""
    RheoFreqData(Gp::Vector{T1}, Gpp::Vector{T2}, ω::Vector{T3}, log::OrderedDict{Any,Any}) where {T1<:Real, T2<:Real, T3<:Real}

`RheoFreqData` contains storage modulus, loss modulus and frequency data.

If preferred, an instance can be generated manually by just providing the three data
vectors in the right order.

# Fields

- `Gp`: storage modulus
- `Gpp`: loss modulus
- `ω`: frequency
- `log`: a log of struct's events, e.g. preprocessing
"""
struct RheoFreqData
    Gp::Vector{RheoFloat}
    Gpp::Vector{RheoFloat}
    ω::Vector{RheoFloat}

    log::Union{RheoLog,Nothing}
end

function RheoFreqData(;Gp::Vector{T1} = RheoFloat[], Gpp::Vector{T2} = RheoFloat[], omega = RheoFloat[], ω::Vector{T3} = omega, comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing)  where {T1<:Real, T2<:Real, T3<:Real}

    typecheck = check_freq_data_consistency(ω,Gp,Gpp)
    RheoFreqData(convert(Vector{RheoFloat},Gp), convert(Vector{RheoFloat},Gpp), convert(Vector{RheoFloat},ω),
    log === nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck,)))]     )
end

function Base.show(io::IO, d::RheoFreqData)
    b = length(d.ω) > 10
    n = b ? 10 : length(d.ω)

    if d.ω != RheoFloat[]
        print("ω =\t")
        for i = 1:n
            print("$(d.ω[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.Gp != RheoFloat[]
        print("Gp =\t")
        for i = 1:n
            print("$(d.Gp[i])\t")
        end
        println( b ? "..." : "")
    end
    if d.Gpp != RheoFloat[]
        print("Gpp =\t")
        for i = 1:n
            print("$(d.Gpp[i])\t")
        end
        println( b ? "..." : "")
    end
end

@enum FreqDataType invalid_freq_data=-1 freq_only=0 with_modulus=1

function check_freq_data_consistency(o,gp,gpp)
    @assert (length(o)>0)  "Freq data empty"

    gpdef=(gp != RheoFloat[])
    gppdef=(gpp != RheoFloat[])
    if (gpdef && gppdef)
        @assert (length(gp)==length(o)) && (length(gpp)==length(o)) "Data length inconsistent"
        return with_modulus
    end

    if ((!gpdef) && (!gppdef))
        return freq_only
    end

    return invalid_freq_data

end

function rheofreqdatatype(d::RheoFreqData)
    return check_freq_data_consistency(d.ω,d.Gp,d.Gpp)
end

"""
    hasfreq(d::RheoFreqData)

returns 'true' if d contains a frequency array
"""
function hasfreq(d::RheoFreqData)
    return (d.ω != RheoFloat[])
end

"""
    hasmodulus(d::RheoFreqData)

returns 'true' if d contains a frequency array
"""
function hasmodulus(d::RheoFreqData)
    (d.Gp != RheoFloat[]) && (d.Gpp != RheoFloat[])
end

"""
    getomega(d::RheoFreqData)

returns the frequency ω vector.
"""
function getfreq(d::RheoFreqData)
    return (d.ω)
end

"""
    getstorage(d::RheoFreqData)

returns the storage modulus vector.
"""
function getstorage(d::RheoFreqData)
    return (d.Gp)
end

"""
    getloss(d::RheoFreqData)

returns the loss modulus vector.
"""
function getloss(d::RheoFreqData)
    return (d.Gpp)
end

# Utilities for in-place functions
#
# Internal - Not be exposed to the user.

function _setfreq!(data::RheoTimeData, ω::Vector{T}) where {T<:Real}
    if hasfreq(data)
      empty!(data.ω)
    end
    append!(data.ω, rheoconvert(ω))
end

function _setmodulus!(data::RheoTimeData, Gp::Vector{T1}, Gpp::Vector{T2}) where {T1<:Real,T2<:Real}
    if hasmodulus(data)
      empty!(data.Gp)
      empty!(data.Gpp)
    end
    append!(data.Gp, rheoconvert(Gp))
    append!(data.Gpp, rheoconvert(Gpp))
end

#=
---------------------------------------
Log-based data generation functionality
---------------------------------------
=#
function rheologrun(rli::RheoLogItem, d=nothing)

      if typeof(rli.action) <: NamedTuple && :type in keys(rli.action)
         type=rli.action.type
         if type==:source && d===nothing
            return(eval(rli.action.funct)(rli.action.params...;rli.action.keywords...))
         elseif type==:process && d!==nothing
            return(eval(rli.action.funct)(d,rli.action.params...;rli.action.keywords...))
         elseif type==:analysis && d!==nothing
            return(eval(rli.action.funct)(d,rli.action.params...;rli.action.keywords...))
         end
      end
   println(rli.info)
end

"""
    rheologrun(log::RheoLog, d::Union{RheoTimeData,RheoFreqData,Nothing} = nothing)

execute all actions from the log. It applies them to the data `d` provided, or use the log's first action to reload/recreate it otherwise.
"""
function rheologrun(arli::RheoLog, d::Union{RheoTimeData,RheoFreqData,Nothing} = nothing)

  # check first item is a source item
  if d === nothing
        @assert typeof(arli[1].action) <: NamedTuple && :type in keys(arli[1].action) && arli[1].action.type == :source
                "Source missing in RheoLog"
        end

  for rli in arli
      if typeof(rli.action) <: NamedTuple && :type in keys(rli.action)
          type=rli.action.type
          if type==:analysis && d!==nothing
              rheologrun(rli, d)
          else  d=rheologrun(rli, d)
          end
      else
          println(rli.info)
      end
  end
  return(d)
end

"""
    showlog(d::Union{RheoTimeData,RheoFreqData})

shows the record of operations (`RheoLog`) on rheological data.
"""
function showlog(d::Union{RheoTimeData,RheoFreqData})
    if d.log !== nothing
        display(d.log)
    else
        println("No log data available")
    end
end

#===================================================
#
#  Section to transfer to new file rheomodel.jl
#
===================================================#

#=
---------------------------
Model related functionality
---------------------------
=#

# Function types for the moduli with free params
const FWScaFree = FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}
const FWVecFree = FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}
# Function types for the moduli with fixed params
const FWScaFixed = FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
const FWVecFixed = FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
# Function types for the constraint function
const FWConstraint = FunctionWrapper{Bool,Tuple{Vector{RheoFloat}}}

# constant NaN functions for default parameters and testing existence of moduli functions
const fwscafreeNaN = ((t,p)->NaN) |> FWScaFree
const fwvecfreeNaN = ((t,p)->fill(RheoFloat(NaN), length(t))) |> FWVecFree
const fwscafixedNaN = (t->NaN) |> FWScaFixed
const fwvecfixedNaN = (t->fill(RheoFloat(NaN), length(t))) |> FWVecFixed

struct _RheoModel{TSca,TVec}

	name::String
	freeparams::Tuple{Vararg{Symbol, N}} where N
	fixedparams::NamedTuple

	# Moduli functions specialised for scalar and array values
	_G::TSca
	_Ga::TVec
	_J::TSca
	_Ja::TVec
	_Gp::TSca
	_Gpa::TVec
	_Gpp::TSca
	_Gppa::TVec

	_constraint::FWConstraint
	_Gramp::Bool

	info::String
	expressions::NamedTuple

end

#
#  Specialised forms of _RheoModel
#  Dediated constructors are defined below
#

const RheoModelClass = _RheoModel{FWScaFree,FWVecFree}

const RheoModel = _RheoModel{FWScaFixed,FWVecFixed}

function Base.show(io::IO, m::_RheoModel)
    print(io, "\nModel name: $(m.name)")
    ps=join([string(s) for s in m.freeparams], ", ", " and ")
    print(io, "\n\nFree parameters: $ps\n")
    if length(m.fixedparams)>0
        print(io, "\n\nFixed parameters: $(m.fixedparams)\n")
    end
    print(io, m.info)
    return
end

# Cool replacement function inspired from
# https://stackoverflow.com/questions/29778698/julia-lang-metaprogramming-turn-expression-into-function-with-expression-depend
# This probably could go in its own module as this is very useful.
expr_replace!(ex, s, v) = ex
expr_replace!(ex::Symbol, s, v) = s == ex ? v : ex

function expr_replace!(ex::Expr, s, v)
    for i=1:length(ex.args)
        ex.args[i] = expr_replace!(ex.args[i], s, v)
    end
    return ex
end

function expr_replace(ex, s, v)
    expr_replace!(copy(ex), s, v)
end

function expr_replace(ex::Expr, nt::NamedTuple)
    e = copy(ex)
    k=keys(nt)
    v=values(nt)
    for i in 1:length(nt)
        expr_replace!(e, k[i], v[i])
    end
    return e
end

function expr_replace(ex::Nothing, nt::NamedTuple)
	return nothing
end

#
#  Look for a specific symbol in an expression
#
expr_search(ex, s) = false
expr_search(ex::Symbol, s) = s == ex ? true : false

@inline function expr_search(ex::Expr, s)
        any( expr_search(e, s) for e in ex.args )
end

@inline function _replace_symbols_with_array(e::Expr, psymbs::Tuple)
    return( expr_replace(e,NamedTuple{psymbs}([Meta.parse(":(p_arr[$i])").args[1] for i in 1:length(psymbs)]) ) )
end

@inline function _expr_cleanup!(ex::Expr)
    filter!(e -> typeof(e)!=LineNumberNode,  ex.args)
    return()
end

@inline function _expr_cleanup(ex::Expr)
    e = copy(ex)
    _expr_cleanup!(e)
    return(e)
end

modulusexists(modulus::FWScaFree) =  (modulus !== fwscafreeNaN)
modulusexists(modulus::FWVecFree) =  (modulus !== fwvecfreeNaN)
modulusexists(modulus::FWScaFixed) =  (modulus !== fwscafixedNaN)
modulusexists(modulus::FWVecFixed) =  (modulus !== fwvecfixedNaN)

function modulusexists(modulus_single_input)
   # 5.0 is not meaningful, just an input to check that the function does not return NaN
    # at a non-zero value, which is the default function output for moduli that have
    # not been defined.
    return !isnan(modulus_single_input(5.0))
end

#
#  Create the scalar and vector function wrappers for the moduli.
#

function _buildmoduli_t(G::Expr, psymbs::Tuple)
    # Replace symbols by elements in array
	Ge = _replace_symbols_with_array(G, psymbs)
    _expr_cleanup!(Ge)
    isexprmultiline = (length(Ge.args) > 1)

    # Detect if expression depends on variable or if it is a constant
    # This is because the @. macro fails to vectorize if there is no vector to iterate on
    # It also fails on multiline expressions
    islaplace, isconstant = expr_search(Ge, :t) ? (false, false) : expr_search(Ge, :s) ? (true, false) : (false, true)

    if isconstant  # constant value returned
        @eval return( (     ((t,p_arr) -> $Ge)      |> FWScaFree,
                            ((t,p_arr) -> fill(rheoconvert($Ge), length(t)))  |> FWVecFree ) )
#    elseif islaplace
#        return()
    elseif isexprmultiline
	    @eval return( (     ((t,p_arr) -> $Ge)      |> FWScaFree,
		    			    ((ta,p_arr) -> [$Ge for t in ta] )  |> FWVecFree ) )
    else
	    @eval return( (     ((t,p_arr) -> $Ge)      |> FWScaFree,
		    			    ((t,p_arr) -> @. $Ge )  |> FWVecFree ) )
    end
end

function _buildmoduli_t(G::Expr)
    Ge = _expr_cleanup(G)
    isexprmultiline = (length(Ge.args) > 1)

    islaplace, isconstant = expr_search(Ge, :t) ? (false, false) : expr_search(Ge, :s) ? (true, false) : (false, true)

    if isconstant     # constant value returned
        @eval return( ( (t -> $Ge)                             |> FWScaFixed,
                        (t -> fill(rheoconvert($Ge), length(t)))  |> FWVecFixed ) )
#    elseif islaplace
#        return()
    elseif isexprmultiline
        @eval return( ( (t -> $Ge)          |> FWScaFixed,
                (ta -> [$Ge for t in ta] )  	|> FWVecFixed ) )
    else
        @eval return( ( (t -> $Ge)          |> FWScaFixed,
                        (t -> @. $Ge )  	|> FWVecFixed ) )
    end
end

function _buildmoduli_ω(G::Expr, psymbs::Tuple)
    # Replace symbols by elements in array
	Ge = _replace_symbols_with_array(G, psymbs)
    _expr_cleanup!(Ge)
    isexprmultiline = (length(Ge.args) > 1)

    # Detect if expression depends on variable or if it is a constant
    # This is because the @. macro fails to vectorize if there is no vector to iterate on
    # It also fails on multiline expressions
    isconstant = !expr_search(Ge, :ω)

    if isconstant  # constant value returned
        @eval return( (     ((ω,p_arr) -> $Ge)      |> FWScaFree,
                            ((ω,p_arr) -> fill(rheoconvert($Ge), length(ω)))  |> FWVecFree ) )
    elseif isexprmultiline
	    @eval return( (     ((ω,p_arr) -> $Ge)      |> FWScaFree,
		    			    ((ωa,p_arr) -> [$Ge for ω in ωa] )  |> FWVecFree ) )
    else
	    @eval return( (     ((ω,p_arr) -> $Ge)      |> FWScaFree,
		    			    ((ω,p_arr) -> @. $Ge )  |> FWVecFree ) )
    end
end

function _buildmoduli_ω(G::Expr)
    Ge = _expr_cleanup(G)
    isexprmultiline = (length(Ge.args) > 1)

    isconstant = !expr_search(Ge, :ω)

    if isconstant     # constant value returned
        @eval return( ( (ω -> $Ge)                             |> FWScaFixed,
                        (ω -> fill(rheoconvert($Ge), length(ω)))  |> FWVecFixed ) )
#    elseif islaplace
#        return()
    elseif true # isexprmultiline
        @eval return( ( (ω -> $Ge)          |> FWScaFixed,
                (ωa -> [$Ge for ω in ωa] )  	|> FWVecFixed ) )
    else
        @eval return( ( (ω -> $Ge)          |> FWScaFixed,
                        (ω -> @. $Ge )  	|> FWVecFixed ) )
    end
end

function _buildmoduli_t(G::Nothing, psymbs::Tuple)
	return( (       fwscafreeNaN, fwvecfreeNaN ) )
end

function _buildmoduli_t(G::Nothing)
	return( (       fwscafixedNaN, fwvecfixedNaN ) )
end

function _buildmoduli_ω(G::Nothing, psymbs::Tuple)
	return( (       fwscafreeNaN, fwvecfreeNaN ) )
end

function _buildmoduli_ω(G::Nothing)
	return( (       fwscafixedNaN, fwvecfixedNaN ) )
end

function _buildconstraint(constraint::Expr, psymbs::Tuple)
    # Replace symbols by elements in array
	constraint_e = expr_replace(constraint,NamedTuple{psymbs}([Meta.parse(":(p_arr[$i])").args[1] for i in 1:length(psymbs)]) )
	@eval return( (     (p_arr -> $constraint_e)      |> FWConstraint ) )
end

function _buildconstraint(constraint::Expr)
	@eval return( (     (p_arr -> $constraint)      |> FWConstraint ) )
end

"""
    RheoModelClass(;name::String, p::Tuple, G,J,Gp,Gpp, constraint, info, G_ramp)

`RheoModelClass` is a complex data structure that contains all the relevant data to
allow RHEOS to fit models and make predictions once parameters are set.
RHEOS possesses already a large number of such data structures to represent common rheological models,
such as Maxwell or KelvinVoigt.

Users are not expected to directly manipulate the content of a RheoModelClass object, but will
pass it to relevant functions for fitting and numerically evaluating visco-elastic moduli.

The model name and its parameters can however be printed and display on the REPL.

# Example
```@example
julia> Maxwell

Model name: maxwell

Free parameters: η and k

                ___
            _____| |________╱╲  ╱╲  ╱╲  ___
                _|_|          ╲╱  ╲╱  ╲╱
                  η                  k

```
"""
function RheoModelClass(;name::String,
        p::Tuple,
        G = nothing,
        J = nothing,
        Gp = nothing,
        Gpp = nothing,
        constraint::Expr = quote true end,
        info="",
        # flag to indicate use of integral forms of the relaxation modulus.
               # using the ramp response instead of the step response helps avoid singularities
        G_ramp::Bool = false,
        )

    # Building expressions tuple to store data provided to constructor
    expressions = (G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint)

    return(RheoModelClass(name, p, NamedTuple{}(),
        _buildmoduli_t(G,p)..., _buildmoduli_t(J,p)...,
        _buildmoduli_ω(Gp,p)..., _buildmoduli_ω(Gpp,p)...,
        _buildconstraint(constraint,p),
        G_ramp, info, expressions) )
end

#
#  Function used to check and convert any named tupple of parameters into an array ordered as expected by the moduli functions of the RheoModelClass
#

function check_and_reorder_parameters(nt::NamedTuple, params_ordered_list::Tuple; err_string::String="calling function")
	# check that every parameter in m exists in the named tuple nt
	@assert all(i->i in keys(nt),params_ordered_list) "Missing parameter(s) in " * err_string
	# check that no extra parameters have been provided
	@assert length(params_ordered_list) == length(nt) "Mismatch number of model parameters and parameters provided in " * err_string

	p = map(i->RheoFloat(nt[i]), params_ordered_list)
    rheoconvert([p...])
end

function _freeze_params(m::RheoModelClass, nt0::NamedTuple)

	# check that every parameter to fix actually exists in the model
	@assert all( i-> i in m.freeparams,keys(nt0)) "A parameter to freeze is not involved in the model"
	# convert values format for consistency
	nt=NamedTuple{keys(nt0)}([RheoFloat(i) for i in nt0])

	# update list of fixed parameters
	fixedparams = merge(m.fixedparams, nt)

	# create array of remaining variables
	freeparams = Tuple( filter(s -> !(s in keys(nt)),m.freeparams) )

	# This section creates moduli expressions with material parameters
	# replaced by specific values.

	G = expr_replace(m.expressions.G, nt)
	J = expr_replace(m.expressions.J, nt)
	Gp = expr_replace(m.expressions.Gp, nt)
	Gpp = expr_replace(m.expressions.Gpp, nt)
	constraint = expr_replace(m.expressions.constraint, nt)

	return( (freeparams=freeparams, fixedparams=fixedparams, G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint) )

end

"""
    freezeparams(m::RheoModelClass, nt0::NamedTuple)

Return a new `RheoModelClass` with some of the parameters frozen to specific values

# Arguments

- `m`: original `RheoModelClass`
- `nt0`: named tuple with values for each parameter to freeze

# Example
```@example
julia> SLS2_mod = freezeparams( SLS2, (G₀=2,η₂=3.5))
[...]

julia> SLS2.G(1,[2,1,2,3,3.5])
3.8796492

julia> SLS2_mod.G(1,[1,2,3])
3.8796492
```

"""
function freezeparams(m::RheoModelClass, nt0::NamedTuple)

	freeparams, fixedparams, G,J,Gp,Gpp,constraint = _freeze_params(m, nt0)

	# Check that some free params remains
	@assert length(freeparams) > 0  "All parameters are set. Build a RheoModel instead."

	# Building expressions tuple to store data provided to constructor
	expressions = (G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint)

	return(RheoModelClass(m.name, freeparams, fixedparams,
		_buildmoduli_t(G,freeparams)..., _buildmoduli_t(J,freeparams)...,
		_buildmoduli_ω(Gp,freeparams)..., _buildmoduli_ω(Gpp,freeparams)...,
		_buildconstraint(constraint,freeparams),
        m._Gramp, m.info, expressions) )
end

function freezeparams(m::RheoModelClass; kwargs...)
    return(freezeparams(m,symbol_to_unicode(kwargs.data)))
end

function freeze_params(m::RheoModelClass, nt0::NamedTuple)
    println("freeze_params is deprecated - use freezeparams instead.")
    return(freezeparams(m,nt0))
end

function freeze_params(m::RheoModelClass; kwargs...)
    println("freeze_params is deprecated - use freezeparams instead.")
    return(freezeparams(m,symbol_to_unicode(kwargs.data)))
end

"""
    RheoModel(m::RheoModelClass, nt0::NamedTuple)

`RheoModel` represents a rheological model with set parameters. They are obtained by fitting a model to data using `modelfit`,
or by specialising the relevant RheoModelClass by prescribing its parameters. Parameters can be provided as a named tuple or keyword arguments.

`RheoModel` objects can then be used to simulate the response to an arbitrary input using `modelpredict`, and access the values of the moduli functions.

# Example
```@example
julia> model = RheoModel(Maxwell, (k=1, η=2.))
[...]

julia> model = RheoModel(Maxwell, k=1, η=2.)
[...]
```
"""
function RheoModel(m::RheoModelClass, nt0::NamedTuple)

	freeparams, fixedparams, G,J,Gp,Gpp,constraint = _freeze_params(m, nt0)

	# Check all free params are set
	@assert length(freeparams) == 0  "Some parameters need to be set: $freeparams"

	# Building expressions tuple to store data provided to constructor
	expressions = (G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint)

	return(RheoModel(m.name, freeparams, fixedparams,
		_buildmoduli_t(G)..., _buildmoduli_t(J)...,
        _buildmoduli_ω(Gp)..., _buildmoduli_ω(Gpp)...,
		_buildconstraint(constraint),
		m._Gramp, m.info, expressions) )
end

function RheoModel(m::RheoModelClass; kwargs...)
    return(RheoModel(m,symbol_to_unicode(kwargs.data)))
end

"""
    getparams(m::RheoModel; unicode=true)

`getparams` return the list of model parameters with their values as a NamedTuple.
If `unicode` is set to `false`, the unicode symbols are converted their text equivalent.

# Example
```@example
julia> model = RheoModel(Maxwell, k=1, η=2.)
[...]

julia> getparams(m)
(k = 1.0, η = 2.0)

julia> getparams(m,unicode=false)
(k = 1.0, eta = 2.0)
```
"""
function getparams(m::RheoModel; unicode=true)
    if unicode
        return(m.fixedparams)
    else
        return(unicode_to_text(m.fixedparams))
    end
end

#=
function freeze_params(m::RheoModelClass, nt0::NamedTuple)

    # check that every parameter in m exists in the named tuple nt
    @assert all( i-> i in m.params,keys(nt0)) "A parameter to freeze is not involved in the model"
    # convert values format for consistency
    nt=NamedTuple{keys(nt0)}([RheoFloat(i) for i in nt0])

    # create array of remaining variables
    p = filter(s -> !(s in keys(nt)),m.params)

    name="$(m.name) with set parameters: $nt"
    info = m.info

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))

    # This section creates moduli expressions with material parameters
    # replaced by specific values.

    G = expr_replace(m.expressions.G, nt)
    J = expr_replace(m.expressions.J, nt)
    Gp = expr_replace(m.expressions.Gp, nt)
    Gpp = expr_replace(m.expressions.Gpp, nt)
    constraint = expr_replace(m.expressions.constraint, nt)

    expressions=NamedTuple{(:G, :J, :Gp, :Gpp, :constraint, :Ga_safe, :Ja_safe)}( ( G, J, Gp, Gpp, constraint, m.expressions.Ga_safe, m.expressions.Ja_safe ) )

    @eval return( RheoModelClass($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((t,params) -> begin $unpack_expr; $J; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((ω,params) -> begin $unpack_expr; $Gp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        ((ω,params) -> begin $unpack_expr; $Gpp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}},
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}},
        (params -> begin $unpack_expr; $constraint; end) |> FunctionWrapper{Bool,Tuple{Vector{RheoFloat}}},
        $info, $expressions)   )
end

struct RheoModel

    _G::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Ga::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
    _J::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Ja::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
    _Gp::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Gpa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
    _Gpp::FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
    _Gppa::FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}

    expressions::NamedTuple
    params::NamedTuple
    info::String    # do we need this?
    # log::OrderedDict{Any,Any}

end

function RheoModel(m::RheoModelClass; kwargs...)
    return(RheoModel(m,kwargs.data))
end

function RheoModel(m::RheoModelClass, nt0::NamedTuple)

    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(nt0, m.params, err_string = "model definition")
    nt=NamedTuple{Tuple(m.params)}(p)

    # string printed
    info = string("\nModel: $(m.name)\n\nParameter values: $nt \n",m.info)

    # This attaches expressions to the models with parameters substituted by values.
    # Future development: this could be disabled with a key word if processing time is an issue
    G = expr_replace(m.expressions.G, nt)
    J = expr_replace(m.expressions.J, nt)
    Gp = expr_replace(m.expressions.Gp, nt)
    Gpp = expr_replace(m.expressions.Gpp, nt)

    expressions=NamedTuple{(:G, :J, :Gp, :Gpp, :Ga_safe, :Ja_safe)}( ( G, J, Gp, Gpp, m.expressions.Ga_safe, m.expressions.Ja_safe ) )

    @eval return( RheoModel(
    (t -> begin $G; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ta -> begin [$G for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    (t -> begin $J; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ta -> begin [$J for t in ta]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    (ω -> begin $Gp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ωa -> begin [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    (ω -> begin $Gpp; end) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}},
    (ωa -> begin [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}},
    $expressions, $nt, $info) )
end

function Base.show(io::IO, m::RheoModel)
    print(io,m.info)
end

=#

#
# These are exported functions enabling users to access the different moduli functions for each model
#

"""
    relaxmod(m[, t, params])

provide access to the relaxation modulus (G) of a given model m. It can be used with a broad range of inputs.

* When a time value is provided (or an array of time values), the function returns the corresponding value(s) of the relaxation modulus.

    relaxmod(m::RheoModel, time (single value or array))

    relaxmod(m::RheoModelClass, time (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    relaxmod(Maxwell, 1, k=1., η=1)

    relaxmod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, relaxmod returns the relaxation function itself.

    relaxmod(m::RheoModel)

    relaxmod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    G = relaxmod(m)
    G([0,1,2])
"""
function relaxmod(m::RheoModel, t::Number)
    m._G(t)
end

function relaxmod(m::RheoModel, ta::Vector{T}) where T <: Number
    m._Ga(ta)
end

relaxmod(m::RheoModel) = x -> relaxmod(m,x)

function relaxmod(m::RheoModelClass, t::Number, params::Vector{T}) where T <: Number
    m._G(t, params)
end

function relaxmod(m::RheoModelClass, ta::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    m._Ga(ta, params)
end

function relaxmod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.freeparams, err_string = "relaxmod")
    relaxmod(m, x, p)
end

function relaxmod(m::RheoModelClass, x; kwargs...)
    relaxmod(m, x, symbol_to_unicode(kwargs.data))
end

relaxmod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> relaxmod(m, x, params)
relaxmod(m::RheoModelClass, params::NamedTuple) =  x -> relaxmod(m, x, params)
relaxmod(m::RheoModelClass; kwargs...) =  x -> relaxmod(m, x, symbol_to_unicode(kwargs.data))

"""
    creepcomp(m[, t, params])

provide access to the creep compliance (J) of a given model m. It can be used with a broad range of inputs.

* When a time value is provided (or an array of time values), the function returns the corresponding value(s) of the creep compliance.

    creepcomp(m::RheoModel, time (single value or array))

    creepcomp(m::RheoModelClass, time (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    creepcomp(Maxwell, 1, k=1., η=1)

    creepcomp(Maxwell, [1,2,3], [1,1])

* When no time value is provided, creepcomp returns the creep compliance function itself.

    creepcomp(m::RheoModel)

    creepcomp(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    J = creepcomp(m)
    J([0,1,2])
"""
function creepcomp(m::RheoModel, t::Number)
    m._J(t)
end

function creepcomp(m::RheoModel, ta::Vector{T}) where T <: Number
    m._Ja(ta)
end

creepcomp(m::RheoModel) = x -> creepcomp(m,x)

function creepcomp(m::RheoModelClass, t::Number, params::Vector{T}) where T <: Number
    m._J(t, params)
end

function creepcomp(m::RheoModelClass, ta::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    m._Ja(ta, params)
end

function creepcomp(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.freeparams, err_string = "relaxmod")
    creepcomp(m, x, p)
end

function creepcomp(m::RheoModelClass, x; kwargs...)
    creepcomp(m, x, symbol_to_unicode(kwargs.data))
end

creepcomp(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> creepcomp(m, x, params)
creepcomp(m::RheoModelClass, params::NamedTuple) =  x -> creepcomp(m, x, params)
creepcomp(m::RheoModelClass; kwargs...) =  x -> creepcomp(m, x, symbol_to_unicode(kwargs.data))

"""
    storagemod(m[, ω, params])

provide access to the storage modulus (G') of a given model m.

* When a frequency value is provided (or an array of frequency values), the function returns the corresponding value(s) of the storage modulus.

    storagemod(m::RheoModel, frequency (single value or array))

    storagemod(m::RheoModelClass, frequency (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    storagemod(Maxwell, 1, k=1., η=1)

    storagemod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, storagemod returns the storage modulus function itself.

    storagemod(m::RheoModel)

    storagemod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    Gp = storagemod(m)
    Gp([0,1,2])
"""
function storagemod(m::RheoModel, ω::Number)
    m._Gp(ω)
end

function storagemod(m::RheoModel, ωa::Vector{T}) where T <: Number
    m._Gpa(ωa)
end

storagemod(m::RheoModel) = x -> storagemod(m,x)

function storagemod(m::RheoModelClass, ω::Number, params::Vector{T}) where T <: Number
    m._Gp(ω, params)
end

function storagemod(m::RheoModelClass, ωa::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
      m._Gpa(ωa, params)
end

function storagemod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.freeparams, err_string = "relaxmod")
    storagemod(m, x, p)
end

function storagemod(m::RheoModelClass, x; kwargs...)
    storagemod(m, x, symbol_to_unicode(kwargs.data))
end

storagemod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> storagemod(m, x, params)
storagemod(m::RheoModelClass, params::NamedTuple) =  x -> storagemod(m, x, params)
storagemod(m::RheoModelClass; kwargs...) =  x -> storagemod(m, x, symbol_to_unicode(kwargs.data))

"""
    lossmod(m[, ω, params])

provide access to the loss modulus (G'') of a given model m.

* When a frequency value is provided (or an array of frequency values), the function returns the corresponding value(s) of the loss modulus.

    lossmod(m::RheoModel, frequency (single value or array))

    lossmod(m::RheoModelClass, frequency (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    lossmod(Maxwell, 1, k=1., η=1)

    lossmod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, lossmod returns the loss modulus function itself.

    lossmod(m::RheoModel)

    lossmod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    Gpp = lossmod(m)
    Gpp([0,1,2])
"""
function lossmod(m::RheoModel, ω::Number)
    m._Gpp(ω)
end

function lossmod(m::RheoModel, ωa::Vector{T}) where T <: Number
    m._Gppa(ωa)
end

lossmod(m::RheoModel) = x -> lossmod(m,x)

function lossmod(m::RheoModelClass, ω::Number, params::Vector{T}) where T <: Number
    m._Gpp(ω, params)
end

function lossmod(m::RheoModelClass, ωa::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    m._Gppa(ωa, params)
end

function lossmod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.freeparams, err_string = "relaxmod")
    lossmod(m, x, p)
end

function lossmod(m::RheoModelClass, x; kwargs...)
    lossmod(m, x, symbol_to_unicode(kwargs.data))
end

lossmod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> lossmod(m, x, params)
lossmod(m::RheoModelClass, params::NamedTuple) =  x -> lossmod(m, x, params)
lossmod(m::RheoModelClass; kwargs...) =  x -> lossmod(m, x, symbol_to_unicode(kwargs.data))

"""
    dynamicmod(m[, ω, params])

provide access to the complex dynamic modulus (G' + i G'') of a given model m.

* When a frequency value is provided (or an array of frequency values), the function returns the corresponding value(s) of the complex dynamic modulus.

    dynamicmod(m::RheoModel, frequency (single value or array))

    dynamicmod(m::RheoModelClass, frequency (single value or array), parameters (array, named tupple or keyword parameters))

Examples:

    dynamicmod(Maxwell, 1, k=1., η=1)

    dynamicmod(Maxwell, [1,2,3], [1,1])

* When no time value is provided, lossmod returns the complex dynamic modulus function itself.

    dynamicmod(m::RheoModel)

    dynamicmod(m::RheoModelClass, parameters (array, named tupple or keyword parameters))

Examples:

    m = RheoModel(Maxwell, k=1., η=1)
    Gpp = dynamicmod(m)
    Gpp([0,1,2])

Note: use abs() and angle() to get the magnitude and phase of the complex modulus.
"""
function dynamicmod(m::RheoModel, ω::Number)
    m._Gp(ω) + m._Gpp(ω) * im
end

function dynamicmod(m::RheoModel, ωa::Vector{T}) where T <: Number
    m._Gpa(ωa) + m._Gppa(ωa) * im
end

dynamicmod(m::RheoModel) = x -> dynamicmod(m,x)

function dynamicmod(m::RheoModelClass, ω::Number, params::Vector{T}) where T <: Number
    m._Gp(ω, params) + m._Gpp(ω, params) * im
end

function dynamicmod(m::RheoModelClass, ωa::Vector{T1}, params::Vector{T2}) where {T1 <: Number, T2 <: Number}
    m._Gpa(ωa, params) + m._Gppa(ωa, params) * im
end

function dynamicmod(m::RheoModelClass, x, params::NamedTuple)
    # check all parameters are provided and create a well ordered named tuple
    p = check_and_reorder_parameters(params, m.freeparams, err_string = "relaxmod")
    dynamicmod(m, x, p)
end

function dynamicmod(m::RheoModelClass, x; kwargs...)
    dynamicmod(m, x, symbol_to_unicode(kwargs.data))
end

dynamicmod(m::RheoModelClass, params::Vector{T}) where T <: Number =  x -> dynamicmod(m, x, params)
dynamicmod(m::RheoModelClass, params::NamedTuple) =  x -> dynamicmod(m, x, params)
dynamicmod(m::RheoModelClass; kwargs...) =  x -> dynamicmod(m, x, symbol_to_unicode(kwargs.data))


text:
##############
# Owner Type #
##############

type Owner <: GitHubType
    typ::Nullable{GitHubString}
    email::Nullable{GitHubString}
    name::Nullable{GitHubString}
    login::Nullable{GitHubString}
    bio::Nullable{GitHubString}
    company::Nullable{GitHubString}
    location::Nullable{GitHubString}
    gravatar_id::Nullable{GitHubString}
    id::Nullable{Int}
    public_repos::Nullable{Int}
    owned_private_repos::Nullable{Int}
    total_private_repos::Nullable{Int}
    public_gists::Nullable{Int}
    private_gists::Nullable{Int}
    followers::Nullable{Int}
    following::Nullable{Int}
    collaborators::Nullable{Int}
    blog::Nullable{HttpCommon.URI}
    url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    updated_at::Nullable{Dates.DateTime}
    created_at::Nullable{Dates.DateTime}
    date::Nullable{Dates.DateTime}
    hireable::Nullable{Bool}
    site_admin::Nullable{Bool}
    avatar_url::Nullable{HttpCommon.URI}
    followers_url::Nullable{HttpCommon.URI}
    following_url::Nullable{HttpCommon.URI}
    gists_url::Nullable{HttpCommon.URI}
    starred_url::Nullable{HttpCommon.URI}
    subscriptions_url::Nullable{HttpCommon.URI}
    organizations_url::Nullable{HttpCommon.URI}
    repos_url::Nullable{HttpCommon.URI}
    events_url::Nullable{HttpCommon.URI}
    received_events_url::Nullable{HttpCommon.URI}
end

Owner(data::Dict) = json2github(Owner, data)
Owner(login::AbstractString, isorg = false) = Owner(Dict("login" => login, "typ" => isorg ? "User" : "Organization"))

namefield(owner::Owner) = owner.login

typprefix(isorg) = isorg ? "orgs" : "users"

#############
# Owner API #
#############

isorg(owner::Owner) = get(owner.typ, "") == "Organization"

owner(owner_obj::Owner; options...) = owner(name(owner_obj), isorg(owner_obj); options...)

function owner(owner_obj, isorg = false; options...)
    result = gh_get_json("/$(typprefix(isorg))/$(name(owner_obj))"; options...)
    return Owner(result)
end

function orgs(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/orgs"; options...)
    return map(Owner, results), page_data
end

function followers(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/followers"; options...)
    return map(Owner, results), page_data
end

function following(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/following"; options...)
    return map(Owner, results), page_data
end

repos(owner::Owner; options...) = repos(name(owner), isorg(owner); options...)

function repos(owner, isorg = false; options...)
    results, page_data = gh_get_paged_json("/$(typprefix(isorg))/$(name(owner))/repos"; options...)
    return map(Repo, results), page_data
end

ext: jl


text:
using Plots

rotation(x,y,theta)= [cos(theta) -sin(theta); sin(theta) cos(theta)]*[x;y]
points=[[0.0;0.0],[1.0;0.0]]
tempo=(points[2]-points[1])/sqrt(2)
point=rotation(tempo[1], tempo[2], pi/4)+points[1]
insert!(points,2, point)

for step in 1:20
    i=1
    while points[i]!= points[end]
        temp=(points[i+1]-points[i])/sqrt(2)
        point1=rotation(temp[1], temp[2], pi/4)+points[i]
        insert!(points, i+1, point1)
        point2=rotation(temp[1], temp[2], -pi/4)+points[i+2]
        insert!(points,i+3, point2)
        i+=4
    end
end
plot(hcat(points...)[1,:] , hcat(points...)[2,:], color=:darkred, legend=false)
savefig("D:\\Terms\\FALL00\\شبیه سازی رایانه ای\\report\\Levy_dragon.png")


text:
using ProceduralNoise
using Test

@testset "ProceduralNoise.jl" begin
    perlin((512, 512), (32, 32))
end


text:
### A Pluto.jl notebook ###
# v0.16.2

using Markdown
using InteractiveUtils

# ╔═╡ c30e6a6e-342c-11ec-0eb5-a553322ce2e7
begin
	using Images
	using LinearAlgebra
	using Plots
	using ImageIO

end

# ╔═╡ Cell order:
# ╠═c30e6a6e-342c-11ec-0eb5-a553322ce2e7


text:
"""Function that corrects chronologies for sudden jumps in time

Some occurrences in the model results can lead the CumDY function
to detect extra year transitions, resulting in sudden jumps in
the shell chronology or a start of the chronology at an age
beyond 1 year. This function removes these sharp transitions
and late onset by adding or subtracting whole years to the age
result.

[resultarray] Array containing the full results of
the optimized growth model

[T_per] The period length of one year (in days)

[agecorrection] Correct for jumps in age (/code{TRUE}) or
only for starting time (/code{FALSE})

[plot] Should the results be plotted? (/code{TRUE/FALSE})

Returns an updated and corrected version of \code{resultarray}

Examples

"""
module age_corr

export age_corr!

using DataFrames

function age_corr!(resultarray, T_per=365, plot=true,  agecorrection=true)
    mean_window_age = DataFrame(
                            1:(length(resultarray[1,,1])-5)
                        )
end

end


text:
"""
      optimize(
            f::Function, # objective function
            bounds::AbstractMatrix,
            method::AbstractAlgorithm = ECA();
            logger::Function = (status) -> nothing,
      )

Minimize a n-dimensional function `f` with domain `bounds` (2×n matrix) using `method = ECA()` by default.

# Example
Minimize f(x) = Σx² where x ∈ [-10, 10]³.

Solution:

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> bounds = [  -10.0 -10 -10; # lower bounds
                    10.0  10 10 ] # upper bounds
2×3 Array{Float64,2}:
 -10.0  -10.0  -10.0
  10.0   10.0   10.0

julia> result = optimize(f, bounds)
+=========== RESULT ==========+
  iteration: 1429
    minimum: 2.5354499999999998e-222
  minimizer: [-1.5135301653303966e-111, 3.8688354844737692e-112, 3.082095708730726e-112]
    f calls: 29989
 total time: 0.1543 s
+============================+
```
"""
function optimize(
        f::Function, # objective function
        bounds::AbstractMatrix,
        method::AbstractAlgorithm = ECA();
        logger::Function = (status) -> nothing,
    )

    #####################################
    # common methods
    #####################################
    # status = method.status
    information = method.information
    options = method.options
    parameters = method.parameters
    ###################################

    problem = Problem(f, Array(bounds); parallel_evaluation=options.parallel_evaluation)
    seed!(options.seed)

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

    start_time = time()

    status = method.status
    options.debug && @info("Initializing population...")
    status = initialize!(status,parameters, problem, information, options)
    method.status = status

    show_status(status, parameters, options)

    status.iteration = 1

    status.start_time = start_time
    convergence = State{typeof(status.best_sol)}[]

    ###################################
    # store convergence
    ###################################
    if options.store_convergence
        update_convergence!(convergence, status)
    end

    options.debug && @info("Starting main loop...")

    logger(status)

    while !status.stop
        status.iteration += 1

        update_state!(status, parameters, problem, information, options)

        # store the number of fuction evaluations
        status.f_calls = problem.f_calls

        show_status(status, parameters, options)

        if options.store_convergence
            update_convergence!(convergence, status)
        end

        status.overall_time = time() - status.start_time
        logger(status)

        # common stop criteria
        status.stop = status.stop ||
        call_limit_stop_check(status, information, options) ||
        iteration_stop_check(status, information, options)  ||
        time_stop_check(status, information, options)

        # user defined stop criteria
        stop_criteria!(status, parameters, problem, information, options)

    end

    status.overall_time = time() - status.start_time

    final_stage!(
                 status,
                 parameters,
                 problem,
                 information,
                 options
                )

    status.convergence = convergence

    return status

end

function show_status(status, parameters, options)
    !options.debug && (return)
    status.final_time = time()
    msg = "Current Status of " * string(typeof(parameters))
    @info msg
    display(status)
end


text:
function get_target_pv(obj::GroundPointing, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getTargetPV", TimeStampedFieldPVCoordinates, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_target_pv(obj::GroundPointing, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getTargetPV", TimeStampedPVCoordinates, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end


text:
module test_22_MatrixOfPlots

using  ModiaPlot

include("test_base/test_22_MatrixOfPlots.jl")

end


text:
part1(input) = nothing

function part2(input)
    medicine = last(split(readchomp(input), "\n\n"))
    c = filter(isuppercase, collect(medicine))
    return length(c) - 2 * (count(==('R'), c) + count(==('Y'), c)) - 1
end


text:
## MutablePaddedViews.jl

struct MutablePaddedView{T,N,I,A} <: AbstractArray{T,N}
    data::PaddedView{T,N,I,A}
end

Base.size(A::MutablePaddedView) = size(A.data)
# Base.ndims(A::MutablePaddedView) = ndims(A.data)

# The change below makes select_region into a writable region
Base.@propagate_inbounds function Base.setindex!(A::MutablePaddedView{T,N}, v, pos::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A.data, pos...)
    if Base.checkbounds(Bool, A.data.data, pos...)
        return setindex!(A.data.data, v, pos... )
    else
        return A.data.fillvalue
    end
end

Base.@propagate_inbounds function Base.getindex(A::MutablePaddedView{T,N}, pos::Vararg{Int,N}) where {T,N}
    getindex(A.data, pos...)
end


text:
# ---
# title: 393. UTF-8 Validation
# id: problem393
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/utf-8-validation/description/>
# hidden: true
# ---
#
# A character in UTF8 can be from **1 to 4 bytes** long, subjected to the
# following rules:
#
#   1. For 1-byte character, the first bit is a 0, followed by its unicode code.
#   2. For n-bytes character, the first n-bits are all one's, the n+1 bit is 0, followed by n-1 bytes with most significant 2 bits being 10.
#
# This is how the UTF-8 encoding would work:
#
#
#
#        Char. number range  |        UTF-8 octet sequence
#           (hexadecimal)    |              (binary)
#        --------------------+---------------------------------------------
#        0000 0000-0000 007F | 0xxxxxxx
#        0000 0080-0000 07FF | 110xxxxx 10xxxxxx
#        0000 0800-0000 FFFF | 1110xxxx 10xxxxxx 10xxxxxx
#        0001 0000-0010 FFFF | 11110xxx 10xxxxxx 10xxxxxx 10xxxxxx
#
#
# Given an array of integers representing the data, return whether it is a valid
# utf-8 encoding.
#
# **Note:**
# The input is an array of integers. Only the **least significant 8 bits** of
# each integer is used to store the data. This means each integer represents
# only 1 byte of data.
#
# **Example 1:**
#
#
#
#     data = [197, 130, 1], which represents the octet sequence: **11000101 10000010 00000001**.
#
#     Return **true**.
#     It is a valid utf-8 encoding for a 2-bytes character followed by a 1-byte character.
#
#
# **Example 2:**
#
#
#
#     data = [235, 140, 4], which represented the octet sequence: **11101011 10001100 00000100**.
#
#     Return **false**.
#     The first 3 bits are all one's and the 4th bit is 0 means it is a 3-bytes character.
#     The next byte is a continuation byte which starts with 10 and that's correct.
#     But the second continuation byte does not start with 10, so it is invalid.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Base: show, getindex, iterate, HasLength, HasEltype, length

struct ConditionalIndependency
  f::SymbolList
  s::SymbolList
  c::SymbolListOrNothing
end

ConditionalIndependency(f::SymbolList, s::SymbolList) =
  ConditionalIndependency(f, s, nothing)

function ConditionalIndependency(a::Array{Symbol,1})
  @assert length(a) > 1 "Length Symbol set too short."
  if length(a) == 2
    return ConditionalIndependency(a[1], a[2])
  end
  ConditionalIndependency(a[1], a[2], a[3:end])
end

import Base.print
function print(io::IO, x::BasisSet)
  str = ""
  str *= "BasisSet["
  for c in bs.bs
    ci = ConditionalIndependency(c)
    if isnothing(ci.c)
      str *= "  :$(ci.f) \u2210 :$(ci.s)"
    else
      str *= "  :$(ci.f) \u2210 :$(ci.s) | $(ci.c)"
    end
  end
  str *= "  ]"
  str
end

function pluto_string(bs::BasisSet)
  str = ""
  str *= "BasisSet["
  for c in bs.bs
    ci = ConditionalIndependency(c)
    if isnothing(ci.c)
      str *= "  :$(ci.f) \u2210 :$(ci.s)"
    else
      str *= "  :$(ci.f) \u2210 :$(ci.s) | $(ci.c)"
    end
  end
  str *= "  ]"
  str
end

function ci_show(io::IO, ci::ConditionalIndependency)
  if isnothing(ci.c)
    println("  :$(ci.f) \u2210 :$(ci.s)")
  else
    println("  :$(ci.f) \u2210 :$(ci.s) | $(ci.c)")
   end
end

show(io::IO, ci::ConditionalIndependency) = ci_show(io, ci)

export
  ConditionalIndependency,
  print,
  pluto_string


text:
macro replace_1_with_x(expr)
   esc(replace_1(expr))
end

replace_1(atom) = atom == 1 ? :x : atom
replace_1(e::Expr) =
    Expr(e.head, map(replace_1, e.args)...)

max_stars_repo_name: ninjaaron/dispatching-design-patterns

max_issues_repo_name: ninjaaron/dispatching-design-patterns

max_forks_repo_name: ninjaaron/dispatching-design-patterns


text:
# https://github.com/jaypmorgan/Adversarial.jl.git
#
# I'm specifically not using PGD, but myPGD instead.
# using Adversarial: FGSM
# import Adversarial

using BSON: @save, @load

# CAUTION: A bug in CuArrays, https://github.com/FluxML/Flux.jl/issues/839
#
# onecold(cpu([1,-2])) # => 1
# onecold(gpu([1,-2])) # => 2
#
# The only scalar operation is the "mapreduce" in onecold. It also has bug that
# causes argmax to give incorrect results, thus, I'm disablng scalar operations,
# and force using CPU to compute onecold. CAUTION I should also check whether
# argmax bug affects my code elsewhere.
CuArrays.allowscalar(false)
# FIXME do I need .data for y?
# FIXME logits.data?
accuracy_with_logits(logits, y) = mean(onecold(cpu(logits)) .== onecold(cpu(y)))
# CuArrays.allowscalar(true)
# accuracy_with_logits(logits, y) = mean(onecold(logits) .== onecold(y))

my_xent(logits, y) = Flux.logitcrossentropy(logits, y)

"""Not using, the same as Adversarial.jl's FGSM.
"""
function myFGSM(model, x, y; ϵ = 0.3, clamp_range = (0, 1))
    # FIXME what if I just don't create Flux.params(model)?
    px, θ = Flux.param(x), Flux.params(model)
    Flux.Tracker.gradient(() -> my_xent(model(px), y), θ)
    x_adv = clamp.(x + (Float32(ϵ) * sign.(px.grad)), clamp_range...)
end

"""
My modifications:

- remove (δ < ϵ) condition. This seems to fix PGD performance problem. This
  condition also has potential problem when attacking in batch: you should
  really not stop when one image in the batch reached the epsilon-ball

- clipped into valid range of epsilon, i.e. clip_eta in cleverhans
"""
function myPGD(model, x, y;
               ϵ = 0.3, step_size = 0.01,
               iters = 40, clamp_range = (0, 1))
    # start from the random point
    # eta = gpu(randn(Float32, size(x)...))
    # eta = clamp.(eta, -ϵ, ϵ)
    # DEBUG uniform sample
    eta = (gpu(rand(Float32, size(x)...)) .- 0.5) * 2 * ϵ
    # FIXME this clamp should not be necessary
    eta = clamp.(eta, -ϵ, ϵ)
    x_adv = x + Float32.(eta)
    x_adv = clamp.(x_adv, clamp_range...)
    # x_adv = clamp.(x + (r * Float32(step_size)), clamp_range...)
    for iter = 1:iters
        x_adv = myFGSM(model, x_adv, y;
                       ϵ = step_size, clamp_range = clamp_range)
        eta = x_adv - x
        eta = clamp.(eta, -ϵ, ϵ)
        x_adv = x + Float32.(eta)
        x_adv = clamp.(x_adv, clamp_range...)
    end
    return x_adv
end

function attack_FGSM(model, x, y)
    x_adv = myFGSM(model, x, y; ϵ = 0.3)
end

function attack_PGD(model, x, y)
    x_adv = myPGD(model, x, y;
                  ϵ = 0.3,
                  step_size = 0.01,
                  iters = 40)
end

function attack_PGD_k(k)
    (model, x, y) -> begin
        x_adv = myPGD(model, x, y;
                      ϵ = 0.3,
                      step_size = 0.01,
                      iters = k)
    end
end

function attack_CIFAR10_PGD_k(k)
    # FIXME the model contains BN layers, and is set to test mode during
    # testing. However, I still need to perform attack during testing, which
    # requires gradient. But testing mode model cannot take gradient?
    function attack(model, x, y)
        # cifar: 8.0/255, 7, 2./255
        myPGD(model, x, y;
              ϵ = 8/255,
              step_size = 2/255,
              iters = k)
    end
end

# something like tf.keras.metrics.Mean
# a good reference https://github.com/apache/incubator-mxnet/blob/master/julia/src/metric.jl
mutable struct MeanMetric
    sum::Float64
    n::Int
    MeanMetric() = new(0.0, 0)
end
function add!(m::MeanMetric, v)
    m.sum += v
    m.n += 1
end
get(m::MeanMetric) = m.sum / m.n
function get!(m::MeanMetric)
    res = m.sum / m.n
    reset!(m)
    res
end
function reset!(m::MeanMetric)
    m.sum = 0.0
    m.n = 0
end

function train!(model, opt, ds;
                train_steps=ds.nbatch,
                from_steps=1,
                print_steps=100,
                augment=(x)->x,
                logger=nothing,
                save_cb=(i)->nothing,
                test_cb=(i)->nothing)
    ps=Flux.params(model)

    loss_metric = MeanMetric()
    acc_metric = MeanMetric()

    @info "Training for $train_steps steps, printing every $print_steps steps .."
    @showprogress 0.1 "Training..." for step in from_steps:train_steps
        x, y = next_batch!(ds)
        # TODO add data augment here
        x = augment(x) |> gpu
        y = gpu(y)

        gs = Flux.Tracker.gradient(ps) do
            logits = model(x)
            loss = my_xent(logits, y)
            add!(loss_metric, loss.data)
            add!(acc_metric, accuracy_with_logits(logits.data, y))
            loss
        end
        Flux.Tracker.update!(opt, ps, gs)

        if step % print_steps == 0
            println()
            loss=get!(loss_metric)
            acc=get!(acc_metric)
            @info "data"  loss acc
            if typeof(logger) <: TBLogger
                log_value(logger, "loss", loss, step=step)
                log_value(logger, "acc", acc, step=step)
            end
        end
        test_cb(step)
        save_cb(step)
    end
end

function advtrain!(model, opt, attack_fn, ds;
                   train_steps=ds.nbatch,
                   from_steps=1,
                   print_steps=50,
                   logger=nothing,
                   λ = 0,
                   save_cb=(i)->nothing,
                   test_cb=(i)->nothing)
    ps=Flux.params(model)

    m_cleanloss = MeanMetric()
    m_cleanacc = MeanMetric()
    m_advloss = MeanMetric()
    m_advacc = MeanMetric()
    @info "Training for $train_steps steps, printing every $print_steps steps .."
    @showprogress 0.1 "Training..." for step in from_steps:train_steps
        x, y = next_batch!(ds) |> gpu
        # this computation won't affect model parameter gradients
        x_adv = attack_fn(model, x, y)

        gs = Flux.Tracker.gradient(ps) do
            # FIXME this will slow down the model twice
            adv_logits = model(x_adv)
            adv_loss = my_xent(adv_logits, y)

            # TODO add l2 regularizers
            #
            # l1(x) = sum(abs.(x))
            # l2(x) = sum(x.^2)
            # l2_decay = sum(l2, weight_params(model))

            clean_logits = model(x)
            clean_loss = my_xent(clean_logits, y)

            add!(m_advloss, adv_loss.data)
            add!(m_advacc, accuracy_with_logits(adv_logits.data, y))
            add!(m_cleanloss, clean_loss.data)
            add!(m_cleanacc, accuracy_with_logits(clean_logits.data, y))

            l = adv_loss + λ * clean_loss
        end
        Flux.Tracker.update!(opt, ps, gs)

        if step % print_steps == 0
            println()
            natloss = get!(m_cleanloss)
            advloss = get!(m_advloss)
            natacc = get!(m_cleanacc)
            advacc = get!(m_advacc)
            @info "data" natloss advloss natacc advacc
            # TODO log training time
            if typeof(logger) <: TBLogger
                log_value(logger, "loss/nat_loss", natloss, step=step)
                log_value(logger, "loss/adv_loss", advloss, step=step)
                log_value(logger, "acc/nat_acc", natacc, step=step)
                log_value(logger, "acc/adv_acc", advacc, step=step)
            end
        end
        test_cb(step)
        save_cb(step)
    end
end

function free_train!(model, opt, ds, ε;
                     train_steps=ds.nbatch,
                     from_steps=1,
                     print_steps=50,
                     logger=nothing,
                     λ = 0,
                     replay_m=8,
                     save_cb=(i)->nothing,
                     test_cb=(i)->nothing)
    ps=Flux.params(model)

    m_cleanloss = MeanMetric()
    m_cleanacc = MeanMetric()
    m_advloss = MeanMetric()
    m_advacc = MeanMetric()
    @info "Training for $train_steps steps, printing every $print_steps steps .."
    @showprogress 0.1 "Training..." for step in from_steps:train_steps
        x, y = next_batch!(ds) |> gpu

        # comment out to use "warm up"
        δ = gpu(zeros(size(x)...))
        # FIXME which is one step?
        for i in 1:replay_m
            xx = x + δ
            px = Flux.param(xx)
            gs = Flux.Tracker.gradient(ps) do
                logits = model(px)
                loss = my_xent(logits, y)
                # record statistics
                add!(m_advloss, loss.data)
                add!(m_advacc, accuracy_with_logits(logits.data, y))
                loss
            end
            Flux.Tracker.update!(opt, ps, gs)
            δ += Float32(ε) * sign.(px.grad)
            δ = clamp.(δ, -Float32(ε), Float32(ε))
        end

        if step % print_steps == 0
            println()
            natloss = get!(m_cleanloss)
            advloss = get!(m_advloss)
            natacc = get!(m_cleanacc)
            advacc = get!(m_advacc)
            @info "data" natloss advloss natacc advacc
            # TODO log training time
            if typeof(logger) <: TBLogger
                log_value(logger, "loss/nat_loss", natloss, step=step)
                log_value(logger, "loss/adv_loss", advloss, step=step)
                log_value(logger, "acc/nat_acc", natacc, step=step)
                log_value(logger, "acc/adv_acc", advacc, step=step)
            end
        end
        test_cb(step)
        save_cb(step)
    end
end

function create_save_cb(model_file, model; save_steps)
    function save_cb(step)
        if step % save_steps == 0
            @info "saving .."
            # FIXME opt cannot be saved
            # FIXME logger cannot be saved
            #
            # FIXME sometimes user interrupt would crash the file, I'm saving it
            # as xxx-tmp, and perform move afterwards, hopefully fix this.
            @time @save model_file*".tmp" model=cpu(model) from_steps=step
            mv(model_file*".tmp", model_file, force=true)
        end
    end
end

function create_advae_test_cb(ae, cnn, test_ds; kwargs...)
    testmode_fn() = Flux.testmode!(ae)
    testmode_exit_fn() = Flux.testmode!(ae, false)
    create_adv_test_cb_impl(Chain(ae, cnn), test_ds,
                            testmode_fn, testmode_exit_fn;
                            kwargs...)
end

function create_adv_test_cb(model, test_ds; kwargs...)
    testmode_fn() = Flux.testmode!(model)
    testmode_exit_fn() = Flux.testmode!(model, false)
    create_adv_test_cb_impl(model, test_ds,
                            testmode_fn, testmode_exit_fn;
                            kwargs...)
end

function create_adv_test_cb_impl(model, test_ds, testmode_fn, testmode_exit_fn;
                                 test_per_steps, test_run_steps, attack_fn, logger=nothing)
    function test_cb(step)
        if step % test_per_steps == 0
            println()
            @info "testing for $test_run_steps steps .."
            m_cleanloss = MeanMetric()
            m_cleanacc = MeanMetric()
            m_advloss = MeanMetric()
            m_advacc = MeanMetric()

            testmode_fn()
            @showprogress 0.1 "Inner testing..." for i in 1:test_run_steps
                x, y = next_batch!(test_ds) |> gpu
                # this computation won't affect model parameter gradients
                x_adv = attack_fn(model, x, y)
                adv_logits = model(x_adv)
                adv_loss = my_xent(adv_logits, y)
                clean_logits = model(x)
                clean_loss = my_xent(clean_logits, y)
                add!(m_advloss, adv_loss.data)
                add!(m_advacc, accuracy_with_logits(adv_logits.data, y))
                add!(m_cleanloss, clean_loss.data)
                add!(m_cleanacc, accuracy_with_logits(clean_logits.data, y))
            end
            testmode_exit_fn()

            natloss = get!(m_cleanloss)
            advloss = get!(m_advloss)
            natacc = get!(m_cleanacc)
            advacc = get!(m_advacc)
            @info "test data" natloss advloss natacc advacc

            # use explicit API to avoid manipulating log_step_increment
            if typeof(logger) <: TBLogger
                @info "logging results .."
                log_value(logger, "loss/nat_loss", natloss, step=step)
                log_value(logger, "loss/adv_loss", advloss, step=step)
                log_value(logger, "acc/nat_acc", natacc, step=step)
                log_value(logger, "acc/adv_acc", advacc, step=step)
            end
        end
    end
end

function create_test_cb(model, test_ds; logger=nothing)
    function test_cb(step)
        test_per_steps = 100
        test_run_steps = 20

        if step % test_per_steps == 0
            println()
            @info "testing for $test_run_steps steps .."
            m_cleanloss = MeanMetric()
            m_cleanacc = MeanMetric()

            # into testing mode
            Flux.testmode!(model)

            @showprogress 0.1 "Inner testing..." for i in 1:test_run_steps
                x, y = next_batch!(test_ds) |> gpu
                # this computation won't affect model parameter gradients
                logits = model(x)
                loss = my_xent(logits, y)
                add!(m_cleanloss, loss.data)
                add!(m_cleanacc, accuracy_with_logits(logits.data, y))
            end

            # back into training node
            Flux.testmode!(model, false)

            loss = get!(m_cleanloss)
            acc = get!(m_cleanacc)
            @info "test data" loss acc

            # use explicit API to avoid manipulating log_step_increment
            if typeof(logger) <: TBLogger
                @info "logging results .."
                log_value(logger, "loss", loss, step=step)
                log_value(logger, "acc", acc, step=step)
            end
        end
    end
end

##############################
## AE train
##############################

my_mse(ŷ, y) = sum((ŷ .- y).*(ŷ .- y)) * 1 // length(y)

function aetrain!(model, opt, ds;
                  train_steps=ds.nbatch,
                  from_steps=1,
                  print_steps=50)
    ps=Flux.params(model)

    loss_metric = MeanMetric()

    @info "Training for $train_steps steps, printing even $print_steps steps .."
    @showprogress 0.1 "Training..." for step in from_steps:train_steps
        x, _ = next_batch!(ds) |> gpu
        gs = Flux.Tracker.gradient(ps) do
            logits = model(x)
            loss = my_mse(logits, x)
            add!(loss_metric, loss.data)
            loss
        end
        Flux.Tracker.update!(opt, ps, gs)

        if step % print_steps == 0
            println()
            @info "data" get!(loss_metric)
        end
    end
end

function advae_train!(ae, cnn, opt, attack_fn, ds;
                      train_steps=ds.nbatch,
                      from_steps=1,
                      β=1,
                      λ=0,
                      γ=0,
                      print_steps=50,
                      logger=nothing,
                      save_cb=(i)->nothing,
                      test_cb=(i)->nothing)
    # update only ae parameter
    ps=Flux.params(ae)

    m_cleanloss = MeanMetric()
    m_cleanacc = MeanMetric()
    m_advloss = MeanMetric()
    m_advacc = MeanMetric()
    m_recloss = MeanMetric()

    full_model = Chain(ae, cnn)

    @info "Training for $train_steps steps, printing every $print_steps steps .."
    @showprogress 0.1 "Training..." for step in from_steps:train_steps
        x, y = next_batch!(ds) |> gpu
        # this computation won't affect model parameter gradients
        x_adv = attack_fn(full_model, x, y)

        gs = Flux.Tracker.gradient(ps) do
            adv_logits = full_model(x_adv)
            adv_loss = my_xent(adv_logits, y)

            clean_logits = full_model(x)
            clean_loss = my_xent(clean_logits, y)

            rec_loss = my_mse(ae(x), x)

            add!(m_advloss, adv_loss.data)
            add!(m_advacc, accuracy_with_logits(adv_logits.data, y))
            add!(m_cleanloss, clean_loss.data)
            add!(m_cleanacc, accuracy_with_logits(clean_logits.data, y))
            add!(m_recloss, rec_loss.data)

            l = adv_loss + λ * clean_loss + γ * rec_loss
        end
        # FIXME clean up cnn parameter gradients?
        Flux.Tracker.update!(opt, ps, gs)

        if step % print_steps == 0
            println()
            natloss = get!(m_cleanloss)
            advloss = get!(m_advloss)
            recloss = get!(m_recloss)
            natacc = get!(m_cleanacc)
            advacc = get!(m_advacc)
            @info "data" natloss advloss natacc advacc recloss
            if typeof(logger) <: TBLogger
                log_value(logger, "loss/nat_loss", natloss, step=step)
                log_value(logger, "loss/adv_loss", advloss, step=step)
                log_value(logger, "loss/rec_loss", recloss, step=step)
                log_value(logger, "acc/nat_acc", natacc, step=step)
                log_value(logger, "acc/adv_acc", advacc, step=step)
            end
        end
        test_cb(step)
        save_cb(step)
    end
end

max_stars_repo_name: lihebi/AdvAE

max_issues_repo_name: lihebi/AdvAE

max_forks_repo_name: lihebi/AdvAE


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =========================================================================================
"""
    mutable struct TessOutput{T}
        result::Union{T, Nothing}
    end

Output object returned when adding an output renderer to a pipeline.  This object usually
holds `nothing` until the [`tess_run_pipeline`](@ref) method completes.

See slso: [`tess_pipeline_text`](@ref)
"""
mutable struct TessOutput{T}
    result::Union{T, Nothing}
    TessOutput(String) = new{String}(nothing)
end

# =========================================================================================
"""
    show(
        io::IO,
        inst::TessOutput
    )::Nothing

Display summary information about the TessOutput instance.

__Arguments:__

| T | Name  | Default | Description
|---| :---- | :------ | :----------
| R | io    |         | The stream to write the information to.
| R | inst  |         | The TessOutput instance to display info about.
"""
function Base.show(
            io::IO,
            inst::TessOutput
        )::Nothing
    if inst.result === nothing
        print(io, "No output available.")
    else
        show(io, inst.result)
    end
    nothing
end

# =========================================================================================
"""
    is_available(
        p::TessOutput
    )::Bool

Test if [`TessOutput`](@ref) object contains something yet.  Returns `false` if it contains
`nothing`.
"""
function is_available(
            p::TessOutput
        )::Bool
    return p.result !== nothing
end

# =========================================================================================
"""
    getindex(
        p::TessOutput{T}
    )::Union{T, Nothing} where T

Provides functionality to retrieve the contents of the  [`TessOutput`](@ref) object by
using [].
"""
function Base.getindex(
            p::TessOutput{T}
        )::Union{T, Nothing} where T
    return p.result
end


text:
function _convert(node::treeregressor.NodeMeta{S}, labels::Array{T}) where {S, T <: Float64}
    if node.is_leaf
        return Leaf{T}(node.label, labels[node.region])
    else
        left = _convert(node.l, labels)
        right = _convert(node.r, labels)
        return Node{S, T}(node.feature, node.threshold, left, right)
    end
end

function _convertOOB(
        node   :: treeregressor.NodeMeta{S},
        labels :: Array{T},
        inds   :: AbstractVector{Int},
        oob    :: AbstractMatrix{S}) where {S, T}

    tree = _convert(node, labels)
    return TreeOOB{S, T}(tree, inds, oob)
end

function build_tree(
        inds               :: AbstractVector{Int},
        labels             :: AbstractVector{T},
        features           :: AbstractMatrix{S},
        n_subfeatures       = 0,
        max_depth           = -1,
        min_samples_leaf    = 5,
        min_samples_split   = 2,
        min_purity_increase = 0.0;
        rng                 = Random.GLOBAL_RNG) where {S, T <: Float64}

    if max_depth == -1
        max_depth = typemax(Int)
    end
    if n_subfeatures == 0
        n_subfeatures = size(features, 2)
    end

    rng = mk_rng(rng)::Random.AbstractRNG
    t = treeregressor.fit(
        X                   = features[inds, :],
        Y                   = labels[inds],
        W                   = nothing,
        max_features        = Int(n_subfeatures),
        max_depth           = Int(max_depth),
        min_samples_leaf    = Int(min_samples_leaf),
        min_samples_split   = Int(min_samples_split),
        min_purity_increase = Float64(min_purity_increase),
        rng                 = rng)

    return _convertOOB(t.root, labels[t.labels], setdiff(collect(1:length(labels)), inds),
        features[setdiff(collect(1:length(labels)), inds), :])
end

function build_forest_oob(
        labels              :: AbstractVector{T},
        features            :: AbstractMatrix{S},
        n_subfeatures       = -1,
        n_trees             = 10,
        partial_sampling    = 0.7,
        max_depth           = -1,
        min_samples_leaf    = 5,
        min_samples_split   = 2,
        min_purity_increase = 0.0;
        rng                 = Random.GLOBAL_RNG) where {S, T <: Float64}

    if n_trees < 1
        throw("the number of trees must be >= 1")
    end
    if !(0.0 < partial_sampling <= 1.0)
        throw("partial_sampling must be in the range (0,1]")
    end

    if n_subfeatures == -1
        n_features = size(features, 2)
        n_subfeatures = round(Int, sqrt(n_features))
    end

    t_samples = length(labels)
    n_samples = floor(Int, partial_sampling * t_samples)

    forest = Vector{TreeOOB{S, T}}(undef, n_trees)

    if rng isa Random.AbstractRNG
        Threads.@threads for i in 1:n_trees
            inds = rand(rng, 1:t_samples, n_samples)
            forest[i] = build_tree(
                inds,
                labels,
                features,
                n_subfeatures,
                max_depth,
                min_samples_leaf,
                min_samples_split,
                min_purity_increase,
                rng = rng)
        end
    elseif rng isa Integer # each thread gets its own seeded rng
        Threads.@threads for i in 1:n_trees
            Random.seed!(rng + i)
            inds = rand(1:t_samples, n_samples)
            forest[i] = build_tree(
                inds,
                labels,
                features,
                n_subfeatures,
                max_depth,
                min_samples_leaf,
                min_samples_split,
                min_purity_increase)
        end
    else
        throw("rng must of be type Integer or Random.AbstractRNG")
    end

    return EnsembleOOB{S, T}(length(labels), features, forest)
end


text:
using Documenter, Rocket

makedocs(
    modules  = [ Rocket ],
    clean    = true,
    sitename = "Rocket.jl",
    pages    = [
        "Home"            => "index.md",
        "Getting started" => "getting-started.md",
        "Manual"     => [
            "Observable"   => "observables/about.md",
            "Actor"        => "actors/about.md",
            "Subscription" => "teardown/about.md",
            "Operator"     => "operators/about.md",
            "Subject"      => "subjects/about.md",
        ],
        "Library" => [
            "Observables" => [
                "Function"          => "observables/types/function.md",
                "Single"            => "observables/types/single.md",
                "Array"             => "observables/types/array.md",
                "Iterable"          => "observables/types/iterable.md",
                "Faulted"           => "observables/types/faulted.md",
                "Never"             => "observables/types/never.md",
                "Completed"         => "observables/types/completed.md",
                "Timer"             => "observables/types/timer.md",
                "Interval"          => "observables/types/interval.md",
                "Proxy"             => "observables/types/proxy.md",
                "File"              => "observables/types/file.md",
                "Combined"          => "observables/types/combined.md",
                "Combined Updates"  => "observables/types/combined_updates.md",
                "Collected"         => "observables/types/collected.md",
                "Race"              => "observables/types/race.md",
                "Connectable"       => "observables/types/connectable.md",
                "Merged"            => "observables/types/merged.md",
                "Concat"            => "observables/types/concat.md",
                "Generate"          => "observables/types/generate.md",
                "Network"           => "observables/types/network.md",
                "Defer"             => "observables/types/defer.md",
                "Zipped"            => "observables/types/zipped.md",
            ],
            "Actors" => [
                "Lambda"        => "actors/types/lambda.md",
                "Logger"        => "actors/types/logger.md",
                "Sync"          => "actors/types/sync.md",
                "Keep"          => "actors/types/keep.md",
                "Buffer"        => "actors/types/buffer.md",
                "Void"          => "actors/types/void.md",
                "Function"      => "actors/types/function.md",
                "Server"        => "actors/types/server.md",
            ],
            "Subscriptions" => [
                "Void"  => "teardown/types/void.md",
            ],
            "Operators" => [
                "All"            => "operators/all.md",
                "Creation"       => "operators/creation/about.md",
                "Transformation" => [
                    "About transformation operators" => "operators/transformation/about.md",
                    "map"                            => "operators/transformation/map.md",
                    "map_to"                         => "operators/transformation/map_to.md",
                    "scan"                           => "operators/transformation/scan.md",
                    "accumulated"                    => "operators/transformation/accumulated.md",
                    "enumerate"                      => "operators/transformation/enumerate.md",
                    "uppercase"                      => "operators/transformation/uppercase.md",
                    "lowercase"                      => "operators/transformation/lowercase.md",
                    "to_array"                       => "operators/transformation/to_array.md",
                    "switch_map"                     => "operators/transformation/switch_map.md",
                    "`switch_map_to`"                => "operators/transformation/switch_map_to.md",
                    "merge_map"                      => "operators/transformation/merge_map.md",
                    "concat_map"                     => "operators/transformation/concat_map.md",
                    "`concat_map_to`"                => "operators/transformation/concat_map_to.md",
                    "exhaust_map"                    => "operators/transformation/exhaust_map.md",
                    "start_with"                     => "operators/transformation/start_with.md",
                    "pairwise"                       => "operators/transformation/pairwise.md",
                ],
                "Filtering" => [
                    "About filtering operators" => "operators/filtering/about.md",
                    "filter"                    => "operators/filtering/filter.md",
                    "filter_type"               => "operators/filtering/filter_type.md",
                    "some"                      => "operators/filtering/some.md",
                    "take"                      => "operators/filtering/take.md",
                    "take_until"                => "operators/filtering/take_until.md",
                    "first"                     => "operators/filtering/first.md",
                    "last"                      => "operators/filtering/last.md",
                    "find"                      => "operators/filtering/find.md",
                    "find_index"                => "operators/filtering/find_index.md",
                    "ignore"                    => "operators/filtering/ignore.md",
                ],
                "Mathematical and Aggregate" => [
                    "About mathematical operators" => "operators/mathematical/about.md",
                    "count"                        => "operators/mathematical/count.md",
                    "max"                          => "operators/mathematical/max.md",
                    "min"                          => "operators/mathematical/min.md",
                    "reduce"                       => "operators/mathematical/reduce.md",
                    "sum"                          => "operators/mathematical/sum.md",
                ],
                "Error handling" => [
                    "About error handling operators" => "operators/errors/about.md",
                    "catch_error"                    => "operators/errors/catch_error.md",
                    "rerun"                          => "operators/errors/rerun.md"
                ],
                "Join" => [
                    "with_latest" => "operators/join/with_latest.md"
                ],
                "Multicasting" => [
                    "About multicasting operators" => "operators/multicasting/about.md",
                    "multicast"                    => "operators/multicasting/multicast.md",
                    "publish"                      => "operators/multicasting/publish.md",
                    "publish_behavior"             => "operators/multicasting/publish.md",
                    "publish_replay"               => "operators/multicasting/publish.md",
                    "share"                        => "operators/multicasting/share.md",
                    "share_replay"                 => "operators/multicasting/share.md",
                ],
                "Utility" => [
                    "About utility operators" => "operators/utility/about.md",
                    "tap"                     => "operators/utility/tap.md",
                    "`tap_on_subscribe`"      => "operators/utility/tap_on_subscribe.md",
                    "`tap_on_unsubscribe`"    => "operators/utility/tap_on_unsubscribe.md",
                    "`tap_on_complete`"       => "operators/utility/tap_on_complete.md",
                    "delay"                   => "operators/utility/delay.md",
                    "safe"                    => "operators/utility/safe.md",
                    "noop"                    => "operators/utility/noop.md",
                    "ref_count"               => "operators/utility/ref_count.md",
                    "async"                   => "operators/utility/async.md",
                    "`default_if_empty`"      => "operators/utility/default_if_empty.md",
                    "`error_if_empty`"        => "operators/utility/error_if_empty.md",
                    "skip_next"               => "operators/utility/skip_next.md",
                    "skip_error"              => "operators/utility/skip_error.md",
                    "skip_complete"           => "operators/utility/skip_complete.md",
                    "discontinue"             => "operators/utility/discontinue.md",
                    "limit_subscribers"       => "operators/utility/limit_subscribers.md",
                ]
            ],
            "Subjects" => [
                "Subject"  => "subjects/types/subject.md",
                "Behavior" => "subjects/types/behavior.md",
                "Replay"   => "subjects/types/replay.md",
                "Pending"  => "subjects/types/pending.md",
                "Recent"   => "subjects/types/recent.md",
            ]
        ],
        "API"          => [
            "Observables"  => "api/observables.md",
            "Actors"       => "api/actors.md",
            "Teardown"     => "api/teardown.md",
            "Operators"    => "api/operators.md",
            "Subjects"     => "api/subjects.md"
        ],
        "TODO"         => "todo.md",
        "Contributing" => "contributing.md",
        "Utils"        => "utils.md"
    ],
    format   = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

if get(ENV, "CI", nothing) == "true"
    deploydocs(
        repo = "github.com/biaslab/Rocket.jl.git"
    )
end


text:
function compute_L(file, species)
    w, (x,y,z), px, py, pz = file["weight/"*species,
                                  "grid/"*species,
                                  "px/"*species,
                                  "py/"*species,
                                  "pz/"*species]

    r = SVector{3}.(x,y,z)
    p = VectorVariable(px, py, pz)

    L = w .* r .× p
end

compute_L(dir::Symbol, file, species) = compute_L(Val(dir), file, species)

function compute_L(::Val{:x}, file, species)
    w, r, py, pz = file["weight/"*species,
                        "grid/"*species,
                        "py/"*species,
                        "pz/"*species]

    y = r[2]
    z = r[3]
    Lx = @. w * (y * pz - z * py)
end


text:
"""

The exported names are
$(EXPORTS)
"""
module CellData

using Test
using DocStringExtensions
using FillArrays

using Gridap.Helpers
using Gridap.Algebra
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry

using NearestNeighbors
using StaticArrays

import Gridap.Arrays: lazy_append
import Gridap.Arrays: get_array
import Gridap.Arrays: evaluate!
import Gridap.Arrays: return_cache
import Gridap.Fields: gradient
import Gridap.Fields: ∇∇
import Gridap.Fields: integrate
import Gridap.Fields: grad2curl
import Gridap.Geometry: num_cells
import Gridap.Geometry: get_triangulation

import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
import LinearAlgebra: det, tr, cross, dot, ⋅
import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj

export gradient, ∇
export ∇∇
export inner, ⊙, outer, ⊗, double_contraction, ⋅¹, ⋅², symmetric_part
export det, tr, cross, ×, dot, ⋅

export DomainStyle
export ReferenceDomain
export PhysicalDomain
export CellDatum
export get_data
export get_triangulation
export change_domain
export test_cell_datum
export CellPoint
export get_cell_points
export CellField
export jump
export mean
export GenericCellField
export CellQuadrature
export Integrand
export ∫
export CellDof
export get_normal_vector
export get_cell_measure
export Interpolable

export make_inverse_table
export compute_cell_points_from_vector_of_points

export DomainContribution
export num_domains
export get_domains
export get_contribution
export add_contribution!
export Measure

export attach_dirichlet
export attach_constraints_rows
export attach_constraints_cols
export identity_constraints
export get_physical_coordinate

export CellState
export update_state!

export DiracDelta

include("CellDataInterface.jl")

include("CellFields.jl")

include("CellQuadratures.jl")

include("CellStates.jl")

include("DomainContributions.jl")

include("DiracDeltas.jl")

include("CellDofs.jl")

include("AttachDirichlet.jl")

include("AttachConstraints.jl")

end # module