text:
import ..CombinedParsers: Repeat_max
## TODO:
# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC2
# affect . ^ $
# alt("(*CR)" => with_doc(CharNotIn('\r'), "carriage return"),
#     "(*LF)" => with_doc(CharNotIn('\r'), "linefeed"),
#     "(*CRLF)" => with_doc("carriage return, followed by linefeed"),
#     "(*ANYCRLF)" => with_doc("any of the three above"),
#     "(*ANY)" => with_doc("all Unicode newline sequences"))

#(*NO_AUTO_POSSESS)
#(*LIMIT_MATCH=d)
#(*LIMIT_RECURSION=d)
whitespace_string = " \t\U0085\U200E\U200F\U2028\U2029"*"\U2029\U000C\U000B"
whitespace_char = CharIn("[:space:]",whitespace_string)
meta_chars = raw"\^$.[|()?*+{"

# The horizontal space characters are:
horizontal_space=(
    '\U0009', # "Horizontal tab (HT)"),
    '\U0020', # "Space"),
    '\U00A0', # "Non-break space"),
    '\U1680', # "Ogham space mark"),
    '\U180E', # "Mongolian vowel separator"),
    '\U2000', # "En quad"),
    '\U2001', # "Em quad"),
    '\U2002', # "En space"),
    '\U2003', # "Em space"),
    '\U2004', # "Three-per-em space"),
    '\U2005', # "Four-per-em space"),
    '\U2006', # "Six-per-em space"),
    '\U2007', # "Figure space"),
    '\U2008', # "Punctuation space"),
    '\U2009', # "Thin space"),
    '\U200A', # "Hair space"),
    '\U202F', # "Narrow no-break space"),
    '\U205F', # "Medium mathematical space"),
    '\U3000' # "Ideographic space"))
)

# The vertical space characters are:
vertical_space=(
    '\U000A', # "Linefeed (LF)"),
    '\U000B', # "Vertical tab (VT)"),
    '\U000C', # "Form feed (FF)"),
    '\U000D', # "Carriage return (CR)"),
    '\U0085', # "Next line (NEL)"),
    '\U2028', # "Line separator"),
    '\U2029') # "Paragraph separator"))

bsr = Atomic(Either("\r\n",
                    CharIn(raw"\n\x0b\f\r\x85", '\n','\x0b','\f','\r','\U0085', '\U2028','\U2029'))); # backslash R (BSR)

newline = bsr
inline = !Repeat(CharNotIn(vertical_space))
whitespace_maybe = !Repeat(CharIn("\\h",horizontal_space))
whitespace_horizontal = !Repeat1(CharIn("\\h",horizontal_space))
whitespace = whitespace_horizontal

export whitespace_newline
whitespace_newline = Repeat1(CharIn("\\h\\v",horizontal_space,vertical_space))

bracket_range(start) =
    with_name(:char_range,Sequence(start,
        skip_whitespace_on(Base.PCRE.EXTENDED_MORE,Repeat),
        '-',
        skip_whitespace_on(Base.PCRE.EXTENDED_MORE,Repeat),
        bracket_char) do v
            if v[1] isa WithOptions && ( v[1].flags & Base.PCRE.CASELESS > 0 )
              cs = convert(Char,v[1]):convert(Char,v[5])
              CharIn("$(v[1])-$(v[5])",unique([ ( lowercase(x) for x in cs )...,
                                                ( uppercase(x) for x in cs )... ]))
            else
                cs = convert(Char,v[1]):convert(Char,v[5])
                CharIn("$(v[1])-$(v[5])",cs)
            end
        end)

skip_whitespace_on(flags, wrap=identity) =
    with_name(:skip_ws,on_options(
        flags,
        wrap(CharIn(whitespace_string...,'\n'))=>Always()))

skip_whitespace_and_comments =
    Repeat(Either(
        skip_whitespace_on(
            Base.PCRE.EXTENDED),
        ## comment
        on_options(
            Base.PCRE.EXTENDED,
            with_name(
                :comment,
                Sequence('#',Repeat(whitespace_char),
                    Repeat_until(
                        AnyChar(),
                        Sequence(Repeat(whitespace_char),Either(bsr,AtEnd())),
                        wrap = JoinSubstring
                    )) do v
                with_log(v[3],Always())
                end
            )),
        with_name(
            :comment,
            Sequence(
                "(?#",Repeat(whitespace_char),
                Repeat_until(
                    AnyChar(),
                    Sequence(Repeat(whitespace_char),')'),
                    wrap = JoinSubstring
                )) do v
        with_log(v[3],Always())
        end))) do v
            [a for a in v if !isa(a,Always)]
        end;

make_control(c) =
    let ui=UInt32(uppercase(c))
        ui > 127 && error("no control-$c character")
        Char(xor(ui, 0x40))
    end

seq_log(f::Function,a...) =
    Sequence(f, ( with_log("$i",e) for (i,e) in enumerate(a) )...)

export UnsupportedError
struct UnsupportedError <: Exception
    message::String
end
Base.showerror(io::IO, e::UnsupportedError) = print(io,"unsupported PCRE syntax ",e.message)

_integer(maxchar=3) =
    Sequence(Optional('-'),integer_base(10,1,maxchar)) do v
        if v[1]===missing
            v[2]
        else
            -v[2]
        end
    end
integer = _integer(Repeat_max)

at_linestart = Either(AtStart(),PositiveLookbehind(bsr))
lineend   = Either(AtEnd(),bsr)
at_lineend   = Either(AtEnd(),PositiveLookahead(bsr))

# pattern alternatives
# circumflex and dollar https://www.pcre.org/original/doc/html/pcrepattern.html#SEC6
pattern = Either{CombinedParser}(
    Any[ on_options(Base.PCRE.MULTILINE, '^' => at_linestart),
         parser('^' => AtStart()),
         on_options(Base.PCRE.MULTILINE, '$' => at_lineend),
         on_options(Base.PCRE.DOLLAR_ENDONLY, '$' => AtEnd()),
         parser('$' => Either(AtEnd(),
                              PositiveLookahead(Sequence(2,'\n',AtEnd()))))
         ]
);

# Either allows adding alternatives qith push!, that themselves use the Either object for recursive parsers.

# https://www.regular-expressions.info/refcharacters.html
# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC4
char =  Either(
    CharNotIn(meta_chars),
    Sequence(2,'\\', CharIn(meta_chars))) do v
        convert(CombinedParser,v)
    end
repeatable = Either{CombinedParser}(Any[char])

# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC5
escape_sequence(stop=AtEnd()) =
    Sequence(2,"\\Q",
        Repeat_until(AnyChar(),
                  Either("\\E",PositiveLookahead(stop)),
                  wrap=JoinSubstring));
push!(pattern,
      map(parser, with_name(:escape_sequence, escape_sequence())));

word_char=CharIn("\\w",UnicodeClass("L","N"),'_')
word = JoinSubstring(Repeat1(word_char)) ## "[[:alpha:] ]+"
words = JoinSubstring(Repeat1(CharIn("\\w\\h", horizontal_space..., word_char))) ## "[[:alpha:] ]+"

non_word_char=CharNotIn("\\W",UnicodeClass("L","N"),'_')
non_word = JoinSubstring(Repeat1(non_word_char)) ## "[[:alpha:] ]+"
non_word_char_ = Either(non_word_char,AtStart(),AtEnd())

@with_names word_boundary = Either(
    Sequence(PositiveLookbehind(word),PositiveLookahead(non_word_char_)),
    Sequence(PositiveLookbehind(non_word_char_),PositiveLookahead(word))
)

@with_names simple_assertion =
    Sequence(2,
        '\\',
        Either(
            'A' => AtStart(),
            map(parser('G')) do v
            @warn "limited \\G support: ignoring pcre2 startoffset"
            AtStart()
            end,
            'z' => AtEnd(),
            'Z' => PositiveLookahead(Sequence(Optional(bsr),AtEnd())),
            'b' => word_boundary,
            'B' => NegativeLookahead(word_boundary)
        ))
push!(pattern,simple_assertion);

push!(pattern,parser( "\\R" => bsr ));

name = JoinSubstring(
    Sequence(CharIn('a':'z','A':'Z','_'),
        Repeat(CharIn('0':'9','a':'z','A':'Z','_'))))

# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC19
@with_names backreference = map(
    Either(
        Sequence(2,'\\',Either(
            _integer(3), ## todo: maybe octal char
            Sequence(2,'g',_integer(3)),
            Sequence(2,"g{",_integer(3),'}'),
            Sequence(2,"g{",name,'}'),
            Sequence(2,"k<",name,'>'),  # perl
            Sequence(2,"k'",name,'\''), #
        )),
        Sequence(2,"(?P=",name,')'))) do v
            Backreference(v) do
                ## todo: backreference, if a capture with number (in decimal) is defined,
                ## escaped_character otherwise (if name/index not found)
                #   \ddd      character with octal code ddd, or back reference
                ## todo: error on \g<ddd>
                v isa Integer || error("capture group $v not found!")
                parse(Sequence(integer_base(8,1,3),
                               Repeat(AnyChar())) do v
                      sSequence(parser(Char(v[1])),
                                v[2]...)
                      end,
                      "$v")
            end
        end;

push!(repeatable,
      on_options(
          Base.PCRE.CASELESS,
          map(p->set_options(Base.PCRE.CASELESS,p),
              backreference)
      ));

push!(repeatable,
      backreference);

@with_names escaped_character =
    Sequence(2, '\\',
             Either(
                 'a' => ('\a'), # alarm, that is, the BEL character (hex 07)
                 Sequence('c',AnyChar()) do v  # \cx "control-x", where x is any ASCII character
                 make_control(v[2])
                 end,
                 'e' => '\e',   #  escape (hex 1B)
                 'f' => '\f',   #  form feed (hex 0C)
                 'n' => '\n',   #  linefeed (hex 0A)
                 'r' => '\r',   #  carriage return (hex 0D)
                 't' => '\t',   #  tab (hex 09)
                 '"' => '"',
                 #   \0dd      character with octal code 0dd
                 Sequence('0',integer_base(8,0,2)) do v; Char(v[2]); end,
                 #   \ddd      character with octal code ddd, or back reference
                 ## Sequence(integer_base(8,3,3), transform=v->(Char(v[1]))),
                 ## see backreference, if a capture with number (in decimal) is defined
                 #   \o{ddd..} character with octal code ddd..
                 Sequence('o','{',integer_base(8),'}') do v; Char(v[3]); end,
                 #   \xhh      character with hex code hh
                 Sequence('x','{',integer_base(16),'}') do v; Char(v[3]); end,
                 #   \x{hhh..} character with hex code hhh.. (non-JavaScript mode)
                 Sequence('x',integer_base(16,0,2)) do v; Char(v[2]); end,
                 #   \uhhhh    character with hex code hhhh (JavaScript mode only)
                 Sequence('h',integer_base(16,4,4)) do v; Char(v[2]); end,
                 CharNotIn('Q','E')
             ))

chartype_letter_parser = Either(
            # "any decimal digit"),
            'd' => CharIn("\\d",'0':'9'),
            # "any character that is not a decimal digit"),
            'D' => CharNotIn("\\D",'0':'9'),
            # "any horizontal white space character"),
            'h' => CharIn("\\h",horizontal_space...),
            # "any character that is not a horizontal white space character"),
            'H' => CharNotIn("\\H",horizontal_space...),
            # "any white space character"),
            's' => CharIn("\\s",horizontal_space...,vertical_space...),
            # "any character that is not a white space character"),
            'S' => CharNotIn("\\S",horizontal_space...,vertical_space...),
            # "any vertical white space character"),
            'v' => CharIn("\\v",vertical_space...),
            # "any character that is not a vertical white space character"),
            'V' => CharNotIn("\\V",vertical_space...),
            # "any "word" character"),
            'w' => word_char,
            # "any "non-word" character"),
            'W' => non_word_char,
        )

@with_names generic_character_type =
    Sequence(2,
        '\\', Either(
            chartype_letter_parser,
            Sequence(2,"p{",
                     Either(Dict(string(k)=>CharIn("\\p{$k}",UnicodeClass(v[3]))
                                 for (k,v) in CombinedParsers.unicode_classes)) ,
                     '}')

        ));
push!(repeatable,generic_character_type);

# https://www.regular-expressions.info/posixbrackets.html#class
bracket_char = let bracket_meta_chars = raw"]\^-"
    Either(
        CharNotIn(bracket_meta_chars),
        "\\b" => '\x08',
        Sequence('\\',integer_base(8,1,3)) do v
        Char(v[2])
        end,
        escaped_character
    )
end;

@with_names bracket=Sequence(
    CombinedParser,
    '[',Optional('^')
    , Repeat(0,1,Either(
        bracket_range(']'),
        ']'=>']'))
    , Repeat(Either(
        Sequence(
            2,
            "[:",
            Either(
                "alnum" => CharIn("[:alnum:]",UnicodeClass("L","N")), # Xan
                "alpha" => CharIn("[:alpha:]",UnicodeClass("L")),
                ##"ascii" => CharIn(UnicodeClass("InBasicLatin")),
                "blank" => CharIn("[:blank:]",UnicodeClass("Zs"),'\t'),
                "cntrl" => CharIn("[:cntrl:]",UnicodeClass("Cc")),
                "digit" => CharIn("[:digit:]",UnicodeClass("Nd")),
                "graph" => CharNotIn("[:space:]",UnicodeClass("Z","C")),
                "lower" => CharIn("[:lower:]",UnicodeClass("Ll")),
                "print" => CharIn("[:print:]",UnicodeClass("C")),
                "punct" => CharIn("[:punct:]",UnicodeClass("P")),
                "space" => CharIn("[:space:]",UnicodeClass("Z"),'\t','\r','\n','\v','\f'),
                "upper" => CharIn("[:upper:]",UnicodeClass("Lu")),
                "word" => CharIn("[:word:]",UnicodeClass("L","Nl","Nd","Pc")),
                "xdigit" => hex_digit,
            ),
            ":]"),
        skip_whitespace_on(Base.PCRE.EXTENDED_MORE,Repeat) => Never(),
        "\\E" => Never(),
        with_name(:escape_sequence,map(v->CharIn(v),escape_sequence())),
        generic_character_type,
        bracket_range(bracket_char),
        map(v->convert(CombinedParser,v),bracket_char),
        '^'=>'^',
        '-'=>'-'))
    , ']') do v
        r = (filter(!(x->isa(x,Never)),v[3])...,
             filter(!(x->isa(x,Never)),v[4])...)
        if v[2]===missing
            CharIn(r...)
        else
            CharNotIn(r...)
        end
    end;
push!(repeatable,bracket);

# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC17
@with_names repetition = Either(
    '+' => 1:Repeat_max,
    '*' => 0:Repeat_max,
    '?' => 0:1,
    Sequence(
        '{',
        integer,
        Optional(Sequence(
            2,',',
            Optional(integer, default=Repeat_max))),
        '}') do v
    if v[3] isa Missing
    v[2]:v[2]
    else
    v[2]:v[3]
    end::UnitRange{Int}
    end
)

@with_names quantified=Sequence(
    CombinedParser,
    repeatable,
    skip_whitespace_and_comments, ## for test 1130, preserve in map?
    Optional(repetition, default=1:1),
    skip_whitespace_and_comments,
    Optional(map(v->convert(Char,v),CharIn('+','?'))), # possessive quantifier, strip option
) do v
    pat = sSequence(v[1],v[2]...)
    result = if v[3] == 1:1
        parser(pat)
    elseif v[3]==0:1
        Optional(pat)
    else
        Repeat(v[3],pat)
    end
    if v[5] === missing
        result
    elseif v[5]=='+'
        Atomic(result)
    elseif v[5]=='?'
        Lazy(result)
    else
        result
    end
end;
push!(pattern, quantified)
throw_unsupported(p) =
    map(v -> throw(UnsupportedError(v)), p;
        throw_empty_union=false)
throw_unsupported(p,s) =
    map(v -> throw(UnsupportedError(s)), p;
        throw_empty_union=false)

pushfirst!(pattern,throw_unsupported(parser("\\K")));

# Sequences and Alternation
@with_names sequence = Repeat(Sequence(
    2,
    skip_whitespace_and_comments,
    pattern,
    skip_whitespace_and_comments)) do v
        length(v) ==1 ? v[1] : Sequence(v...)
    end;

alternations = Sequence(
    sequence, Repeat(Sequence(2, '|',sequence))) do v
        CombinedParser[v[1],v[2]...]
    end;

#  Options apply to subpattern,
#  (a(?i)b|c)
#  matches "ab", "aB", and "c".
#
#  Note, in PCRE, "
#  changes made in one alternative do carry on into
#  subsequent branches within the same subpattern. For
#  example,
#
#  (a(?i)b|c)
#
# matches "ab", "aB", "c", and "C", even though when
# matching "C" the first branch is abandoned before the
# option setting. This is because the effects of option
# settings happen at compile time. There would be some
# very weird behaviour otherwise."
@with_names pcre_option_start = Sequence(
    2,
    "(?",
    Either(Sequence(Optional('^'),
                    Either(Sequence(pcre_options,
                                    Optional(Sequence(2, '-',pcre_options), default=UInt32(0))),
                           Sequence(Tuple{UInt32,UInt32},'-',pcre_options) do v
                           (UInt32(0),v[2])
                           end)
           ) do v
    if v[1]===missing
    # The two "extended" options are not independent; unsetting either one cancels the effects of both of them.
    affects_extended = !iszero((v[2][1] | v[2][2]) & ( Base.PCRE.EXTENDED | Base.PCRE.EXTENDED_MORE ))
    v[2][1], affects_extended ? (v[2][1] | ( Base.PCRE.EXTENDED | Base.PCRE.EXTENDED_MORE )) : v[2][2]
    else
    (v[2][1],Base.PCRE.CASELESS | Base.PCRE.MULTILINE | Base.PCRE.NO_AUTO_CAPTURE | Base.PCRE.DOTALL| Base.PCRE.EXTENDED | Base.PCRE.EXTENDED_MORE | v[2][2])
    end
    end,
    '^' => (UInt32(0),Base.PCRE.CASELESS | Base.PCRE.MULTILINE | Base.PCRE.NO_AUTO_CAPTURE | Base.PCRE.DOTALL | Base.PCRE.EXTENDED  | Base.PCRE.EXTENDED_MORE )
));

options_alternations = after(
    Sequence(pcre_option_start,Optional(')')),
    Vector{CombinedParser}) do l
        set_options(l[1]..., l[2] === missing ?  Sequence(1, alternations,')') : alternations)
    end;

option_sequences = map(
    CombinedParser,
    Sequence(
        alternations,
        Repeat(options_alternations))) do v
            r = Any[ CombinedParser[e] for e in v[1] ]
            ro = v[2]
            for i in 1:length(ro)
                length(ro[i])>0 && push!(r[end],popfirst!(ro[i]))
                for x in ro[i]
                    ## if length(ro[i])>0
                    ## @show r[end],x
                    push!(r,CombinedParser[ x ])
                end
            end
            sEither( ( sSequence(x...) for x in r)... )
        end;
alternation = option_sequences;

# Atomic groups
# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC18
@with_names atomic_group=Sequence("(",Either("?>","*atomic:"),alternation,")") do v
    Atomic(v[3])
end;
push!(repeatable,atomic_group);

@with_names captured=Sequence("(",
             Either(Sequence(2,"?<",name,'>'),
                 Sequence(2,"?P<",name,'>'),
                 Sequence(2,"?'",name,"'"),
                 ""),
             alternation,
             ")") do v
                 with_name(v[2],Capture(Symbol(v[2]),v[3]))
             end;
push!(repeatable,captured);

@with_names subpattern=Sequence(
    2,
    "(?:",alternation,")");
push!(repeatable,subpattern);

@with_names lookahead=Sequence(
    2,
    "(",
    Either(Sequence(v -> Lookahead(true,v[2]),
                    Either("?=","*positive_lookahead:","*pla:"),alternation),
           Sequence(v -> Lookahead(false,v[2]),
                    Either("?!","*negative_lookahead:","*nla:"),alternation)),
    ")");
push!(repeatable,lookahead);

@with_names lookbehind=Sequence(
    2,
    "(",
    Either(Sequence(v -> Lookbehind(true,v[2]),
            Either("?<=","*positive_lookbehind:","*plb:"),alternation),
        Sequence(v -> Lookbehind(false,v[2]),
            Either("?<!","*negative_lookbehind:","*nlb:"),alternation)),
    ")");
push!(repeatable,lookbehind);

# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC19
@with_names subroutine = Sequence(
    2,"(?",
    Either(Sequence(Either('+','-',""),
            integer) do v
        Subroutine(nothing,Symbol(v[1]),v[2])
        end,
        Sequence(Either('&',"P>"), name) do v
        Subroutine(Symbol(v[2]),Symbol(""),-1)
        end),
    ')');
push!(repeatable,subroutine);

@with_names resetting_capture_numbers = Sequence(
    "(?|",
    alternation,
    ")") do v
        DupSubpatternNumbers(v[2])
    end;
push!(repeatable,resetting_capture_numbers);

@with_names condition = Either(
    Sequence(
        2,
        '(',
        Either(
            integer,
            "DEFINE",
            throw_unsupported(
                Sequence(
                    'R',
                    Either(
                        integer,
                        Sequence(2,'&',name),
                        Always())),
                "checking for pattern recursion"),
            Sequence(2,'\'',name,'\''),
            Sequence(2,'<',name,'>'),
            name),
        ')'),
    lookbehind,
    lookahead)

@with_names conditional = map(
    Sequence("(?",condition,
        sequence,
        Optional(Sequence(2,"|",sequence), default=Always()),
        ")")) do v
            c = v[2]
            if c=="DEFINE"
                Atomic(Either(Always(),v[3])) ## ignore in match
            elseif c isa Union{Integer,AbstractString}
                Conditional(Backreference(c) do
                            c == "R" && return Subroutine()
                            c isa Integer ? Backreference(()->error("?"),nothing, c) : error("no capture group $c")
                            end,
                            v[3],v[4])
                                      elseif c isa CombinedParser
                Conditional(c,v[3],v[4])
            else
                Conditional(Subroutine(c[2]),v[3],v[4])
            end
        end
push!(repeatable, conditional);

# https://www.regular-expressions.info/refbasic.html
dot = Either(
    on_options(Base.PCRE.DOTALL,'.') => AnyChar(), ## todo: allow \n matching context
    '.' => CharNotIn('\n'), ## todo: allow \n matching context
    "\\N" => CharNotIn('\n')
)
push!(repeatable,dot);

push!(repeatable,map(parser,escaped_character));

# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC13
@with_names sequence_with_options = after(
    Sequence(1,pcre_option_start,':'),CombinedParser) do v
        Sequence(1,set_options(v..., alternation),')')
    end
push!(repeatable,sequence_with_options);

# https://www.pcre.org/original/doc/html/pcrepattern.html#SEC27
@with_names backtrack_control = Sequence(
    2,"(*",
    Either(
        throw_unsupported(
            Sequence("ACCEPT",Optional(Sequence(2,":",JoinSubstring(Repeat_stop(AnyChar(),parser(')')))))),
            "ACCEPT"),
        Sequence(Either("FAIL","F"),Optional(Sequence(2,":",JoinSubstring(Repeat_stop(AnyChar(),parser(')')))))) do v; Never(); end,
        throw_unsupported(
            Sequence("PRUNE",Optional(Sequence(2,":",JoinSubstring(Repeat_stop(AnyChar(),parser(')')))))),
            "PRUNE"),
        throw_unsupported(
            Sequence("SKIP",Optional(Sequence(2,":",JoinSubstring(Repeat_stop(AnyChar(),parser(')')))))),
            "SKIP"),
        Sequence(Optional(parser("MARK")),':',
                 JoinSubstring(Repeat_stop(AnyChar(),parser(')')))) do v;
        with_log(v[3],Always());
        end,
        throw_unsupported(
            Sequence("COMMIT",Optional(Sequence(2,":",JoinSubstring(Repeat_stop(AnyChar(),parser(')')))))),
            "COMMIT"),
        throw_unsupported(
            Sequence("THEN",Optional(Sequence(2,":",JoinSubstring(Repeat_stop(AnyChar(),parser(')')))))),
            "THEN")),
    ")");
push!(pattern,backtrack_control);

export pcre_parser, @re_str, Regcomb
pcre_parser = Sequence(AtStart(),alternation,AtEnd()) do v
    ParserWithCaptures(v[2])
end

export Regcomb
function Regcomb(x)
    try
        r=parse(pcre_parser,x)
        r === nothing && error("invalid regex")
        r
    catch e
        if e isa UnsupportedError
            println(x,": ",e)
            throw(UnsupportedError(e.message))
        else
            rethrow(e)
        end
    end
end

"""
    Regcomb(x::AbstractString[, flags=""])

Syntax for `flags` in [`@re_str`](@ref).
"""
function Regcomb(x::AbstractString,flags::AbstractString)
    o = parse_options(flags)
    Regcomb(with_options(o...,x))
end

"""
    @re_str(x,flags)

Construct a `ParserWithCaptures` from PCRE regex syntax, such as `re"^[a-z]*\$"`, without interpolation and unescaping (except for
quotation mark `"` which still has to be escaped).
Plug-in replacement for PCRE string macro @r_str.

The regex also accepts one or more flags, listed after the ending quote, to change its behaviour:

- `i` enables case-insensitive matching
- `m` treats the `^` and `\$` tokens as matching the start and end of individual lines, as
  opposed to the whole string.
- `s` allows the `.` modifier to match newlines.
- `x` enables "comment mode": whitespace is ignored except when escaped with `\\`, and `#`
  is treated as starting a comment.
- `a` disables `UCP` mode (enables ASCII mode). By default `\\B`, `\\b`, `\\D`, `\\d`, `\\S`,
  `\\s`, `\\W`, `\\w`, etc. match based on Unicode character properties. With this option,
  these sequences only match ASCII characters.
- `xx` enables "extended comment mode": whitespace in bracket character matchers are ignored.

```jldoctest
julia> re"a|c"i
|🗄... Either
├─ [Aa] CharIn
└─ [cC] CharIn
::Char

julia> re"a+c"
🗄 Sequence
├─ a+  |> Repeat
└─ c
::Tuple{Array{Char,1},Char}
```

See also [`Regcomb`](@ref), [`parse_options`](@ref).
"""
macro re_str(x)
    esc(quote
        Regcomb($x)
        end)
end

macro re_str(x,flags)
    esc(quote
        Regcomb($x,$flags)
        end)
end

macro test_pcre(pattern,seq,log=false,flags="")
    quote
        let name = string($seq)
            @testset "$name" begin
                pcre=Regex($pattern,$flags)
                pc  =Regcomb($pattern,$flags)
                pcre_m = match(pcre,$seq)
                pc_m = match(pc,$seq)
                if $log
                    @info "testing r\"$($pattern)\" on \"$($seq)\"" pc_m pcre_m
                end
                @test pcre_m == pc_m
            end
        end
    end |> esc
end


text:
@safetestset "2D Ellipsoid Volume" begin
    using Test, QHull, Random, Distributions, LinearAlgebra
    using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface
    Random.seed!(633)
    n = 2
    μ = fill(0.0, n)
    for _ in 1:5
        x = randn(n, n)
        cov_mat = x' * x
        fun = sample_ellipsoid_surface
        x = map(_ -> fun(μ, n, cov_mat), 1:10_000)
        points = reduce(vcat, transpose.(x))

        v1 = volume_ellipsoid(n, cov_mat)

        hull = chull(points)
        v2 = hull.volume

        @test v1 ≈ v2 rtol = .03
    end
end

@safetestset "3D Ellipsoid Volume" begin
    using Test, QHull, Random, Distributions, LinearAlgebra
    using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface
    Random.seed!(584)
    n = 3
    μ = fill(0.0, n)
    for _ in 1:5
        x = randn(n, n)
        cov_mat = x' * x
        fun = sample_ellipsoid_surface
        x = map(_ -> fun(μ, n, cov_mat), 1:10_000)
        points = reduce(vcat, transpose.(x))

        v1 = volume_ellipsoid(n, cov_mat)

        hull = chull(points)
        v2 = hull.volume

        @test v1 ≈ v2 rtol = .03
    end
end

@safetestset "4D Ellipsoid Volume" begin
    using Test, QHull, Random, Distributions, LinearAlgebra
    using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface
    Random.seed!(28805)
    n = 4
    μ = fill(0.0, n)

    x = randn(n, n)
    cov_mat = x' * x
    fun = sample_ellipsoid_surface
    x = map(_ -> fun(μ, n, cov_mat), 1:10_000)
    points = reduce(vcat, transpose.(x))

    v1 = volume_ellipsoid(n, cov_mat)

    hull = chull(points)
    v2 = hull.volume

    @test v1 ≈ v2 rtol = .03
end

@safetestset "3D Odd Shape Volume" begin

    using Test, Distributions, ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Random, DataFrames, StatsBase
    include("volume_functions.jl")

    Random.seed!(1495)

    c = (
        # number of shapes
        n_shapes = 2,
        # number of simulations for hit or miss
        n_sim = 10_000,
        # number of dimensions
        n_dims = 3,
        # number of partitions per dimension
        n_part = 5,
        # number of cells for each shapes
        n_cells = [10,20],
        # number of starting points for the algorithm
        n_start = 1,
    )

    ratios = map(_ -> volume_sim(c), 1:10)

    tv1 = c.n_cells[1] * (1 / c.n_part)^c.n_dims
    tv2 = c.n_cells[2] * (1 / c.n_part)^c.n_dims
    true_ratio = tv1 / tv2
    ratio = mean(ratios)
    @test ratio ≈ true_ratio rtol = .2
end

@safetestset "5D Odd Shape Volume" begin

    using Test, Distributions, ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Random, DataFrames, StatsBase
    include("volume_functions.jl")

    Random.seed!(82544)

    c = (
        # number of shapes
        n_shapes = 2,
        # number of simulations for hit or miss
        n_sim = 10_000,
        # number of dimensions
        n_dims = 5,
        # number of partitions per dimension
        n_part = 5,
        # number of cells for each shapes
        n_cells = [10,20],
        # number of starting points for the algorithm
        n_start = 1,
    )

    ratios = map(_ -> volume_sim(c), 1:10)

    tv1 = c.n_cells[1] * (1 / c.n_part)^c.n_dims
    tv2 = c.n_cells[2] * (1 / c.n_part)^c.n_dims
    true_ratio = tv1 / tv2
    ratio = mean(ratios)
    @test ratio ≈ true_ratio rtol = .5
end

@safetestset "Volume 5D Polytope" begin

    using Test, ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Random, DataFrames, Distributions
    Random.seed!(8492)

    # dimensions of the hypbercue
    n_dims = 5
    # number of partitions
    n_part = 5
    # number of starting points
    n_start = 1

    #points = [round.(rand(n_dims), digits=3) for _ in 1:n_part]

    points = [
        [0.659, 0.547, 0.842, 0.349, 0.667],
        [0.121, 0.758, 0.141, 0.4, 0.955],
        [0.232, 0.696, 0.999, 0.655, 0.257],
        [0.312, 0.057, 0.977, 0.025, 0.154],
        [0.191, 0.782, 0.051, 0.697, 0.653],
    ]

    # partition boundaries
    polytopes = Polytope.(points)
    bounds = fill((0, 1), n_dims)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    options = Options(;
        radius = .1,
        bounds,
        n_iters = 10_000,
        parallel = false,
        init_parms,
        λ = .025,
        t_rate = .3,
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        polytopes
    )

    parm_names = Symbol.("p", 1:n_dims)
    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => parm_names
    )

    groups = groupby(df, :chain_id)
    mean_accept = combine(groups, :acceptance=>mean=>:mean)

    groups = groupby(df, :pattern)

    df_volume = combine(
        groups,
        x -> estimate_volume(
            model,
            p_fun,
            x,
            bounds,
            polytopes;
            parm_names,
        )
    )

    df_volume.volume = df_volume.x1 / sum(df_volume.x1)
    std(df_volume.volume)
    # based on 10M SMC samples
    true_volume = [
        0.40331
        0.10566
        0.17438
        0.0848
        0.23185
    ]
    true_df = DataFrame(pattern = 1:n_part, volume = true_volume)
    abs_error = abs.(true_volume .- df_volume.volume)
    max_abs_error = maximum(abs_error)
    mean_abs_error = mean(abs_error)
    sort!(df_volume, :volume)
    sort!(true_df, :volume)
    r = cor(df_volume.pattern, true_df.pattern)
    @test max_abs_error < .03
    @test mean_abs_error < .015
    @test r > .8
end

ext: jl

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


text:
using AugmentedGaussianProcesses
using Distributions
using Distances
using Random
ν=3.0

testfunction(x) = (1+x/ν)^(-0.5*(ν+1))
testfunction(x) = sech(sqrt(x/2))
# testfunction(x) = (1+sqrt(3*x))*exp(-sqrt(3*x))
# testfunction(x) = (1+sqrt(5*x)+5/3*x)*exp(-sqrt(5*x))
testdist = AGP.LaplaceTransformDistribution(testfunction,0.0,AGP.BromwichInverseLaplace(19,1,10,5))
@btime rand($testdist)
@btime rand(Gamma((ν+1)/2,1/ν))
rand(testdist)
# @btime AGP.draw(AugmentedGaussianProcesses.PolyaGammaDist(),1,0)
evaluate(KLDivergence(),rand(testdist,10000),rand(Gamma((ν+1)/2,1/ν),10000))
# laptrans(testdist,n=1)[1]
p = histogram(rand(testdist,10000),bins=range(0,3,length=100),normalize=true)|>display
histogram!(rand(Gamma((ν+1)/2,1/ν),10000),alpha=0.5,bins=range(0,10,length=100),normalize=true)
# histogram!([AGP.draw(AugmentedGaussianProcesses.PolyaGammaDist(),1,0) for i in 1:10000],alpha=0.5,bins=range(0,10,length=100),normalize=true)


text:
"""
	igetsimport(path_input,name_input,ext_input,load_interval)
Import data in ggp/eterna format downloaded from IGETS database

**Input**
* path_input: input path, e.g. "f:/data/igets/ba009/Level1"
* name_input: file prefix, e.g. "IGETS-SG-MIN-ba009-"
* ext_input: file extension, e.g. "00.ggp"
* load_interval: time interval to be loaded, e.g. DateTime(2000,1,1):Dates.Month(1):DateTime(2001,2,1)
* nanval: flag used for NaNs (NaN will be returned)

**Output**
dataframe containing whole loaded time interval

**Example**
Go to [IGETS](http://isdc.gfz-potsdam.de/igets-data-base/) database and download
the whole folder with data, e.g. all years from /cantley/ba009/Level1
```
# set the folder on your local machine with the data (without year)
path_input = "f:/mikolaj/data/sites/cantley/grav/sg/ca012/igets/Level1"
# set the file prefix
name_input = "IGETS-SG-MIN-ca012-"# without "YYYYMMXX.ggp" extension
# file extension
ext_input = "00.ggp"; # 00 is the file code
# Time interval to be loaded
load_interval = DateTime(2011,1,1):Dates.Month(1):DateTime(2012,2,1);
dataloaded = igetsimport(path_input,name_input,ext_input,
				load_interval,nanval=99.000000)
```
"""
function igetsimport(path_input::String,name_input::String,ext_input::String,
					load_interval;nanval::Float64=9999.999)
	dataout = [];
	for i in load_interval
		file_load = joinpath(path_input,Dates.format(i,"yyyy"),
							name_input*Dates.format(i,"yyyymm")*ext_input)
		if isfile(file_load)
			temp = readggp(file_load,nanval=nanval);
			if isempty(dataout)
				dataout = temp;
			else
				dataout = vcat(dataout,temp);
			end
		end
	end
	return dataout
end

"""
	igetsexport(datain,path_out,name_out,ext_out,exp_interval;
				units,decimal,channels,header,file_format,blockinfo,flagval)
Export data to IGETS format

**Input**
* datain: output DataFrame
* path_out: output path, e.g. "f:/data/igets/ba009/Level1"
* name_out: file prefix, e.g. "IGETS-SG-MIN-ba009-"
* ext_out: file extension, e.g. "00.ggp"
* exp_interval: time interval to be exported, e.g. DateTime(2000,1,1):Dates.Month(1):DateTime(2001,3,1)
* See `?writeggp` for the file settings

**Example**
```
timegrav = collect(DateTime(2016,1,1):Dates.Hour(1):DateTime(2016,3,31,23));
datain = DataFrame(datetime=timegrav,
					grav=rand(Float64,length(timegrav)),
					pres=rand(Float64,length(timegrav)));
path_out = "f:/mikolaj/data/sites/djougon/grav/sg/dj060/igets/Level2";
name_out = "TEST-IGETS-SG-CORMIN-dj060-"
ext_out = "00.ggp";
# Set time with one month longer than the 'timegrav' variable (to include also last month)
exp_interval = DateTime(DateTime(2016,1,1)):Dates.Month(1):DateTime(DateTime(2016,4,1));
unitsout = Dict(:grav=>"V",:pres=>"hPa")
channelsout = [:grav,:pres];
headerout = Dict("Station"=>"Djougon",
              "Instrument"=>"dj060");
igetsexport(datain,path_out,name_out,ext_out,exp_interval,
			units=unitsout,decimal=[2,2],channels=channelsout,
			header=headerout,file_format="preterna",
			flagval="9999.999");
```
"""
function igetsexport(datain::DataFrame,path_out::String,
					name_out::String,ext_out::String,
					exp_interval;
					units=[],decimal=[2],channels=[],header=[],
					file_format="preterna",flagval="9999.999")::Nothing
	for i in 1:length(exp_interval)-1
		path_yyyy = joinpath(path_out,Dates.format(exp_interval[i],"yyyy"))
		file_out = name_out*Dates.format(exp_interval[i],"yyyymm")*ext_out;
		r = findall(x->x<exp_interval[i] || x>=exp_interval[i+1],datain[:datetime])
		dataout = deepcopy(datain);
		!isempty(r) ? deleterows!(dataout,r) : nothing
		if !isempty(dataout)
			!isdir(path_yyyy) ? mkdir(path_yyyy) : nothing
			writeggp(dataout,joinpath(path_yyyy,file_out),units=units,header=header,decimal=decimal,
					channels=channels,flagval=flagval);
		end
	end
end


text:
export element

"""

Core types and functions for generic markup elements and attributes

"""

# Composite type for markup constructs that begin with < and ends with >

# Example: <svg id="root">content</svg>
struct StartAttributesEndTag
  name::String
  attribute_list::String
  function StartAttributesEndTag(name,attribute_list)
    return Dict("start" => string("<",name, attribute_list,">"), "end" => string("</",name,">"))
  end
end

# Example: <svg id="root"/>
struct StartAttributeTag
  name::String
  attribute_list::String
  function StartAttributeTag(name,attribute_list)
    return Dict("start" => string("<",name, attribute_list,"/>"))
  end
end

# Example: <svg></svg> or <svg>content</svg>
struct StartEndTag
  name::String
  function StartEndTag(name)
    return Dict("start" => string("<",name,">"), "end" => string("</",name,">"))
  end
end

# Example: <svg/>
struct StartTag
  name::String
  function StartTag(name)
    return Dict("start" => string("<",name," />"))
  end
end

# Composite type for xml attributes
struct Attribute
  key::String
  value::String
  function Attribute(key,value)
     return string(" ",key,"=\"",value,"\"")
  end
end

# Composite type for core xml elements.
struct Xml
    name::String
    attributes::Dict
    children::Bool
    el::Dict
    function Xml(name::String,attributes::Dict=Dict(),hasChildren::Bool=false)
      if !isnothing(attributes)
         attribute_list = string()
         for a in keys(attributes)
           attribute_list = string(attribute_list, Attribute(a,attributes[a]))
         end
         el = hasChildren == true ? StartAttributesEndTag(name,attribute_list) : StartAttributeTag(name,attribute_list)
      else
         el = hasChildren == true ? StartEndTag(name) : StartTag(name)
      end
      return el
     end
  end

"""
Multiple dispatch: Wrapper functions for various use cases.

"""
# Create element without attributes and children. Default element name is: "xml"
function element(name::String="xml")
  return string(Xml(name)["start"])
end

# Create element without attributes but with children
function element(name::String,child::String)
  return string(Xml(name,Dict(),true)["start"], child, Xml(name,Dict(),true)["end"])
end

# Create element with attributes but without children
function element(name::String,attributes::Dict)
  return string(Xml(name,attributes)["start"])
end

# Create element with attributes and children
function element(name::String,attributes::Dict,child::String)
  return string(Xml(name,attributes,true)["start"], child, Xml(name,attributes,true)["end"])
end


text:
using FileIO
import Base.write, Base.show

function read_data(gg::HDF5.HDF5Group)::NWBData
	_ancestry = read(gg,"ancestry")
	#Check the last entry
	_datatype = _ancestry[end]
	_data = read(gg["data"])
	if "timestamps" in names(gg)
		_timestamps = read(gg["timestamps"])
    start_time = first(_timestamps)
    if _rate in names(gg)
      _rate = read(gg, "rate")
    else
      _ate = 1.0/(start_time[2] - start_time[1])
    end
	elseif "start_time" in names(gg) && "rate" in names(gg)
		_start_time = read(gg,"start_time")*1.0u"s"
		_rate = read(gg, "rate")*1.0u"Hz"
    if isinteger(_rate)
      #pass
    else
		_timestamps = range(_start_time,(1/_rate),size(_data,1))
    end
	else
		warn("No timestanps found. Using dummy timestamps")
		_timestamps = range(1.0u"s", 1.0u"s", size(_data,1))
	end
	_help = read(gg["help"])
	_resolution = read(gg["resolution"])
	_name = read(gg["name"])
	_unit = read(gg["unit"])
	unit = eval(parse("1.0Unitful.$(_unit)"))
	if _datatype == "ElectricalSeries"
		_electrode_idx = read(gg["electrode_idx"])
		return ElectricalSeries(map(x->x*unit,_data),_help, _resolution, _electrode_idx, _name, _timestamps,_rate,_start_time)
	elseif _datatype == "SpatialSeries"
		return SpatialSeries(map(x->x*unit,_data),_help, _name, _timestamps)
  elseif _datatype == "SpikeEventSeries"
    return SpikeEventSeries(map(x->x*unit, _data), _help, _name, _timestamps)
	end
end

function FileIO.load(file::File{DataFormat{:NWB}})::Array{NWBData,1}
	#TODO: Read the hdf5 file, iterating through each group and classifying the datasets in that group according to the identified neurodata_type
	#e.g. neurodata_type == TimeSeries indicate that the datasets in this groups describes a time series object
	file = HDF5.h5open(file.filename)
	rdata = Array(NWBData,0)
	if "acquisition" in names(file)
		g = file["acquisition"]
		if "timeseries" in names(g)
			gg = g["timeseries"]
			for ggg in gg
				_data = read_data(ggg)
				push!(rdata, _data)
			end
		end
	end
	return rdata
end

function save(f::File{format"NWB"},data::NWBData)
	if isfile(f.filename)
		_mode = "r+"
	else
		_mode ="w"
	end
	h5open(f.filename, _mode) do s
		write(s, data)
	end
end

function write(s::HDF5.DataFile, data::TimeSeries)
	_path = "/acquisition/timeseries/$(data.name)"
  if exists(s, _path)
    return nothing
  end
	write(s, "$(_path)/name", data.name)
	write(s, "$(_path)/help", data.help)
	if "resolution" in fieldnames(data)
		_resolution = data.resolution
	else
		_resolution = NaN
	end
	write(s, "$(_path)/resolution", _resolution)
	write(s, "$(_path)/unit",string(unit(first(data.data))))
	if typeof(data) <: ElectricalSeries
		write(s, "$(_path)/electrode_idx", data.electrode_idx)
	end
  write(s, "$(_path)/rate", data.rate)
	if typeof(data.timestamps) <: Range
		write(s, "$(_path)/start_time", first(data.timestamps).val)
	else
    _timestamps = [t.val for t in data.timestamps+0.0u"s"] #convert o seconds
		write(s, "$(_path)/timestamps", _timestamps)
	end
	write(s, "$(_path)/ancestry", ["TimeSeries", split(string(typeof(data).name),".")[end]])
end

function write(s::HDF5.DataFile, data::SpikeEventSeries)
  invoke(write,(HDF5.DataFile, TimeSeries),s, data) #call write for the general TimeSries type first
	_path = "/acquisition/timeseries/$(data.name)"
  if !exists(s, "$(_path)/data")
    dd = HDF5.d_create(s, "$(_path)/data", HDF5.datatype(eltype(first(data.data).val)), HDF5.dataspace(size(data.data)))
    for j in 1:size(data.data,2)
      for i in 1:size(data.data,1)
        for k in 1:size(data.data,3)
          dd[i,j,k] = data.data[i,j,k].val
        end
      end
    end
  end
end

function write(s::HDF5.DataFile, data::Union{ElectricalSeries,SpatialSeries})
  invoke(write,(HDF5.DataFile, TimeSeries), s,data) #call write for the general TimeSries type first
	_path = "/acquisition/timeseries/$(data.name)"
  if !exists(s, "$(_path)/data")
    #write the data
    dd = HDF5.d_create(s, "$(_path)/data", HDF5.datatype(eltype(first(data.data).val)), HDF5.dataspace(size(data.data)))
    for j in 1:size(data.data,2)
      for i in 1:size(data.data,1)
        dd[i,j] = data.data[i,j].val
      end
    end
  end
end

"""
Return a dictionary of all datasets for each NWB datatype
"""
function get_path!(pth::Dict, o::Union{HDF5.HDF5File, HDF5.HDF5Group})
  if exists(o, "name")
    if exists(o, "ancestry")
      _an = read(o, "ancestry")
      _datatype = _an[end]
      if !(_datatype in keys(pth))
        pth[_datatype] = Array{String}(0)
      end
      push!(pth[_datatype], read(o, "name"))
    end
  else
    #get down one level
    for g in o
      get_path!(pth, g)
    end
  end
end

function get_path(tfile::FileIO.File{format"NWB"})
  pth = Dict()
	h5open(tfile.filename, "r") do ff
    get_path!(pth, ff)
  end
  pth
end

function Base.show(io::IO, X::TimeSeries)
  print(io, "$(typeof(X).name) with data ")
  println(io, "$(join(map(string, size(X.data)), "x")) Array{$(unit(first(X.data))),$(ndims(X.data))}")
  nothing
end


text:
function _find_sparse(s,eta)
    # Relative "eta" based sparsity parameter
    s ./= sqrt(sum(s.^2))
    goodies = abs.(s) .- eta*maximum(abs.(s))
    s = goodies .* sign.(s)
    goodies = findall(goodies .> 0)

    return(s,goodies)

end

function sparsify(s,eta,method="relative")

    if method == "relative"
        (s,goodies) = _find_sparse(s,eta)
    elseif method == "absolute"
        (s,goodies) = _find_sparse_l(s,eta)
    else
        throw(DomainError("Methods accepted are `relative` or `absolute`"))
    end

    elimvars = setdiff(1:length(s),goodies)
    s[elimvars] .= 0

    return((s,goodies,elimvars))

end

function _find_sparse_l(s,lambda)

    s ./= sqrt(sum(s.^2))
    goodies = abs.(s) .- lambda/2
    s = goodies .* sign.(s)
    goodies = findall(goodies .> 0)

    return(s,goodies)

end

@doc """

    Autoscale a data matrix

    Either pass on functions of location and scale estimators,
    or strings (e.g. Statistics.median or "median") or "none"
    if no centring and/or scaling needed.

    Output is a NamedTuple containing the autoscaled data, and the location
    and scale estimates.

    """ ->
function autoscale(X::Matrix,locest,scalest)

    if typeof(locest)==String
        if locest in ["mean","median"]
            locest = getfield(Statistics,Symbol(locest))
        else
            # other location estimators to be included
            if locest != "none"
                @warn("Only supported strings for median:" * "\n" *
                    "'mean', 'median', 'none'" * "\n" *
                    "Alternatively pass a function")
            end
        end
        pre_estimated_loc = false
    else
        if typeof(locest) in [Array{Float64,1},Array{Int64,1},Array{Float32,1},Array{Int32,1},
                              Array{Float64,2},Array{Int64,2},Array{Float32,2},Array{Int32,2}]
            pre_estimated_loc = true
        else
            pre_estimated_loc = false
        end
    end

    if typeof(scalest)==String
        if scalest == "std"
            scalest = getfield(Statistics,Symbol(scalest))
        else
            # Further scale estimates to be included
            if scalest != "none"
                # Further scale estimates to be included
                @warn("Only supported strings for scale:" * "\n" *
                    "'std', 'none'" * "\n" *
                    "Alternatively pass a function")
            end
        end
        pre_estimated_sca = false
    else
        if typeof(scalest) in [Array{Float64,1},Array{Int64,1},Array{Float32,1},Array{Int32,1},
                            Array{Float64,2},Array{Int64,2},Array{Float32,2},Array{Int32,2}]
            pre_estimated_sca = true
        else
            pre_estimated_sca = false
        end
    end

    (n,p) = size(X)

    if pre_estimated_loc
        x_loc_ = locest
        if length(size(x_loc_))==1
            x_loc_ = x_loc_'
        end
    else
        if locest == "none"
            x_loc_ = zeros(1,p)
        else
            x_loc_ = mapslices(x -> locest(x),X,dims=1)
        end
    end

    if pre_estimated_sca
        x_sca_ = scalest
        if length(size(x_sca_))==1
            x_sca_ = x_sca_'
        end
    else
        if scalest == "none"
            x_sca_ = ones(1,p)
        else
            x_sca_ = mapslices(x -> scalest(x),X,dims=1)
        end
    end

    if (locest == "none" && scalest == "none")
        X_as_ = X
    else
        X_as_ = (X - ones(n,1)*x_loc_) ./ (ones(n,1)*x_sca_)
    end
    return(X_as_ = X_as_, col_loc_ = x_loc_, col_sca_ = x_sca_)

end

@doc """

    Autoscale a data frame

    Either pass on functions of location and scale estimators,
    or strings (e.g. Statistics.median or "median") or "none"
    if no centring and/or scaling needed.

    Output is a NamedTuple containing the autoscaled data, and the location
    and scale estimates.

    """ ->
function autoscale(X::DataFrame,locest,scalest)

    if typeof(locest)==String
        if locest in ["mean","median"]
            locest = getfield(Statistics,Symbol(locest))
        else
            # other location estimators can be included
            if locest != "none"
                @warn("Only supported strings for median:" * "\n" *
                    "'mean', 'median', 'none'" * "\n" *
                    "Alternatively pass a function")
            end
        end
        pre_estimated_loc = false
    else
        if typeof(locest) in [Array{Float64,1},Array{Int64,1},Array{Float32,1},Array{Int32,1},
                            Array{Float64,2},Array{Int64,2},Array{Float32,2},Array{Int32,2}]
            pre_estimated_loc = true
        else
            pre_estimated_loc = false
        end
    end

    if typeof(scalest)==String
        if scalest == "std"
            scalest = getfield(Statistics,Symbol(scalest))
        else
            if scalest != "none"
                # Further scale estimates to be included
                @warn("Only supported strings for scale:" * "\n" *
                    "'std', 'none'" * "\n" *
                    "Alternatively pass a function")
            end
        end
        pre_estimated_sca = false
    else
        if typeof(scalest) in [Array{Float64,1},Array{Int64,1},Array{Float32,1},Array{Int32,1},
                            Array{Float64,2},Array{Int64,2},Array{Float32,2},Array{Int32,2}]
            pre_estimated_sca = true
        else
            pre_estimated_sca = false
        end
    end

    (n,p) = size(X)
    x_names = names(X)

    if pre_estimated_loc
        x_loc_ = locest
        if length(size(x_loc_))>1
            x_loc_ = x_loc_'
        end
    else
        if locest == "none"
            x_loc_ = zeros(p,1)
        else
            x_loc_ = colwise(locest,X)
        end
    end

    if pre_estimated_sca
        x_sca_ = scalest
        if length(size(x_sca_))>1
            x_sca_ = x_sca_'
        end
    else
        if scalest == "none"
            x_sca_ = ones(p,1)
        else
            x_sca_ = colwise(scalest,X)
        end
    end
    if (locest == "none" && scalest == "none")
        X_as_ = X
    else
        X_as_ = (Matrix(X) - ones(n,1)*x_loc_') ./ (ones(n,1)*x_sca_')
        X_as_ = DataFrame(X_as_)
    end
    names!(X_as_,x_names)
    return(X_as_ = X_as_, col_loc_ = x_loc_, col_sca_ = x_sca_)

end

@doc """

    Autoscale a data vector

    Either pass on functions of location and scale estimators,
    or strings (e.g. Statistics.median or "median") or "none"
    if no centring and/or scaling needed.

    Output is a NamedTuple containing the autoscaled data, and the location
    and scale estimates.

    """ ->
function autoscale(X::Vector,locest,scalest)

    if typeof(locest)==String
        if locest in ["mean","median"]
            locest = getfield(Statistics,Symbol(locest))
        else
            # other location estimators to be included
            if locest != "none"
                @warn("Only supported strings for median:" * "\n" *
                    "'mean', 'median', 'none'" * "\n" *
                    "Alternatively pass a function")
            end
        end
        pre_estimated_loc = false
    else
        if typeof(locest) in [Array{Float64,1},Array{Int64,1},Array{Float32,1},Array{Int32,1},Float64,Int64,Float32,Int32]
            pre_estimated_loc = true
        else
            pre_estimated_loc = false
        end
    end

    if typeof(scalest)==String
        if scalest == "std"
            scalest = getfield(Statistics,Symbol(scalest))
        else
            if scalest != "none"
                # Further scale estimates to be included
                @warn("Only supported strings for scale:" * "\n" *
                    "'std', 'none'" * "\n" *
                    "Alternatively pass a function")
            end
        end
        pre_estimated_sca = false
    else
        if typeof(scalest) in [Array{Float64,1},Array{Int64,1},Array{Float32,1},Array{Int32,1},Float64,Int64,Float32,Int32]
            pre_estimated_sca = true
        else
            pre_estimated_sca = false
        end
    end

    if pre_estimated_loc
        x_loc_ = locest
    else
        if locest == "none"
            x_loc_ = 0
        else
            x_loc_ = locest(X)
        end
    end

    if pre_estimated_sca
        x_sca_ = scalest
    else
        if scalest == "none"
            x_sca_ = 1
        else
            x_sca_ = scalest(X)
        end
    end
    if (locest == "none" && scalest == "none")
        X_as_ = X
    else
        X_as_ = (X .- x_loc_) ./ x_sca_
    end
    return(X_as_ = X_as_, col_loc_ = x_loc_, col_sca_ = x_sca_)

end

max_stars_repo_name: SvenSerneels/spls

max_issues_repo_name: SvenSerneels/spls

max_forks_repo_name: SvenSerneels/spls


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using Test
using MXNet

# run test in the whole directory, latest modified files
# are run first, this makes waiting time shorter when writing
# or modifying unit-tests
function test_dir(dir)
  jl_files = sort(
    filter(x -> occursin(r".*\.jl$", x), readdir(dir)),
           by = fn -> stat(joinpath(dir, fn)).mtime)
  foreach(reverse(jl_files)) do file
    include("$dir/$file")
  end
end

@info "libmxnet version => $(mx.LIB_VERSION[])"

const BASEDIR = joinpath(@__DIR__, "..")

include(joinpath(@__DIR__, "common.jl"))
@testset "MXNet Test" begin
  test_dir(joinpath(@__DIR__, "unittest"))

  # run the basic MNIST mlp example
  if haskey(ENV, "CONTINUOUS_INTEGRATION")
    @testset "MNIST Test" begin
      include(joinpath(BASEDIR, "examples", "mnist", "mlp-test.jl"))
    end
  end
end

max_stars_repo_name: allisonpaily/incubator-mxnet

max_issues_repo_name: ptrendx/mxnet

max_forks_repo_name: ptrendx/mxnet


text: _is_incremental(::Type{Gurobi.Optimizer}) = true

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

export SparsePSHDegree

@doc raw"""
`struct PSH1pBasisFcn` : 1-particle basis function specification
for bases of the type ``P \otimes Y`` with `P::ScalarBasis` and `Y::SHBasis`
"""
struct PSH1pBasisFcn <: OnepBasisFcn
   n::Int
   l::Int
   m::Int
   z::AtomicNumber
end

function PSH1pBasisFcn(t::VecOrTup)
   if length(t) == 3
      return PSH1pBasisFcn(t[1], t[2], t[3], 0)
   elseif length(t) == 4
      return PSH1pBasisFcn(t...)
   end
   error("`PSH1pBasisFcn(t::VecOrTup)` : `t` must have length 3 or 4")
end

Base.show(io::IO, b::PSH1pBasisFcn) = print(io, "znlm[$(b.z.z)|$(b.n),$(b.l),$(b.m)]")

write_dict(b::PSH1pBasisFcn) =
   Dict("__id__" => "ACE_PSH1pBasisFcn",
        "nlmz" => [ b.n, b.l, b.m, Int(b.z) ] )

read_dict(::Val{:SHIPs_PSH1pBasisFcn}, D::Dict) =
   read_dict(Val{:ACE_PSH1pBasisFcn}(), D)

read_dict(::Val{:ACE_PSH1pBasisFcn}, D::Dict) =
   PSH1pBasisFcn(D["nlmz"]...)

scaling(b::PSH1pBasisFcn, p) = b.n^p + b.l^p

abstract type AbstractPSHDegree <: AbstractDegree end

get_maxn(d::AbstractPSHDegree, maxdeg, species::Symbol) =
   get_maxn(d, maxdeg, (species,))

function get_maxn(d::AbstractPSHDegree, maxdeg, species)
   maxn = 0
   for s1 in species, s2 in species
      n = 1
      z  = AtomicNumber(s1)
      z0 = AtomicNumber(s2)
      while degree(d, PSH1pBasisFcn(n, 0, 0, z), z0) < maxdeg
         n += 1
      end
      maxn = max(maxn, n)
   end
   return maxn
end

@doc raw"""
`SparsePSHDegree` : A general sparse-grid type degree definition for
``Pr \otimes Y`` type basis functions
```math
{\rm deg}(\{n_i\}, \{l_i\})
= c_{\rm sp} \sum_i (n_i + w_{\rm L} l_i)
  + c_{\rm hc} \prod_i \max(a_{\rm hc}, b_{\rm hc} + n_i + w_{\rm L} * l_i)
```

### Constructor
```julia
SparsePSHDegree(wL = 1.5, csp = 1.0, chc = 0.0, ahc = 0.0, bhc = 0.0)
```
"""
@with_kw struct SparsePSHDegree <: AbstractPSHDegree
   wL::Float64   = 1.5
   csp::Float64  = 1.0
   chc::Float64  = 0.0
   ahc::Float64  = 0.0
   bhc::Float64  = 0.0
end

degree(d::SparsePSHDegree, phi::PSH1pBasisFcn, z0=nothing) =
      phi.n + d.wL * phi.l

function degree(d::SparsePSHDegree, pphi::VecOrTup, z0=nothing)
   if length(pphi) == 0
      return 0
   else
      return (
         d.csp * sum(  d(phi) for phi in pphi ) +
         d.chc * prod( max(d.ahc, d.bhc + d(phi)) for phi in pphi )
      )
   end
end

@doc raw"""
`SparsePSHDegreeM` : A general sparse-grid type degree definition for
``P_r \otimes Y`` type basis functions, which gives more freedom to adjust
the weights across species and correlation-orders. For simplicity, this doesn't
admit hyperbolic-cross type constructions but only the classical sparse grid.
```math
{\rm deg}(\{n_i\}_{i=1}^N, \{l_i\}_{i=1}^N)
= \sum_{i=1}^N (w^{\rm n}_i n_i + w^{\rm l}_i l_i)
```
where ``w^{\rm n}_i, w^{\rm l}_i`` may now depend on ``z_i, z_0, N``.

### Constructor

```julia
SparsePSHDegreeM(wn_fun, wl_fun)
```
where `wn_fun, wlfun` are functions that must take the arguments `(N, zi, z0)`.

### A More Practical Constructor

The functional constructor is very awkward to use, so there is an alternative
constructor, which constructs the functions wn_fun, wl_fun from a few
dictionaries.
```
SparsePSHDegreeM(Dn::Dict, Dl::Dict, Dd::Dict)
```
which will construct the functions `wn_fun, wl_fun` as follows:
```
wn_fun(N, zi, z0) = Dn[(N, zi, z0)] / Dd[(N, z0)]
wl_fun(N, zi, z0) = Dl[(N, zi, z0)] / Dd[(N, z0)]
```
If the key `(N, zi, z0)` does not exist then it will instead look for
- a key `(zi, z0)`
- or a key `"default"`
If none exist, an error is thrown. Similarly if the key `(N, z0)` does not
exist, then it will look for
- a key `N`
- a key `z0`
- a key "default".
If none of these exist, then it will throw an error.
"""
struct SparsePSHDegreeM <: AbstractPSHDegree
   wNfun
   wLfun
end

degree(d::SparsePSHDegreeM, phi::PSH1pBasisFcn, z0::AtomicNumber) =
      (    d.wNfun(1, phi.z, z0) * phi.n
         + d.wLfun(1, phi.z, z0) * phi.l )

function degree(d::SparsePSHDegreeM, b::PIBasisFcn)
   pphi = b.oneps
   z0 = b.z0
   N = length(pphi)
   if N == 0; return 0; end
   return sum( (  d.wNfun(N, phi.z, z0) * phi.n
                + d.wLfun(N, phi.z, z0) * phi.l)    for phi in pphi )
end

function _finddegree(D::Dict, N, z0)
   if haskey(D, (N, z0))
      return D[(N, z0)]
   elseif haskey(D, N)
      return D[N]
   elseif haskey(D, z0)
      return D[z0]
   elseif haskey(D, "default")
      return D["default"]
   end
   error("can't find a valid degree")
end

function _findweight(D::Dict, N, zi, z0)
   if haskey(D, (N, zi, z0))
      return D[(N, zi, z0)]
   elseif haskey(D, (zi, z0))
      return D[(zi, z0)]
   elseif haskey(D, "default")
      return D["default"]
   end
   error("SparsePSHDegreeM: no valid key found for argument $(args)")
end

SparsePSHDegreeM(DN::Dict, DL::Dict, DD::Dict) = SparsePSHDegreeM(
         (N, zi, z0) -> _findweight(DN, N, zi, z0) / _finddegree(DD, N, z0),
         (N, zi, z0) -> _findweight(DL, N, zi, z0) / _finddegree(DD, N, z0) )


text:
export GPUBackend

macro defkernels(kernels...)
  kernel_names = map(string, kernels)
  field_defs = map(kernels) do ker
    :($ker :: CUDA.CuFunction)
  end
  type_body = Expr(:block, field_defs...)

  field_inits = map(kernels, kernel_names) do ker, ker_name
    :(kernels.$ker = CUDA.CuFunction(mod, $(ker_name)))
  end
  field_init_block = Expr(:block, field_inits...)

  esc(quote
    type MochaKernels
      mod :: CUDA.CuModule

      $type_body

      MochaKernels() = begin
        mod_dir = joinpath(dirname(@__FILE__), "kernels")
        mod_path = joinpath(mod_dir, "kernels.ptx")

        # check that our module is up-to-date
        if !isfile(mod_path)
          error("Mocha CUDA kernels not found, see the documents of BACKEND on how to compile the kernels")
        else
          mod_mtime = stat(mod_path).mtime
          impl_files = glob(mod_dir, r".*.impl$")
          for i = 1:length(impl_files)
            if stat(joinpath(mod_dir, impl_files[i])).mtime > mod_mtime
              error("Mocha CUDA kernels not up-to-date. Please re-compile (see documents of BACKEND)")
            end
          end
        end

        mod = CUDA.CuModule(mod_path)
        kernels = new(mod)

        $field_init_block

        return kernels
      end
    end
  end)
end

@defkernels(
  logistic_loss_forward_float,
  logistic_loss_forward_double,
  binary_cross_entropy_loss_forward_float,
  binary_cross_entropy_loss_forward_double,
  binary_cross_entropy_loss_backward_float,
  binary_cross_entropy_loss_backward_double,
  softmax_loss_backward_float,
  softmax_loss_backward_double,
  hinge_loss_forward_float,
  hinge_loss_forward_double,
  hinge_loss_backward_float,
  hinge_loss_backward_double,
  relu_forward_float,
  relu_forward_double,
  relu_backward_float,
  relu_backward_double,
  lrelu_forward_float,
  lrelu_forward_double,
  lrelu_backward_float,
  lrelu_backward_double,
  sigmoid_forward_float,
  sigmoid_forward_double,
  sigmoid_backward_float,
  sigmoid_backward_double,
  tanh_forward_float,
  tanh_forward_double,
  tanh_backward_float,
  tanh_backward_double,
  accuracy_forward_float,
  accuracy_forward_double,
  binary_accuracy_forward_float,
  binary_accuracy_forward_double,
  argmax_forward_float,
  argmax_forward_double,
  index2onehot_forward_float,
  index2onehot_forward_double,

  add_scal_float,
  add_scal_double,
  mul_scal_float,
  mul_scal_double,
  elem_add_float,
  elem_add_double,
  elem_mul_float,
  elem_mul_double,
  elem_sub_float,
  elem_sub_double,
  elem_div_float,
  elem_div_double,
  elem_div2_float,
  elem_div2_double,
  elem_pow_fi,
  elem_pow_di,
  elem_pow_ff,
  elem_pow_dd,
  elem_log_double,
  elem_log_float,
  elem_exp_double,
  elem_exp_float,

  max_channel_pooling_forward_float,
  max_channel_pooling_forward_double,
  max_channel_pooling_backward_float,
  max_channel_pooling_backward_double,

  dense_to_padded_float,
  dense_to_padded_double,
  padded_to_dense_float,
  padded_to_dense_double,

  copy_to_shifted_float,
  copy_to_shifted_double,
  copy_from_shifted_float,
  copy_from_shifted_double,

  dropout_init,
  dropout_alloc_size,
  dropout_forward_float,
  dropout_forward_double,
  dropout_backward_float,
  dropout_backward_double,

  l1_forward_float,
  l1_forward_double,
  l1_backward_float,
  l1_backward_double,

  stdnormal_init,
  stdnormal_alloc_size,
  stdnormal_forward_float,
  stdnormal_forward_double,
)

function shutdown(mocha :: MochaKernels)
  CUDA.unload(mocha.mod)
end

type GPUBackend <: AbstractGPUBackend
  param_registry :: ParameterRegistry
  initialized    :: Bool
  cu_ctx         :: CUDA.CuContext
  cublas_ctx     :: CuBLAS.Handle
  cudnn_ctx      :: CuDNN.Handle

  mocha          :: MochaKernels

  GPUBackend() = new(ParameterRegistry(), false) # everything will be initialized later
end

function init(backend::GPUBackend)
  @assert backend.initialized == false

  m_info("Initializing CuDNN backend...")
  CUDA.init()
  dev = CUDA.CuDevice(Config.cuda_dev_id)
  backend.cu_ctx = CUDA.create_context(dev)
  backend.cublas_ctx = CuBLAS.create()
  backend.cudnn_ctx = CuDNN.create()
  backend.mocha = MochaKernels()
  backend.initialized = true
  m_info("CuDNN backend initialized!")
end

function shutdown(backend::GPUBackend)
  @assert backend.initialized == true

  m_info("Shutting down CuDNN backend...")
  # NOTE: destroy should be in reverse order of init
  shutdown(backend.mocha)
  CuDNN.destroy(backend.cudnn_ctx)
  CuBLAS.destroy(backend.cublas_ctx)
  CUDA.destroy(backend.cu_ctx)
  backend.initialized = false
  m_info("CuDNN Backend shutdown finished!")
end


text:
###################################
#    Backpropagation test
# El erro debe de ir decreciendo conforme avanzan las
# evaluaciones
###################################
using Test
@testset "Backpropagation" begin
    n = 3 # número de neuronas
    η = 0.005  # queremos que reduzca sin pasarse, de ahí que sea ""pequeño"" el learning rate
    tol = 0.5  # rango de error que permitimos ya que puede existir casos en los que el η sea demasiado grande
    data_set_size = n
    cosin(x,y)=cos(x)+sin(y) # funcion ideal
    h = RandomWeightsNN(2,n, 1) # 2 dimensión entrada 1 dimensión de salida
    X_train = (rand(Float64, (data_set_size, 2)))*3
    Y_train = map(v->cosin(v...),eachrow(X_train))
    disminuye_error = 0.0
    error_0 = error_in_data_set(
            X_train,
            Y_train,
            x->forward_propagation(h,RampFunction,x)
        )
    for i in 1:n
        backpropagation!(h, X_train, Y_train, RampFunction, derivativeRampFunction, n)

        error_1 = error_in_data_set(
            X_train,
            Y_train,
            x->forward_propagation(h,RampFunction,x)
        )
       disminuye_error +=  (error_1 < tol  + error_0) ? 1 : 0  # tolerancia
       error_0 = error_1
    end
    # Debe tenerse en cuenta de que a pesar de estar la
    # tolerancia, el η introduce cierta probabilidad de aumentar el error,
    # es por ello que introducimos esta heurísticas
     @test disminuye_error >= ceil(0.9*(n-1)) # más de noveinta porciento de los casos disminuye el error

end


text:
# Verify that using F works with both methods
SF = findstochasticflexibility(m, num_pts = 100, seed = 42)
@test abs(findstochasticflexibility(m, num_pts = 100, use_flexibility_index = true, seed = 42) - SF) <= 1e-5
@test abs(findstochasticflexibility(m, num_pts = 100, use_vulnerability_model = true, use_flexibility_index = true, seed = 42) - SF) <= 1e-5

# Verify using F works with hyperbos set
dev = ones(flex_data.numRVs) * 10
setuncertaintyset(m, :Hyperbox, [[dev]; [dev]])
solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 0.5) <= 1e-2
SF = findstochasticflexibility(m, num_pts = 100, seed = 42)
@test abs(findstochasticflexibility(m, num_pts = 100, use_flexibility_index = true, seed = 42) - SF) <= 1e-5
@test abs(findstochasticflexibility(m, num_pts = 100, use_vulnerability_model = true, use_flexibility_index = true, seed = 42) - SF) <= 1e-5

# Verify usign F works with PNorm set
setuncertaintyset(m, :PNorm, Inf)
solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 5) <= 1e-2
SF = findstochasticflexibility(m, num_pts = 100, seed = 42)
@test abs(findstochasticflexibility(m, num_pts = 100, use_flexibility_index = true, seed = 42) - SF) <= 1e-5
@test abs(findstochasticflexibility(m, num_pts = 100, use_vulnerability_model = true, use_flexibility_index = true, seed = 42) - SF) <= 1e-5


text:
module mymodule
  export helloworld

  function helloworld()
    println("Hello World!")
  end

end

mymodule.helloworld()

max_stars_repo_name: startprogress/Hello-Julia

max_issues_repo_name: startprogress/Hello-Julia

max_forks_repo_name: startprogress/Hello-Julia


text:
##############################################################################
##
## table()
##
##############################################################################

function table{T}(d::AbstractDataVector{T})
    counts = Dict{Union(T, NAtype), Int}(0)
    for i = 1:length(d)
        if has(counts, d[i])
            counts[d[i]] += 1
        else
            counts[d[i]] = 1
        end
    end
    return counts
end

##############################################################################
##
## paste()
##
##############################################################################

const letters = convert(Vector{ASCIIString}, split("abcdefghijklmnopqrstuvwxyz", ""))
const LETTERS = convert(Vector{ASCIIString}, split("ABCDEFGHIJKLMNOPQRSTUVWXYZ", ""))

# Like string(s), but preserves Vector{String} and converts
# Vector{Any} to Vector{String}.
_vstring{T <: String}(s::T) = s
_vstring{T <: String}(s::AbstractVector{T}) = s
_vstring(s::AbstractVector) = String[_vstring(x) for x in s]
_vstring(s::Any) = string(s)

function paste(s...)
    s = map(vcat * _vstring, {s...})
    sa = {s...}
    N = max(length, sa)
    res = fill("", N)
    for i in 1:length(sa)
        Ni = length(sa[i])
        k = 1
        for j = 1:N
            res[j] = string(res[j], sa[i][k])
            if k == Ni   # This recycles array elements.
                k = 1
            else
                k += 1
            end
        end
    end
    res
end

function paste_columns(d::AbstractDataFrame, sep)
    res = fill("", nrow(d))
    for j in 1:ncol(d)
        for i in 1:nrow(d)
            res[i] *= string(d[i,j])
            if j != ncol(d)
                res[i] *= sep
            end
        end
    end
    res
end
paste_columns(d::AbstractDataFrame) = paste_columns(d, "_")

##############################################################################
##
## cut()
##
##############################################################################

function cut{S, T}(x::Vector{S}, breaks::Vector{T})
    if !issorted(breaks)
        sort!(breaks)
    end
    min_x, max_x = min(x), max(x)
    if breaks[1] > min_x
        unshift!(breaks, min_x)
    end
    if breaks[end] < max_x
        push!(breaks, max_x)
    end
    refs = fill(POOLED_DATA_VEC_REF_CONVERTER(0), length(x))
    for i in 1:length(x)
        if x[i] == min_x
            refs[i] = 1
        else
            refs[i] = searchsortedfirst(breaks, x[i]) - 1
        end
    end
    n = length(breaks)
    from = map(x -> sprint(showcompact, x), breaks[1:(n - 1)])
    to = map(x -> sprint(showcompact, x), breaks[2:n])
    pool = Array(ASCIIString, n - 1)
    if breaks[1] == min_x
        pool[1] = string("[", from[1], ",", to[1], "]")
    else
        pool[1] = string("(", from[1], ",", to[1], "]")
    end
    for i in 2:(n - 1)
        pool[i] = string("(", from[i], ",", to[i], "]")
    end
    PooledDataArray(refs, pool)
end
cut(x::Vector, ngroups::Int) = cut(x, quantile(x, [1 : ngroups - 1] / ngroups))


text:
"""
Training loops on PyTorch models
"""

md"""

Is it possible to do everything via callback?

Only thing user would have to do is:
- pass in a `PyObject` model
- bass in a `PyTorchBackend(; device)` callback
- (optionally) pass in `PyObject` optimiser and loss functions (these can be
converted automatically)

The callback would do:

- `Zygote.gradient(::PyTorchParams)` should be overwritten to not calculate
    gradients, but just return the loss
- `BackwardBegin` should call `state.loss.backward()`
- `BackwardEnd` should set `state.loss = state.loss.item()` so that callbacks like
    `Metrics` will just see a number and it can be GC'd on Python side.
- `Flux.Optimise.update!(::PyTorchOptimiser, ::PyTorchParams, _)` should be overloaded
    to call `optim.step()` and `optim.zero_grad()`
"""

"""
    PyTorchBackend()

Callback that allows you to train PyTorch models with FluxTraining.jl.
To use, pass this callback to `Learner` and pass a `PyObject` PyTorch `model`.
"""
struct PyTorchBackend <: FluxTraining.Callback
    device::Any
end

PyTorchBackend() = PyTorchBackend(torch.cuda.is_available() ? "cuda" : "cpu")

import FluxTraining: init!, on, stateaccess, Read, Write, resolveconflict, Unresolvable

# During `init!`, the `Learner`'s `model`, `params`, `optimizer` and `lossfn`
# are, where possible, replaced with PyTorch versions.
function FluxTraining.init!(cb::PyTorchBackend, learner)
    learner.model = convertmodel(learner.model; device = cb.device)
    learner.params = Flux.params(learner.model)
    learner.optimizer = convertoptimiser(learner.optimizer, learner.params)
    learner.lossfn = convertlossfunction(learner.lossfn)
    return nothing
end

stateaccess(::PyTorchBackend) = (;
    model = Write(),
    params = Write(),
    optimizer = Write(),
    lossfn = Write(),
    step = Write(),
)

# `PyTorchBackend` should not be used with `ToGPU`/`ToDevice`:
resolveconflict(::PyTorchBackend, ::ToDevice) = Unresolvable()
resolveconflict(::ToDevice, ::PyTorchBackend) = Unresolvable()

function FluxTraining.on(::EpochBegin, ::AbstractTrainingPhase, cb::PyTorchBackend, learner)
    Flux.trainmode!(learner.model)
    learner.params = Flux.params(learner.model)
    learner.optimizer = convertoptimiser(learner.optimizer, learner.params)
end
function FluxTraining.on(::EpochBegin, ::AbstractValidationPhase, cb::PyTorchBackend, learner)
    Flux.testmode!(learner.model)
end

# At the beginning of every step, convert data to tensors on the correct
# device:
function FluxTraining.on(::StepBegin, ::Phase, cb::PyTorchBackend, learner)
    learner.step.xs = totensor(learner.step.xs, device=cb.device)
    learner.step.ys = totensor(learner.step.ys, device=cb.device)
end

# To make PyTorch training work with the default training loop, we'll
# implement a custom method for `Zygote.gradient` that dispatches on
# `PyTorchParams`.
Zygote.gradient(f, ::PyTorchParams) = f()

# Instead, when the `BackwardBegin` event is triggered, `loss.backward()`
# is called to do the backward pass. Importantly, this is done only during
# `AbstractTrainingPhase`s.
function FluxTraining.on(
    ::BackwardBegin,
    ::AbstractTrainingPhase,
    ::PyTorchBackend,
    learner,
)
    learner.step.loss.backward()
    learner.step.loss = learner.step.loss.item()
    # no gradient object to return
    return nothing
end

# On `StepEnd` during `AbstractValidationPhase`s, `learner.state.loss` is
# replaced by its scalar value so that callbacks like `Metrics` and can
# read it and things can be garbagecollected on the Python side. For
# `AbstractTrainingPhase`s this already happens during ´BackwardBegin`.

function FluxTraining.on(
    ::StepEnd,
    ::AbstractValidationPhase,
    ::PyTorchBackend,
    learner,
)
    learner.step.loss = learner.step.loss.item()
end

# Finally, `Flux.update!` is overloaded for `PyTorchOptimiser`

function Flux.Optimise.update!(o::PyTorchOptimiser, ::PyTorchParams, gs)
    o.optim.step()
    o.optim.zero_grad()
end


text:
abstract type distribution_hyper_params end
#Suff statistics must contain N which is the number of points associated with the cluster
abstract type sufficient_statistics end
abstract type distribution_sample end
import Base.copy

struct model_hyper_params
    distribution_hyper_params::distribution_hyper_params
    α::Float32
    total_dim::Int64
end

mutable struct cluster_parameters
    hyperparams::distribution_hyper_params
    distribution::distribution_sample
    suff_statistics::sufficient_statistics
    posterior_hyperparams::distribution_hyper_params
end

mutable struct splittable_cluster_params
    cluster_params::cluster_parameters
    cluster_params_l::cluster_parameters
    cluster_params_r::cluster_parameters
    lr_weights::AbstractArray{Float32,1}
    splittable::Bool
    logsublikelihood_hist::AbstractArray{Float32,1}
end

mutable struct thin_cluster_params{T <: distribution_sample}
    cluster_dist::T
    l_dist::T
    r_dist::T
    lr_weights::AbstractArray{Float32,1}
end

mutable struct thin_suff_stats
    cluster_suff::sufficient_statistics
    l_suff::sufficient_statistics
    r_suff::sufficient_statistics
end

mutable struct local_cluster
    cluster_params::splittable_cluster_params
    total_dim::Int64
    points_count::Int64
    l_count::Int64
    r_count::Int64
end

mutable struct local_group
    model_hyperparams::model_hyper_params
    points::AbstractArray{Float32,2}
    labels::AbstractArray{Int64,1}
    labels_subcluster::AbstractArray{Int64,1}
    local_clusters::Vector{local_cluster}
    weights::Vector{Float32}
end

mutable struct pts_less_group
    model_hyperparams::model_hyper_params
    labels::AbstractArray{Int64,1}
    labels_subcluster::AbstractArray{Int64,1}
    local_clusters::Vector{local_cluster}
    weights::Vector{Float32}
end

mutable struct local_group_stats
    labels::AbstractArray{Int64,1}
    labels_subcluster::AbstractArray{Int64,1}
    local_clusters::Vector{local_cluster}
end

mutable struct dp_parallel_sampling
    model_hyperparams::model_hyper_params
    group::local_group
end

function copy_local_cluster(c::local_cluster)
    return deepcopy(c)
end

function create_pts_less_group(group::local_group)
    return pts_less_group(
        group.model_hyperparams,
        Array(group.labels),
        Array(group.labels_subcluster),
        group.local_clusters,
        group.weights,
    )
end

function create_model_from_saved_data(
    group::pts_less_group,
    points::AbstractArray{Float32,2},
    model_hyperparams::model_hyper_params,
)
    group = local_group(
        group.model_hyperparams,
        points,
        group.labels,
        group.labels_subcluster,
        group.local_clusters,
        group.weights,
    )
    return dp_parallel_sampling(model_hyperparams, group)
end


text:
###############################################################################
##### Theta functions
###############################################################################

struct Theta
    nome        ::ComplexF64
    nome_powers ::Array{ComplexF64, 2}
    powers_num  ::Int64
    ϵ           ::Float64
end

function Theta(ω1::ComplexF64,
               ω3::ComplexF64,
               ϵ ::Float64)
    q = nome(ω1, ω3)
    nome_powers = precompute_nome_powers(q, ϵ, ω1+ω3)
    powers_num = size(nome_powers, 2)

    Theta(q, nome_powers, powers_num, ϵ)
end

"""
  nome(ω1::ComplexF64, ω3::ComplexF64)

Return the nome for the giving pair of lattice generators.
"""
function nome(
        ω1::ComplexF64,
        ω3::ComplexF64
        )
    #####
    τ = ω3 / ω1
    q = exp(1im * π * τ)
end

"""
  precompute_nome_powers(; q::ComplexF64, ϵ::Float64, z0::ComplexF64, max_n_terms ::Integer = 10)

Create array of precomputed nome powers for evaluating theta functions with given tolerance.

# RETURN:
- Array{ComplexF64, 2}(2, n)
"""
function precompute_nome_powers(
        # nome
        q           ::ComplexF64,
        # tolerance
        ϵ           ::Float64,
        # point for tolerance checking
        z0          ::ComplexF64,
        # maximum number of terms
        max_n_terms ::Integer = 10
        )
    #####
    n = 0
    while true
        # theta_1(z,q) = 2 sum_{n=0}^infty (-1)^n q^{(n+0.5)^2} sin((2n+1)z)
        dt1 = 2*(2n+1)*cos((2n+1)*z0)*q^((n+0.5)^2)
        # theta_2(z,q) = 2 sum_{n=0}^infty q^{(n+0.5)^2} cos((2n+1)z)
        dt2 = 2*(2n+1)*sin((2n+1)*z0)*q^((n+0.5)^2)
        if (abs(dt1) < ϵ) && (abs(dt1) < ϵ)
            break
        end
        n += 1
        if n > max_n_terms
            error("Series don't converge")
        end
    end
    nome_powers = Array{ComplexF64, 2}(undef, 2, n+1)
    for k in 0 : n
        nome_powers[1, k+1] = q^((k+0.5)^2)
        nome_powers[2, k+1] = nome_powers[1, k+1] * (2k+1)
    end

    return nome_powers
end

@doc """
  theta(th::Theta, z::ComplexF64)

Compute θ_1(z), θ_2(z), θ_1'(z), θ_2'(z) with given precomputed array of nome powers.
"""
function theta(
        th ::Theta,
        z  ::ComplexF64
        )
    #####
    t1  :: ComplexF64 = 0
    t2  :: ComplexF64 = 0
    dt1 :: ComplexF64 = 0
    dt2 :: ComplexF64 = 0
    sign = 1.0

    N = th.powers_num-1
    for n in 0 : N
        z_n = (2n+1)*z
        sin_z = sin(z_n)
        cos_z = cos(z_n)
        # theta_1(z,q) = 2 sum_{n=0}^N (-1)^n q^{(n+0.5)^2} sin((2n+1)z)
        t1 += 2 * sign * th.nome_powers[1, n+1] * sin_z
        # theta_2(z,q) = 2 sum_{n=0}^N q^{(n+0.5)^2} cos((2n+1)z)
        t2 += 2 * th.nome_powers[1, n+1] * cos_z
        # theta_1'(z,q) = 2 sum_{n=0}^N (-1)^n (2n+1) q^{(n+0.5)^2} cos((2n+1)z))
        dt1 += 2 * sign * th.nome_powers[2, n+1] * cos_z
        # theta_2'(z,q) = -2 sum_{n=0}^N (2n+1) q^{(n+0.5)^2} sin((2n+1)z)
        dt2 += -2 * th.nome_powers[2, n+1] * sin_z

        sign = -sign
    end

    return t1, t2, dt1, dt2
end

""""
  theta(θ::Theta; th_k::Integer, d_n::Integer,z::ComplexF64)

Compute ``\\frac{d^n}{dz^n} \\theta_k(z, q)`` with tolerance `ϵ`.
"""
function theta(
        # Theta data struct
        th ::Theta;
        # theta-function number
        th_k ::Integer,
        # derivative number
        d_n  ::Integer,
        # point for computation
        z    ::ComplexF64,
        )
   #####
    q = th.nome
    ϵ = th.ϵ

   if (th_k < 1) || (th_k > 4)
       error("Invalid theta function number")
   end

   rem_d = d_n % 4
   d_sin_sign = (rem_d == 0) || (rem_d == 1) ? 1 : -1
   d_cos_sign = (rem_d == 0) || (rem_d == 3) ? 1 : -1

   res :: ComplexF64 = 0
   if (th_k == 1) || (th_k == 2)
       res = 0
       n = 0
   else
       n = 1
       if d_n == 0
           res = 1
       else
           res = 0
       end
   end

   sign = (-1)^n
   while true
       # theta_1(z,q) = 2 sum_{n=0}^infty (-1)^n q^{(n+0.5)^2} sin((2n+1)z)
       if th_k == 1
           term = 2 * sign * (2n+1)^d_n * q^((n+0.5)^2) * d_sin_sign
           if d_n % 2 == 0
               term *= sin((2n+1)*z)
           else
               term *= cos((2n+1)*z)
           end
       end
       # theta_2(z,q) = 2 sum_{n=0}^infty q^{(n+0.5)^2} cos((2n+1)z)
       if th_k == 2
           term = 2 * (2n+1)^d_n * q^((n+0.5)^2) * d_cos_sign
           if d_n % 2 == 0
               term *= cos((2n+1)*z)
           else
               term *= sin((2n+1)*z)
           end
       end
       # theta_3(z,q) = 1 + 2 sum_{n=1}^infty q^{n^2} cos(2nz)
       if th_k == 3
           term = 2 * (2n)^d_n * q^(n^2) * d_cos_sign
           if d_n % 2 == 0
               term *= cos(2n*z)
           else
               term *= sin(2n*z)
           end
       end
       # theta_4(z,q) = 1 + 2 sum_{n=1}^infty (-1)^n q^{n^2} cos(2nz)
       if th_k == 4
           term = 2 * sign * (2n)^d_n * q^(n^2) * d_cos_sign
           if d_n % 2 == 0
               term *= cos(2n*z)
           else
               term *= sin(2n*z)
           end
       end

       res += term
       n += 1
       sign = -sign

       if abs(term) < ϵ
           break
       end
   end

   return res
end


text:
using ProgressMeter
using FFTW, LinearAlgebra
using Test

include("../src/dependencies.jl")
include("../models/WhithamGreenNaghdiGPU.jl")

function run_cpu( param )

    mesh = Mesh(param) # construct mesh of collocation points, Fourier modes, etc.
    η = exp.(-(mesh.x) .^ 2)
    v = zero(η)  # Initial data
    init = Init(mesh, η, v)
    times = Times(param.dt, param.T)
	model = WhithamGreenNaghdi(param; iterate=true)
    solver = RK4(param;model= model)
    data  = Data(mapto(model, init))
    U = copy(last(data.U))
    dt = times.dt

	problem = Problem(model, init, param)

	solve!( problem )

    real(problem.data.U)

end

function run_gpu( param )

    mesh = Mesh(param) # construct mesh of collocation points, Fourier modes, etc.
    η = exp.(-(mesh.x) .^ 2)
    v = zero(η)  # Initial data
    init = Init(mesh, η, v)
    times = Times(param.dt, param.T)
	model = WhithamGreenNaghdiGPU(param)
    solver = RK4(param;model= model)
    data  = Data(mapto(model, init))
    U = copy(last(data.U))
    dt = times.dt

	problem = Problem(model, init, param)

	solve!( problem )

    real(problem.data.U)

end

param = (
        μ = 1,
        ϵ = 1,
        N = 1024,       # number of collocation points
        L = 10π,        # mesh of size 2L
        T = 0.002,        # final time
        dt = 0.0001,    # timestep
        ns = 1,         # stores data every ns time steps
      )

# trigger compilation
@time U_cpu = run_cpu(param)
@time U_gpu = run_gpu(param)

@test U_gpu ≈ U_cpu


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libarb"
version = v"0.0.0-6c3738555d00b8b8b24a1f5e0065ef787432513c"

# Collection of sources required to build libarb
sources = [
    "https://github.com/fredrik-johansson/arb.git" =>
    "6c3738555d00b8b8b24a1f5e0065ef787432513c",

]

# Bash recipe for building across all platforms
script = raw"""
if [ $target != "x86_64-w64-mingw32" ]; then
cd $WORKSPACE/srcdir
cd arb/
./configure --prefix=$prefix --disable-static --enable-shared --with-gmp=$prefix --with-mpfr=$prefix --with-flint=$prefix
make -j${nproc}
make install

else
cd $WORKSPACE/srcdir
cd arb/
./configure --prefix=$prefix --disable-static --enable-shared --with-gmp=$prefix --with-mpfr=$prefix --with-flint=$prefix
if [ ! -f $prefx/lib/libflint-13.dll ]; then cp $prefix/bin/libflint-13.dll $prefix/lib/; fi
if [ ! -f $prefx/lib/libflint.dll ]; then cp $prefix/bin/libflint.dll $prefix/lib/; fi
#cp -n $prefix/bin/libflint-13.dll $prefix/bin/libflint.dll $prefix/lib/;
cp $prefix/bin/libflint-13.dll $prefix/bin/libflint.dll $prefix/lib/
make -j${nproc}
make install
rm $WORKSPACE/destdir/bin/libflint-13.dll
rm $WORKSPACE/destdir/bin/libflint.dll
cd $prefix/lib
mv libarb.so libarb.dll
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Windows(:x86_64),
    MacOS(:x86_64),
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libarb", :libarb)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/GMP-v6.1.2-1/build_GMP.v6.1.2.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/MPFR-v4.0.2-1/build_MPFR.v4.0.2.jl",
    "https://github.com/thofma/Flint2Builder/releases/download/c58523/build_libflint.v0.0.0-c5852387025bf144f32c0593f0ecc906c81266f1.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
module FeatureSelectors

using DataFrames
using Random

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.Utils

import AutoMLPipeline.AbsTypes: fit!, transform!
export fit!, transform!
export FeatureSelector, CatFeatureSelector, NumFeatureSelector, CatNumDiscriminator

export feature_test

# generic way to extract num/cat features by specifying their columns
mutable struct FeatureSelector <: Transformer
    name::String
    model::Dict
    args::Dict

    function FeatureSelector(args::Dict = Dict())
	default_args = Dict(
			    :name => "featureselector",
			    :columns => Int[],
			    )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

function FeatureSelector(cols::Vector{Int})
    FeatureSelector(Dict(:columns => cols))
end

function fit!(ft::FeatureSelector, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
	error("empty dataframe")
    end
    ft.model = ft.args
end

function transform!(ft::FeatureSelector, features::DataFrame)
    if features == DataFrame()
	error("empty dataframe")
    end
    return features[:,ft.model[:columns]]
end

# ----------
# automatically extracts cat features based on their inferred element non-number types
mutable struct CatFeatureSelector <: Transformer
    name::String
    model::Dict
    args::Dict

    function CatFeatureSelector(args::Dict = Dict())
	default_args = Dict(
			    :name => "catfeatsel",
			    )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

function fit!(ft::CatFeatureSelector, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
        error("empty dataframe")
    end
    catcols,_ = find_catnum_columns(features)

    # create model
    ft.model = Dict(
		    :nominal_columns => catcols
		    )
end

function transform!(ft::CatFeatureSelector, features::DataFrame)
    catcols = ft.model[:nominal_columns]
    return features[:,catcols]
end

# ---------
# automatically extracts numeric features based on their inferred element types
mutable struct NumFeatureSelector <: Transformer
    name::String
    model::Dict
    args::Dict

    function NumFeatureSelector(args::Dict = Dict())
	default_args = Dict(
			    :name => "numfeatsel"
			    )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

function fit!(ft::NumFeatureSelector, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
        error("empty dataframe")
    end
    _,realcols = find_catnum_columns(features)

    # create model
    ft.model = Dict(
		    :real_columns => realcols
		    )
end

function transform!(ft::NumFeatureSelector, features::DataFrame)
    realcols = ft.model[:real_columns]
    return features[:,realcols]
end

# ---------
# convert numeric categories to string based on count of unique elements
mutable struct CatNumDiscriminator <: Transformer
    name::String
    model::Dict
    args::Dict

    function CatNumDiscriminator(args::Dict = Dict())
	default_args = Dict(
			    :name => "catnumdisc",
			    # default max categories for numeric-encoded categories
			    :maxcategories => 24,
			    )
	cargs=nested_dict_merge(default_args,args)
	cargs[:name] = cargs[:name]*"_"*randstring(3)
	new(cargs[:name],Dict(),cargs)
    end
end

function CatNumDiscriminator(maxcat::Int)
    CatNumDiscriminator(Dict(:maxcategories=>maxcat))
end

function fit!(ft::CatNumDiscriminator, features::DataFrame, labels::Vector=[])
    if features == DataFrame()
        error("empty dataframe")
    end
    catcols,realcols = find_catnum_columns(features,ft.args[:maxcategories])

    # create model
    ft.model = Dict(
		    :real_columns => realcols,
		    :nominal_columns => catcols
		    )
end

function transform!(ft::CatNumDiscriminator, features::DataFrame)
    catcols = ft.model[:nominal_columns]
    features[!,catcols] .= features[!,catcols] .|> string
    return features
end

function feature_test()
    data = getiris()
    X = data[:,1:5]
    X[!,5] = X[!,5] .|> string
    catfeat = FeatureSelector([5])
    numfeat = FeatureSelector([1,2,3,4])
    autocat = CatFeatureSelector()
    autonum = NumFeatureSelector()
    @assert (fit_transform!(catfeat,X) .== X[:,5]) |> Matrix |> sum == 150
    @assert (fit_transform!(numfeat,X) .== X[:,1:4]) |> Matrix |> sum == 600
    @assert (fit_transform!(autocat,X) .== X[:,5]) |> Matrix |> sum == 150
    @assert (fit_transform!(autonum,X) .== X[:,1:4]) |> Matrix |> sum == 600
    catnumdata = hcat(X,repeat([1,2,3,4,5],30))
    catnum = CatNumDiscriminator()
    fit_transform!(catnum,catnumdata)
    @assert eltype(catnumdata[:,6]) <: String
    catnumdata = hcat(X,repeat([1,2,3,4,5],30))
    catnum = CatNumDiscriminator(0)
    fit_transform!(catnum,catnumdata)
    @assert eltype(catnumdata[:,6]) <: Int
end

end


text:
using LinearAlgebra: length
using ACE, ACEbase, Test, ACE.Testing
using ACE: evaluate, SymmetricBasis, PIBasis, O3, State, val
using StaticArrays
using ChainRules
import ChainRulesCore: rrule, NoTangent, ZeroTangent
using Zygote
using Zygote: @thunk
using Printf, LinearAlgebra #for the fdtestMatrix

##

@info("differentiable model test")

# construct the basis
maxdeg = 6
ord = 3
Bsel = SimpleSparseBasis(ord, maxdeg)
B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)
φ = ACE.Invariant()
basis = SymmetricBasis(φ, B1p, O3(), Bsel)

# generate a random configuration
nX = 10
cfg = ACEConfig([State(rr = rand(SVector{3, Float64})) for _ in 1:nX])

#initialize the model
np = 2
c_m = rand(SVector{np,Float64}, length(basis))
model = ACE.LinearACEModel(basis, c_m, evaluator = :standard)

evaluate(model, cfg)

##

θ = randn(np * length(basis)) ./ (1:(np*length(basis))).^2
c = reinterpret(SVector{2, Float64}, θ)
ACE.set_params!(model, c)

# FS = props -> sum([ 0.77^n * (1 + props[n]^2)^(1/n) for n = 1:length(props) ] )
FS = props -> sum( (1 .+ val.(props).^2).^0.5 )
fsmodel = cfg -> FS(evaluate(model, cfg))

# @info("check the model and gradient evaluate ok")
fsmodel(cfg)
g = Zygote.gradient(fsmodel, cfg)[1]

##

@info("Check the AD Forces for an FS-like model")
Us = randn(SVector{3, Float64}, length(cfg))
F = t -> fsmodel(ACEConfig(cfg.Xs + t * Us))
dF = t -> sum( dot(u, g.rr) for (u,g) in zip(Us, Zygote.gradient(fsmodel, cfg)[1]) )
dF(0.0)

ACEbase.Testing.fdtest(F, dF, 0.0, verbose=true)

##

mat2svecs(M::AbstractArray{T}) where {T} =
      collect(reinterpret(SVector{np, T}, M))
svecs2vec(M::AbstractVector{<: SVector{N, T}}) where {N, T} =
      collect(reinterpret(T, M))

@info("Check grad w.r.t. Params of FS-like model")

fsmodelp = θ -> ( ACE.set_params!(model, mat2svecs(θ));
                  FS(evaluate(model, cfg)) )
grad_fsmodelp = θ -> (
         ACE.set_params!(model, mat2svecs(θ));
         Zygote.gradient(model -> FS(evaluate(model, cfg)), model)[1] |> svecs2vec )
grad_fsmodelp(θ)

ACEbase.Testing.fdtest(fsmodelp, grad_fsmodelp, θ)

## THIS TEST CURRENTLY THROWS A SEGFAULT

# @info("Check AD for a second partial derivative w.r.t cfg and params")

# fsmodel1 = (model, cfg) -> FS(evaluate(model, cfg))
# grad_fsmodel1 = (model, cfg) -> Zygote.gradient(x -> fsmodel1(model, x), cfg)[1]

# y = randn(SVector{3, Float64}, length(cfg))
# loss1 = model -> sum(sum(abs2, g.rr - y)
#                      for (g, y) in zip(grad_fsmodel1(model, cfg), y))

# # check that loss and gradient evaluate ok
# loss1(model)
# g = Zygote.gradient(loss1, model)[1]  # SEGFAULT IN THIS LINE ON J1.7!!!

# # wrappers to take derivatives w.r.t. the vector or parameters
# F1 = θ -> ( ACE.set_params!(model, mat2svecs(θ));
#             loss1(model) )

# dF1 = θ -> ( ACE.set_params!(model, mat2svecs(θ));
#             Zygote.gradient(loss1, model)[1] |> svecs2vec  )

# F1(θ)
# dF1(θ)
# ACEbase.Testing.fdtest(F1, dF1, θ; verbose=true)


text:
using Gadfly, DataFrames

plot(readtable(joinpath(dirname(@__FILE__), "issue120.csv")), x=:x1, y=:x2)


text:
function get_norm_rc!(v, AT_colptr, AT_rowval, AT_nzval, n, ax)
    T = eltype(v)
    v .= zero(T)
    @inbounds @simd for j=1:n
        for i=AT_colptr[j] : (AT_colptr[j+1]-1)
            k = ax == :row ? AT_rowval[i] : j
            if abs(AT_nzval[i]) > v[k]
                v[k] = abs(AT_nzval[i])
            end
        end
    end

    v .= sqrt.(v)
    @inbounds @simd for i=1:length(v)
        if v[i] == zero(T)
            v[i] = one(T)
        end
    end
end

function mul_AT_D1_D2!(AT_colptr, AT_rowval, AT_nzval, d1, d2, r, c)
    @inbounds @simd for j=1:length(c)
        for i=AT_colptr[j]: (AT_colptr[j+1]-1)
            AT_nzval[i] /= r[AT_rowval[i]] * c[j]
        end
    end
    d1 ./= c
    d2 ./= r
end

function mul_AT_D3!(AT_colptr, AT_rowval, AT_nzval, n, d3)
    @inbounds @simd for j=1:n
        for i=AT_colptr[j]: (AT_colptr[j+1]-1)
            AT_nzval[i] *= d3[AT_rowval[i]]
        end
    end
end

function mul_Q_D!(Q_colptr, Q_rowval, Q_nzval, d, c)
    @inbounds @simd for j=1:length(d)
        for i=Q_colptr[j]: (Q_colptr[j+1]-1)
            Q_nzval[i] /= c[Q_rowval[i]] * c[j]
        end
    end
    d ./= c
end

function mul_Q_D2!(Q_colptr, Q_rowval, Q_nzval, d2)
    @inbounds @simd for j=1:length(d2)
        for i=Q_colptr[j]: (Q_colptr[j+1]-1)
            Q_nzval[i] *= d2[Q_rowval[i]] * d2[j]
        end
    end
end

function scaling_Ruiz!(fd_T0 :: QM_FloatData{T}, id :: QM_IntData, ϵ :: T;
                       max_iter :: Int = 100) where {T<:Real}

    d1, d2 = ones(T, id.ncon), ones(T, id.nvar)
    r_k, c_k = zeros(T, id.nvar), zeros(T, id.ncon)
    # r (resp. c) norm of rows of AT (resp. cols)
    # scaling: D2 * AT * D1
    get_norm_rc!(r_k, fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, id.ncon, :row)
    get_norm_rc!(c_k, fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, id.ncon, :col)
    convergence = maximum(abs.(one(T) .- r_k)) <= ϵ && maximum(abs.(one(T) .- c_k)) <= ϵ
    mul_AT_D1_D2!(fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, d1, d2, r_k, c_k)
    k = 1
    while !convergence && k < max_iter
        get_norm_rc!(r_k, fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, id.ncon, :row)
        get_norm_rc!(c_k, fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, id.ncon, :col)
        convergence = maximum(abs.(one(T) .- r_k)) <= ϵ && maximum(abs.(one(T) .- c_k)) <= ϵ
        mul_AT_D1_D2!(fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, d1, d2, r_k, c_k)
        k += 1
    end

    mul_Q_D2!(fd_T0.Q.colptr, fd_T0.Q.rowval, fd_T0.Q.nzval, d2)
    fd_T0.b .*= d1
    fd_T0.c .*= d2
    fd_T0.lvar ./= d2
    fd_T0.uvar ./= d2

    # scaling Q (symmetric)
    d3 = ones(T, id.nvar)
    r_k .= zero(T) # r_k is now norm of rows of Q
    get_norm_rc!(r_k, fd_T0.Q.colptr, fd_T0.Q.rowval, fd_T0.Q.nzval, id.nvar,:row)
    convergence = maximum(abs.(one(T) .- r_k)) <= ϵ
    mul_Q_D!(fd_T0.Q.colptr, fd_T0.Q.rowval, fd_T0.Q.nzval, d3, r_k)
    k = 1
    while !convergence && k < max_iter
        get_norm_rc!(r_k, fd_T0.Q.colptr, fd_T0.Q.rowval, fd_T0.Q.nzval, id.nvar,:row)
        convergence = maximum(abs.(one(T) .- r_k)) <= ϵ
        mul_Q_D!(fd_T0.Q.colptr, fd_T0.Q.rowval, fd_T0.Q.nzval, d3, r_k)
        k += 1
    end

    mul_AT_D3!(fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, fd_T0.AT.n, d3)
    fd_T0.c .*= d3
    fd_T0.lvar ./= d3
    fd_T0.uvar ./= d3

    return fd_T0, d1, d2, d3
end

function div_D2D3_Q_D3D2!(Q_colptr, Q_rowval, Q_nzval, d2, d3, n)
    @inbounds @simd for j=1:n
        for i=Q_colptr[j]: (Q_colptr[j+1]-1)
            Q_nzval[i] /= d2[Q_rowval[i]] * d2[j] * d3[Q_rowval[i]] * d3[j]
        end
    end
end

function div_D1_A_D2D3!(AT_colptr, AT_rowval, AT_nzval, d1, d2, d3, n)
    @inbounds @simd for j=1:n
        for i= AT_colptr[j]: (AT_colptr[j+1]-1)
            AT_nzval[i] /= d1[j] * d2[AT_rowval[i]] * d3[AT_rowval[i]]
        end
    end
end

function post_scale(d1 :: Vector{T}, d2 :: Vector{T}, d3 :: Vector{T}, pt :: Point{T}, res :: Residuals{T},
                    fd_T0 :: QM_FloatData{T}, id :: QM_IntData, Qx :: Vector{T}, ATy :: Vector{T},
                    Ax :: Vector{T}, cTx :: T, pri_obj :: T, dual_obj :: T, xTQx_2 :: T) where {T<:Real}

    pt.x .*= d2 .* d3
    div_D2D3_Q_D3D2!(fd_T0.Q.colptr, fd_T0.Q.rowval, fd_T0.Q.nzval, d2, d3, id.nvar)
    Qx = mul!(Qx, Symmetric(fd_T0.Q, :U), pt.x)
    xTQx_2 =  dot(pt.x, Qx) / 2
    div_D1_A_D2D3!(fd_T0.AT.colptr, fd_T0.AT.rowval, fd_T0.AT.nzval, d1, d2, d3, id.ncon)
    pt.y .*= d1
    ATy = mul!(ATy, fd_T0.AT, pt.y)
    Ax = mul!(Ax, fd_T0.AT', pt.x)
    fd_T0.b ./= d1
    fd_T0.c ./= d2 .* d3
    cTx = dot(fd_T0.c, pt.x)
    pri_obj = xTQx_2 + cTx + fd_T0.c0
    fd_T0.lvar .*= d2 .* d3
    fd_T0.uvar .*= d2 .* d3
    pt.s_l ./= @views d2[id.ilow] .* d3[id.ilow]
    pt.s_u ./= @views d2[id.iupp] .* d3[id.iupp]
    dual_obj = dot(fd_T0.b, pt.y) - xTQx_2 + dot(pt.s_l, view(fd_T0.lvar, id.ilow)) -
                    dot(pt.s_u, view(fd_T0.uvar, id.iupp)) + fd_T0.c0
    res.rb .= Ax .- fd_T0.b
    res.rc .= ATy .-Qx .- fd_T0.c
    res.rc[id.ilow] .+= pt.s_l
    res.rc[id.iupp] .-= pt.s_u
    #         rcNorm, rbNorm = norm(rc), norm(rb)
    res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)

    return pt, pri_obj, res
end