text:
for p in ("Knet","ArgParse","ImageMagick","MAT","Images", "Compat", "QuartzImageIO")
    Pkg.installed(p) == nothing && Pkg.add(p)
end
using Knet
!isdefined(:VGG) && include(Knet.dir("examples","vgg.jl"))
!isdefined(:Tokenizer) && include("./tokenizer.jl")
"""
julia lrcn.jl image-file-or-url

This example implements the Long-term recurrent convolutional network model from

Donahue, Jeffrey, et al. "Long-term recurrent convolutional networks for visual recognition and description."
Proceedings of the IEEE conference on computer vision and pattern recognition. 2015.

* Paper url: https://arxiv.org/pdf/1411.4389.pdf
* Project page: https://github.com/ekinakyurek/lrcn
"""
module LRCN
using Knet, AutoGrad, ArgParse, Compat, MAT, Images;
using Tokenizer;

const vggurl = "http://www.vlfeat.org/matconvnet/models/imagenet-vgg-verydeep-16.mat"
const imgurl = "https://github.com/BVLC/caffe/raw/master/examples/images/cat.jpg"
const Flickr30k_captions = "data/Flickr30k/results_20130124.token"
const MsCoCo_captions = "data/MsCoCo/captions_train2014.json"
const MsCoCo_validation = "data/MsCoCo/captions_val2014.json"
const LAYER_TYPES = ["conv", "relu", "pool", "fc", "prob"]
const cnnout = 4096;
function main(args=ARGS)
    s = ArgParseSettings()
    s.description = "LRCN.jl (c) Ekin Akyürek, 2017. Long-term Recurrent Convolutional Networks for Visual Recognition and Description"
    @add_arg_table s begin
        ("image"        ; default=imgurl; help="Image file or URL.")
        ("--model"      ; default=Knet.dir("data","imagenet-vgg-verydeep-16.mat"); help="Location of the model file")
        ("--datafiles"  ; nargs='+'; help="If provided, use first file for training, second for dev, others for test.")
        ("--loadfile"   ; help="Initialize model from file")
        ("--savefile"   ; help="Save final model to file")
        ("--generate"   ; arg_type=Int; default=0; help="If non-zero generate given number of characters.")
        ("--hidden"     ; nargs='+'; arg_type=Int; default=[1000,1000]; help="Sizes of one or more LSTM layers.")
        ("--embed"      ; arg_type=Int; default=1000; help="Size of the embedding vector.")
        ("--epochs"     ; arg_type=Int; default=10; help="Number of epochs for training.")
        ("--capnumber"  ; arg_type=Int; default=1000; help="Number of captions will be generated.")
        ("--batchsize"  ; arg_type=Int; default=25; help="Number of senteces to train on in parallel.")
        ("--lr"         ; arg_type=Float64; default=0.1; help="Initial learning rate.")
        ("--gclip"      ; arg_type=Float64; default=5.0; help="Value to clip the gradient norm at.")
        ("--seed"       ; arg_type=Int; default=-1; help="Random number seed.")
        ("--atype"      ; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--train"      ; action=:store_true; help="skip loss printing for faster run")
        ("--cnn"        ; action=:store_true; help="loads cnn weights")
        ("--extfeatures";action=:store_true; help="extract features")
        ("--flickr"     ;action=:store_true; help="works on flickr30k dataset")
        ("--coco"       ;action=:store_true; help="works on mscoco 2014 dataset")
        ("--beam_width" ;arg_type=Int;default=3;help="width of beam search")
        #TODO ("--dropout"; arg_type=Float64; default=0.0; help="Dropout probability.")
    end
    println(s.description)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    o[:seed] > 0 && srand(o[:seed])
    o[:atype] = eval(parse(o[:atype]))

    if any(f->(o[f]!=nothing), (:loadfile, :savefile, :bestfile))
        Pkg.installed("JLD")==nothing && Pkg.add("JLD") # error("Please Pkg.add(\"JLD\") to load or save files.")
        eval(Expr(:using,:JLD))
    end

    isempty(o[:datafiles]) && (o[:flickr] && push!(o[:datafiles],Flickr30k_captions))
    isempty(o[:datafiles]) && (o[:coco] && push!(o[:datafiles],MsCoCo_captions, MsCoCo_validation))

    gpu()>=0 && (KnetArray(Float32,1,1) * KnetArray(Float32,1,1)) #To initialize cudablas
    global batchsize = o[:batchsize];

    if !isempty(o[:datafiles])
      println("Tokenization starts")
      #vocab = {"words1"=>1,"words2"=>2, ....}
      #captions_dicts = [((id1, ["word1","word2",...]),length1),((id2, ["word1","word2",...]),length2),...]
      global vocab = Dict{String, Int}();
      caption_dicts = Array{Array{Tuple{Tuple{Int64,Array{String,1}},Int64},1},1}()
      Tokenizer.tokenize(vocab, caption_dicts; data_files=o[:datafiles])
      println("Tokenization finished")
      global vocab_size = length(vocab);
    end

    println("Intializing LSTM weigths")
    if o[:loadfile]==nothing
      model = initweights(o[:atype], o[:hidden], vocab_size, o[:embed])
    else
      info("Loading model from $(o[:loadfile])")
      model = map(p->convert(o[:atype],p), load(o[:loadfile], "model"))
      global vocab = load(o[:loadfile], "vocab")
      global vocab_size = length(vocab);
    end
    if o[:train]
      optim = initparams(model);
    end
    println("LSTM is initialized")

    println("$vocab_size unique words")

    global convnet = nothing; global averageImage = nothing;
    if o[:cnn]
      println("Intializing cnn_weights")
      info("Reading $(o[:model])");
      if !isfile(o[:model])
       println("Should I download the VGG model (492MB)? Enter 'y' to download, anything else to quit.")
       readline() == "y\n" || return
       download(vggurl,o[:model])
      end
      vgg = matread(o[:model]);
      params = get_params_cnn(vgg)
      global convnet = get_convnet(params...);
      global averageImage = convert(Array{Float32},vgg["meta"]["normalization"]["averageImage"])
      println("Cnn is initialized")
  end

    if o[:train] || (o[:generate]>0 && !o[:cnn])
      println("Loading existing features to train")
      o[:flickr] && (global feats = load("./data/Flickr30k/featsn.jld", "features"))
      o[:coco]   && (global feats = load("./data/MsCoCo/train2014/train_featsn.jld","features"));
      o[:coco]   && (global featsvl = load("./data/MsCoCo/val2014/val_featsn.jld", "features"));
      println("Features loaded")
    end

    if o[:generate] > 0
      if o[:cnn]
          image = read_image_data(o[:image], averageImage)
          generate(model,initstate(model,1), image, vocab, o[:generate], o[:beam_width], nothing)
      else
          if o[:flickr]
            out = open("./eval/candidates_flickr","w")
            in_out = open("./eval/candidate_ids_flickr","w")
            dict = caption_dicts[3];
          end
         if o[:coco]
            out = open("./eval/candidates.txt","w")
            in_out = open("./eval/candidate_ids.txt","w")
            dict = caption_dicts[2];
        end
          unique_ids = Dict{Int,Bool}();
          limit = o[:capnumber]

          for item in shuffle(dict)
              get!(unique_ids,item[1][1],true)
              if length(unique_ids)==limit
                 break
              end
          end

          for (id,_) in unique_ids
              o[:flickr] && generate(model,initstate(model, 1), id, vocab, o[:generate], o[:beam_width],feats; out=out, in_out=in_out);
              o[:coco] && generate(model,initstate(model, 1), id, vocab, o[:generate], o[:beam_width], featsvl; out=out, in_out=in_out);
          end
          close(out);
          close(in_out);
      end
      return;
   end

   if o[:extfeatures]
     println("extracting image features for all dataset for once")
     if o[:flickr]
       extract_features(caption_dicts[1],  "./data/Flickr30k/","feats2", "");
     elseif o[:coco]
       extract_features(caption_dicts[1], "./data/MsCoCo/train2014/","train_feats", "COCO_train2014_");
       #extract_features(caption_dicts[2], "./data/MsCoCo/val2014/","val_feats","COCO_val2014_");
     end
     println("image features extracted")
     return;
   end

  if !isempty(caption_dicts)
    println("Batching starts")
    sequence = map(t->minibatch(t, vocab, o[:batchsize]), caption_dicts)
    caption_dicts = 0; gc();
    println("Batching finished")
    println("Training starts:....")
    train!(model, optim, sequence, vocab, o)
  end

  if o[:savefile] != nothing
      info("Saving last model to $(o[:savefile])")
      save(o[:savefile], "model", model, "vocab", vocab)
  end

end

function extract_features(seq, dataset,filename, header)
  #feats = Dict{Int,Array{Float32,2}}();
  feats =   load("$(dataset)$(filename).jld", "features"); gc();
  ids = Dict{Int,Int}();
  for t=1:length(seq)
      get!(ids,seq[t][1][1],length(ids)+1)
  end
  ids = collect(keys(ids))
  l = length(ids);
  println(l);
  prefix = dataset*header
  k=1;
  for id in ids
       if !haskey(feats,id)
        println(id)
        image = read_image_data(prefix*dec(id,12)*".jpg", averageImage)
        image = convnet(image);
        feats[id] = convert(Array{Float32}, image)
       end

      if k%1000 == 1
            println(length(feats));
      end
    #  if k%10000 == 1
            #save(dataset*filename*".jld", "features", feats);
            #println(length(feats), ". image saved");
      #end
     k+=1;
  end
  println("It is finished. I do last saving")
  save(dataset*filename*".jld", "features", feats);
end

function train!(model, optim, sequence, vocab, o)
    s0 = initstate(model, o[:batchsize])
    if o[:train]
         for epoch=1:o[:epochs]
               train1(model, optim, s0, sequence[1]; batch_size=o[:batchsize], lr=o[:lr], gclip=o[:gclip], pdrop=0.4)
               if o[:savefile] != nothing
                   info("Saving last model to $(o[:savefile])")
                   save(o[:savefile], "model", model, "vocab", vocab)
               end
               losses = zeros(Float32, 2);
               losses[1] = average_loss(model, sequence[1],feats; batch_size=o[:batchsize], pdrop=0.0)
               losses[2] = average_loss(model, sequence[2], featsvl; batch_size=o[:batchsize], pdrop=0.0)
               #losses = map(d->average_loss(model, d ; batch_size=o[:batchsize], cnnout=o[:cnnout], pdrop=0.0), sequence)
               println((:epoch,epoch,:loss,losses...))
               datasheet = open("coco_e750_h700750_p_0.0.out","a+");
               println(datasheet,(:epoch,epoch,:loss,losses...));
               close(datasheet);
          end
           gpu()>=0 && Knet.cudaDeviceSynchronize()
        return
    end
    losses = map(d->average_loss(model, d ; batch_size=o[:batchsize], pdrop=0.0), sequence)
    println((:epoch,epoch,:loss,losses...))
end

function initeosbos(batch_size;atype=KnetArray{Float32})
  global eos = zeros(Int,batch_size)
  eos[:] = 1
  global bos = zeros(Int,batch_size)
  bos[:] = 2
  global unk = zeros(Int,batch_size)
  unk[:] = 3
end

function minibatch(caption_dict, word_to_index, batch_size)
    println("Minibatching starts")
    #to minibatch the train data
    if length(caption_dict) > 30000
      println("Deleting unbatchable captions: $(length(caption_dict))")
      delete_unbatchable_captions!(caption_dict, batch_size)
      println("Unbatchable captions is deleted: $(length(caption_dict))")
    else
      batch_size = 10
      println("Deleting unbatchable captions: $(length(caption_dict))")
      delete_unbatchable_captions!(caption_dict, batch_size)
      println("Unbatchable captions is deleted: $(length(caption_dict))")
    end

    initeosbos(batch_size)
    #length for one caption  caption[i][2]
    lengths = map(t->t[2], caption_dict)
    #total word count
    nbatch = div(sum(lengths), batch_size)

    #initialize output sequence, input sequence, and range of the LSTM network for each batch
    sequence = [zeros(Int,batch_size) for i=1:nbatch];
    input_ids = [zeros(Int,batch_size) for i=1:batch_size:length(lengths)]
    index = 1; l = 1; input_index = 1;

    for i=1:batch_size:length(lengths)
        l = lengths[i]
        for j=i:i+batch_size-1
          (id,words), _ = caption_dict[j]
          input_ids[input_index][j-i+1] = id
          for k=index:index+l-1
             sequence[k][j-i+1] = get(word_to_index,words[k-index+1],unk[1])
          end
        end
        index = index+l;
        input_index += 1;
    end

    println("Minibatching finished")
    return sequence,input_ids,lengths
end

function delete_unbatchable_captions!(caption_dict,batch_size)
    lengths = map(t->t[2], caption_dict)
    limit = length(lengths)-batch_size+1;
    max_length = maximum(lengths)
    current_length = lengths[1]; current_index = 1;
    ranges = Int64[];
    while current_index < limit
        if lengths[current_index+batch_size-1]==current_length
            current_index += batch_size;
        else
            old_index = current_index;
            current_index = 0;
            while current_index==0
                current_length+=1;
                if current_length > max_length
                  break
                end
                current_index = findfirst(lengths,current_length)
            end
            append!(ranges,collect(old_index:current_index-1))
        end
        if current_index >= limit
            append!(ranges,collect(current_index:length(lengths)))
            break;
        end
    end
    deleteat!(caption_dict,ranges);
    return caption_dict;
end

function train1(param, optim, state, seq; batch_size=20, lr=2.0, gclip=0.0, pdrop=0.0, atype=KnetArray{Float32})
  Knet.gc(); gc();
  sequence = seq[1];input_ids = seq[2];lengths = seq[3]

  initeosbos(batchsize);

  start_indices = ones(Int,length(sequence));
  index = 1; count = 1;
  for t = 1:batch_size:length(lengths)
        l=lengths[t];
        start_indices[count] = index
        index = index+l; count+=1;
  end
  total_sequence_size = sum(lengths);

  # index_to_char = Array(String, length(vocab))
  # for (k,v) in vocab; index_to_char[v] = k; end

  index = 1; l=1; input_index = 1; trained = 0;
  input = KnetArray(Float32, batch_size, cnnout);

 for t in shuffle(1:batch_size:length(lengths))
    l = lengths[t]
    if l>28
       continue;
    end
    input_index = Int((t-1)/batch_size + 1);
    index = start_indices[input_index];

    #   print(input_ids[input_index][1],":")
    #   print(index_to_char[bos[1]],' ')
    #   for i=index:index+l-1
    #      expected = sequence[i][1]
    #      print(index_to_char[expected],' ');
    #   end
    #  print(index_to_char[eos[1]],' ')
    #  println();

    #println("length of the senteces in the batch: ", l);
    for i=1:batch_size
           id = input_ids[input_index][i]
           feature_index = get(feats,id,nothing)
           if feature_index == nothing
             println(id)
           end
           input[i,:] = convert(atype,feature_index);
    end
    #input = convert(KnetArray{Float32},input);
    gloss = lossgradient(param, copy(state), input, sequence, index:index+l-1; pdrop = pdrop)

     trained = trained + l
     if t%2000 == 1
       println((trained*batch_size)/total_sequence_size);
       println("loss in sentence: ", loss(param, copy(state), input, sequence, index:index+l-1; pdrop = pdrop))
       Knet.knetgc();gc();
     end
    #  if gclip > 0
    #     gnorm = sqrt(mapreduce(sumabs2, +, 0, gloss))
    #      if gnorm > gclip
    #        for k=1:length(gloss)
    #          gloss[k] = (gloss[k] * gclip) / gnorm
    #        end
    #      end
    #  end
     update!(param, gloss, optim);
     isa(state,Vector{Any}) || error("State should not be Boxed.")
    end
end

function initparams(model)
  prms = Array(Any,length(model))
  for k=1:length(prms)
    prms[k] = Adam()
  end
  return prms
end

function average_loss(param, seq, featsms; batch_size=20, atype=KnetArray{Float32}, pdrop=0.0)
  sequence = seq[1]
  input_ids = seq[2]
  lengths = seq[3]

  batch_size = size(sequence[1],1)
  state = initstate(param, batch_size)
  initeosbos(batch_size)

  start_indices = ones(Int,length(sequence));
  index = 1; i = 1;
  for t = 1:batch_size:length(lengths)
        l=lengths[t];
        start_indices[i] = index
        index = index+l; i+=1;
  end

  total_sequence_size = sum(lengths);

  #index_to_char = Array(String, length(vocab))
  #for (k,v) in vocab; index_to_char[v] = k; end

  index = 1; l=1; input_index = 1; calculated = 0;
  total = 0.0; count = 0

  #for t = 1:batch_size:length(lengths)
  input = KnetArray(Float32, batch_size, cnnout);
  Knet.gc(); gc();
  for t in shuffle(1:batch_size:length(lengths))

     l = lengths[t]
     if l>28
       continue;
     end

     input_index = Int((t-1)/batch_size + 1);
     index = start_indices[input_index];

     #println("length of the senteces in the batch: ", l);
     for i=1:batch_size
           id = input_ids[input_index][i]
           feature_index = get(featsms,id,nothing)
           input[i,:] = convert(atype,feature_index);
     end

    cnn_input = input*param[end-3];
    lstm_input = param[end-2][bos,:];

    new_state = copy(state);
    for t in  index:index+l-1
        ypred = lrcn(param,new_state,cnn_input,lstm_input;pdrop=pdrop)
        ynorm = logp(ypred,2) # ypred .- log(sum(exp(ypred),2))
        index = similar(sequence[t])
        @inbounds for i=1:batch_size
            index[i] = i + (sequence[t][i]-1)*batch_size;
        end
        total += sum(ynorm[index])
        count += batch_size;
        lstm_input = param[end-2][sequence[t],:];
    end

    ypred = lrcn(param,new_state,cnn_input,lstm_input; pdrop=pdrop)
    ynorm = logp(ypred,2) # ypred .- log(sum(exp(ypred),2))
    index = similar(eos)
    @inbounds for i=1:batch_size
        index[i] = i + (eos[i]-1)*batch_size
    end
    total += sum(ynorm[index])
    count += batch_size;

    calculated = calculated + l;
    if t%100 == 1
      println(100*(calculated*batch_size/total_sequence_size), " % of training completed: ");
      println("current loss: ", -total/count);
      Knet.knetgc();gc();
    end
  end
  println(-total / count)
  return -total / count
end

function initweights(atype, hidden, vocab, embed)
  init(d...)=atype(xavier(d...))
  bias(d...)=atype(zeros(d...))
  model = Array(Any, 2*length(hidden)+5)
  X = embed
  for k = 1:length(hidden)
    H = hidden[k]
    if k==2
        X = hidden[end]
    end
    model[2k-1]    = init(X+H, 4H)
    model[2k]      = bias(1, 4H)
    model[2k][1:H] = 1 # forget gate bias = 1
    X = H
  end
  model[end-4] = init(hidden[end-1],ceil(Int,hidden[end]/2))
  model[end-3] = init(cnnout,ceil(Int,hidden[end]/2))
  model[end-2] = init(vocab,embed)
  model[end-1] = init(hidden[end],vocab)
  model[end]   = bias(1,vocab)
  return model;
end

let blank = nothing; global initstate
function initstate(model, batch)
    nlayers = div(length(model)-3,2)
    state = Array(Any, 2*nlayers)
    for k = 1:nlayers
        bias = model[2k]
        hidden = div(length(bias),4)
        if typeof(blank)!=typeof(bias) || size(blank)!=(batch,hidden)
            blank = fill!(similar(bias, batch, hidden),0)
        end
        state[2k-1] = state[2k] = blank
    end
    return state
end
end

function lstm(weight,bias,hidden,cell,input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

function lrcn(w, s, x_cnn, x_lstm; pdrop=0.0)
    x = x_lstm;
    x = dropout(x,pdrop);
    s[1],s[2] = lstm(w[1],w[2],s[1],s[2],x);
    x = s[1];
    x = x * w[end-4];
    x = hcat(x,x_cnn);
    x = dropout(x,pdrop);
    s[3],s[4] = lstm(w[3],w[4],s[3],s[4],x);
    x = s[3];
    return x * w[end-1] .+ w[end]
end

function loss(param,state,input,sequence,range; pdrop=0.0)
    total = 0.0; count = 0;
    #Beginning of sentence is multiplied by embedding matrix
    lstm_input = param[end-2][bos,:]
    #CNN output for the batch of sentences is multiplied with one parameter
    input = input * param[end-3];
    #Training starts for all words
    for t in range
        ypred = lrcn(param,state,input,lstm_input;pdrop=pdrop)
        ynorm = logp(ypred,2) # ypred .- log(sum(exp(ypred),2))
        index = similar(sequence[t])
        @inbounds for i=1:batchsize
             index[i] = i + (sequence[t][i]-1)*batchsize
        end
        total += sum(ynorm[index])
        count += batchsize;
        lstm_input = param[end-2][sequence[t],:]
    end
    #We expect from model to predict end of sentence
    ypred = lrcn(param,state,input,lstm_input; pdrop=pdrop)
    ynorm = logp(ypred,2) # ypred .- log(sum(exp(ypred),2))
    index = similar(eos)
    @inbounds for i=1:batchsize
        index[i] = i + (eos[i]-1)*batchsize
    end
    total += sum(ynorm[index])
    count += batchsize
    return -total / count
end

lossgradient = grad(loss);

function generate(param, state, input, vocab, nword, beam_width, val_feats; out=STDOUT, in_out=STDOUT)
    print("Generating Starts:");
    #Create index to word array
    initeosbos(1)
    index_to_char = Array(String, length(vocab))
    for (k,v) in vocab; index_to_char[v] = k; end
    #Initialize eos, bos and unk tokens for batchsize 1.

    if typeof(input) != Int
      # Get CNN features with convnet
      input = convnet(input);
      # Normalizing features
      input = input/sum(input);
    else
      #Loading from existing features
      println(in_out,input)
      input = get(val_feats,input,nothing);
      if input == nothing
          error("misssing features!!!!!!")
          return;
      end
      input = convert(KnetArray,reshape(input,1,cnnout));
    end
    #Beginning of sentence is multiplied by embedding matrix
    lstm_input = param[end-2][bos,:];
    #CNN output for the batch of sentences is multiplied with one parameter
    input = input * param[end-3];

    # for i=1:nword
    #     ypred = lrcn(param,copy(state),input,lstm_input)
    #     ynorm = exp(logp(ypred,2));
    #     index = sample(ynorm,100)
    #     lstm_input = param[end-2][[index],:];
    #     if index == eos[1]
    #      break;
    #    end
    #    print(out,index_to_char[index], " ");
    # end
    # println(out, ".");

    word_indices = Array{Tuple{Array{Int,1},Float32},1}();
    states = Array{typeof(state),1}();
    for i=1:beam_width
        push!(word_indices,(bos,1.0))
        push!(states,copy(state))
    end
    #Select most probable sequence in the beam
    word_indices = beam_search(word_indices,states, input, param, nword,1)
    word_indices = word_indices[1][1]
    for i=2:length(word_indices)
       if word_indices[i] == eos[1]
           break
       end
       print(out,index_to_char[word_indices[i]], " ");
    end
    println(out,".");
    println("Generating Done")
end

function beam_search(x,states,input,param,nword,current)
  new_x = Array{Tuple{Array{Int,1},Float32},1}();

  for i=1:length(x)
       current_index = x[i][1][end]
       current_probability = x[i][2]
       lstm_input = param[end-2][current_index:current_index,:];
       ypred = lrcn(param,states[i],input,lstm_input)
       ynorm = exp(logp(ypred,2));
       ynorm = convert(Array{Float32},ynorm);
       ynorm = reshape(ynorm,vocab_size);
       xmaxes = sortperm(ynorm, rev=true);
       xmaxes = xmaxes[1:length(x)]
       pmaxes = ynorm[xmaxes]*current_probability
       for j=1:length(x)
         xm = [x[i][1];xmaxes[j]]
         push!(new_x,(xm,pmaxes[j]));
       end
       if current==1
         break;
       end
  end

  sorted = sortperm(new_x, by = tuple -> last(tuple), rev=true)
  xs = new_x[sorted[1:length(x)]];
  #println(xs)
  if xs[1][1][end]==eos[1] || current>nword ##eos
    return xs;
  end
  new_states = similar(states);
  for i=1:length(xs)
      new_states[i] = copy(states[ceil(Int,sorted[i]/length(x))]);
  end
  beam_search(xs, new_states, input, param, nword, current+1)
end

function sample(p)
    p = convert(Array,p)
    r = rand()
    for c = 1:length(p)
        r -= p[c]
        r < 0 && return c
    end
end

function sample(p,N)
    indices = randperm(length(p))[1:N]
    p = convert(Array,p)
    return findmax(p[indices])[2]
end

#CNN Parameters from VLFEAT VGG.NET
function get_params_cnn(CNN; last_layer="fc7")
    layers = CNN["layers"]
    weights, operations, derivatives = [], [], []

    for l in layers
        get_layer_type(x) = startswith(l["name"], x)
        operation = filter(x -> get_layer_type(x), LAYER_TYPES)[1]
        push!(operations, operation)
        push!(derivatives, haskey(l, "weights") && length(l["weights"]) != 0)

        if derivatives[end]
            w = l["weights"]
            if operation == "conv"
                w[2] = reshape(w[2], (1,1,length(w[2]),1))
            elseif operation == "fc"
                w[1] = transpose(mat(w[1]))
            end
            push!(weights, w)
        end

        last_layer != nothing && get_layer_type(last_layer) && break
    end

    map(w -> map(KnetArray, w), weights), operations, derivatives
end

# convolutional network operations
convx(x,w)=conv4(w[1], x; padding=1, mode=1) .+ w[2]
relux = relu
poolx = pool
probx(x) = x
fcx(x,w) = w[1]*mat(x) .+ w[2];
tofunc(op) = eval(parse(string(op, "x")))
forw(x,op) = tofunc(op)(x)
forw(x,op,w) = tofunc(op)(x,w)

function get_convnet(weights, operations, derivatives)
    function convnet(xs)
        i, j = 1, 1
        num_weights, num_operations = length(weights), length(operations)
        while i <= num_operations && j <= num_weights
            if derivatives[i]
                xs = forw(xs, operations[i], weights[j])
                j += 1
            else
                xs = forw(xs, operations[i])
            end
            i += 1
        end
        return transpose(xs);
    end
end

function read_image_data(img, averageImage;atype=KnetArray{Float32})
     if contains(img,"://")
         info("Downloading $img")
         img = download(img)
     end
     a0 = load(img)
     new_size = ntuple(i->div(size(a0,i)*224,minimum(size(a0))),2)
     a1 = Images.imresize(a0, new_size)
     i1 = div(size(a1,1)-224,2)
     j1 = div(size(a1,2)-224,2)
     b1 = a1[i1+1:i1+224,j1+1:j1+224]
     cdim = length(b1[1])
     if cdim != 3
      c1 = cat(3,channelview(b1),channelview(b1),channelview(b1))
      c1 = permutedims(c1,[2,1,3])
     else
        c1 = permutedims(channelview(b1), (3,2,1))
     end
     d1 = convert(Array{Float32}, c1)
     e1 = reshape(d1[:,:,1:3], (224,224,3,1))
     f1 = (255 * e1 .- averageImage)
     g1 = permutedims(f1, [2,1,3,4])
     x1 = KnetArray(g1)
end

# To be able to load/save KnetArrays:
if Pkg.installed("JLD") != nothing
    import JLD: writeas, readas
    type KnetJLD; a::Array; end
    writeas(c::KnetArray) = KnetJLD(Array(c))
    readas(d::KnetJLD) = KnetArray(d.a)
end

!isinteractive() && !isdefined(Core.Main,:load_only) && main(ARGS)
end

max_issues_repo_name: ekinakyurek/LRCN

max_forks_repo_name: ekinakyurek/LRCN


text:
## rand.jl.  Can we generate random GMMs and sample from them?
## (c) 2013--2014 David A. van Leeuwen

## This function initializes a random GMM, with random means and random covariances
## The variances are somewhat arbitrarily chosen.  This can certainly be improved.
function Base.rand(::Type{GMM}, ng::Int, d::Int; sep=2.0, kind=:full)
    μ = sep * randn(ng, d)
    if kind==:diag
        Σ = hcat([rand(Chisq(1.0), ng) for i=1:d]...)
    elseif kind == :full
        Σ = Array{eltype(FullCov{Float64})}(undef, ng)
        for i=1:ng
            T = randn(d,d)
            Σ[i] = cholinv(T' * T / d)
        end
    else
        error("Unknown kind")
    end
    w = ones(ng)/ng
    hist = History(@sprintf("Generated random %s covariance GMM with %d Gaussians of dimension %d",
                            kind, ng, d))
    GMM(w, μ, Σ, [hist], 0)
end

## local helper
function binsearch(x::T, a::Vector{T}) where {T}
    issorted(a) || error("Array needs to be sorted")
    mi = 1
    ma = length(a)
    if x < a[mi]
        return 0
    elseif x >= a[ma]
        return ma
    end
    while ma - mi > 1
        h = mi + div(ma-mi,2)
        if x > a[h]
            mi = h
        else
            ma = h
        end
    end
    return mi
end

forcesymmetric(c::Matrix) = (Symmetric(c))

## This function samples n data points from a GMM.  This is pretty slow, probably due to the array assignments.
function Base.rand(gmm::GMM, n::Int)
    x = Array{Float64}(undef, n, gmm.d)
    ## generate indices distriuted according to weights
    index::Vector{Int} = mapslices(findall, rand(Multinomial(1, gmm.w), n).!=0,
                                   dims=1)[:]
    gmmkind = kind(gmm)
    for i=1:gmm.n
        ind = findall(index.==i)
        nₓ = length(ind)
        if gmmkind == :diag
            x[ind,:] .= (vec(gmm.μ[i,:]) .+ vec(sqrt.(gmm.Σ[i,:])) .* randn(gmm.d, nₓ))' ## v0.5 arraymageddon
        elseif gmmkind == :full
            x[ind,:] .= rand(MvNormal(vec(gmm.μ[i,:]), forcesymmetric(covar(gmm.Σ[i]))), nₓ)'
        else
            error("Unknown kind")
        end
    end
    x
end

ext: jl


text:
let
    SolverLogger{Int}()
    SolverLogger{Float64}()
    logger = SolverLogger{NTuple{3,Float64}}()

    verbosity(logger)
    debug(logger)

	# logfile = "tmp_out.log"
	# errfile = "tmp_err.log"
	# redirect_to_files(logfile,errfile) do
	# 	limit = -1
	# 	msg = "This message should be printed"
	# 	@log_info(limit,verbosity(logger),msg)
	# 	@log_info(limit,verbosity(logger),split(msg," ")...)
	# 	@log_info(limit,verbosity(logger))
	# end
	# @test countlines(logfile) == 3
	# rm(logfile)
	# rm(errfile)

end
let
    SolverLogger{Int}()
    SolverLogger{Float64}()
    logger = SolverLogger{NTuple{3,Float64}}()

    verbosity(logger)
    debug(logger)

	set_time_out_status!(logger,true)
	@test time_out_status(logger)
	set_iteration_max_out_status!(logger,true)
	@test iteration_max_out_status(logger)
	reset_solver!(logger)
	@test !time_out_status(logger)
	@test !iteration_max_out_status(logger)

    iters = iterations(logger)
    increment_iteration_count!(logger)
    @test iterations(logger) == iters + 1

    lim = 100
    set_iteration_limit!(logger,lim)
    @test iteration_limit(logger) == lim

    CRCBS.set_max_iterations!(logger,lim)
    @test max_iterations(logger) == lim

    t = 10.0
    CRCBS.set_start_time!(logger,t)
    @test start_time(logger) == t

    set_runtime_limit!(logger,t)
    @test runtime_limit(logger) == t

    set_deadline!(logger,t)
    @test deadline(logger) == t

    set_lower_bound!(logger,(1.0,0.0,0.0))
	@test lower_bound(logger) == (1.0,0.0,0.0)
    set_lower_bound!(logger,2.0)
	@test lower_bound(logger) == (2.0,0.0,0.0)
    set_lower_bound!(logger,3)
	@test lower_bound(logger) == (3.0,0.0,0.0)

    set_best_cost!(logger,(2.0,0.0,0.0))
	@test best_cost(logger) == (2.0,0.0,0.0)
    set_best_cost!(logger,3.0)
	@test best_cost(logger) == (3.0,0.0,0.0)
    set_best_cost!(logger,4)
	@test best_cost(logger) == (4.0,0.0,0.0)

    reset_solver!(logger)
    hard_reset_solver!(logger)
end
let
    logger = SolverLogger{Tuple{Int,Float64,Float64}}()
    set_lower_bound!(logger,(1,0.0,0.0))
	@test lower_bound(logger) == (1,0.0,0.0)
    set_lower_bound!(logger,2.0)
	@test lower_bound(logger) == (2,0.0,0.0)
    set_lower_bound!(logger,3)
	@test lower_bound(logger) == (3,0.0,0.0)

    set_best_cost!(logger,(2,0.0,0.0))
	@test best_cost(logger) == (2,0.0,0.0)
    set_best_cost!(logger,3.0)
	@test best_cost(logger) == (3,0.0,0.0)
    set_best_cost!(logger,4)
	@test best_cost(logger) == (4,0.0,0.0)
end


text:
# Autogenerated wrapper script for OpenMPI_jll for aarch64-apple-darwin-libgfortran5
export generate_compile_time_mpi_constants, libload_time_mpi_constants, libmpi, mpiexec

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("OpenMPI")
JLLWrappers.@declare_executable_product(generate_compile_time_mpi_constants)
JLLWrappers.@declare_library_product(libload_time_mpi_constants, "@rpath/libload_time_mpi_constants.dylib")
JLLWrappers.@declare_library_product(libmpi, "@rpath/libmpi.40.dylib")
JLLWrappers.@declare_executable_product(mpiexec)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        generate_compile_time_mpi_constants,
        "bin/generate_compile_time_mpi_constants",
    )

    JLLWrappers.@init_library_product(
        libload_time_mpi_constants,
        "lib/libload_time_mpi_constants.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmpi,
        "lib/libmpi.40.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mpiexec,
        "bin/mpiexec",
    )

    JLLWrappers.@generate_init_footer()
    ENV["OPAL_PREFIX"] = artifact_dir

end  # __init__()


text:
include("BioAFMResult.jl")

function plotPredictResult(filenameList, afmDatas::Array{Any,1}, result::Array{BioAFMResult, 1})
    Plots.gr()
    n = size(result, 1)
    plts = Plots.Plot[]

    for i in 1:n
        push!(plts, heatmap(afmDatas[i], xlabel="$(filenameList[i])"))
        push!(plts, bar(result[i].Posteriors, xlabel="PDB = $(result[i].MostPlausiblePDB), Quate = $(result[i].MostPlausibleQuate)"))
    end

    plot(plts..., layout=Plots.GridLayout(n, 2), size=(1000, 400*n))
end

function plotPredictResultWithAfmize(filenameList, afmDatas, result, pdbList, quateList)
    Plots.gr()
    n = size(result, 1)
    plts = Plots.Plot[]

    for i in 1:n
        pdb = pdbList[result[i].MostPlausiblePDB]
        pdb_rotated = MDToolbox.rotate(pdb, quateList[result[i].MostPlausibleQuate, :])
        calAfm = afmize_beta(pdb_rotated, defaultConfig())

        push!(plts, heatmap(afmDatas[i], xlabel="$(filenameList[i])"))
        push!(plts, heatmap(calAfm, xlabel="calculated"))
        push!(plts, bar(result[i].Posteriors, xlabel="PDB = $(result[i].MostPlausiblePDB), Quate = $(result[i].MostPlausibleQuate)"))
    end

    plot(plts..., layout=Plots.GridLayout(n, 3), size=(1500, 400*n))
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
@forwardfun bfmi
@forwardfun geweke
@forwardfun ess
@forwardfun rhat
@forwardfun mcse


text:
# ## クロスバリデーション
# まず線形回帰の場合のクロスバリデーションの関数を定義します。
# ```julia
# function cv_linear(X::AbstractArray,y::Vector,K::Int)
#     n = length(y); m = round(Int,n/K)
#     S = 0
#     for j in 1:K
#         test = j*m-m+1:j*m
#         train = setdiff(1:n,test)
#         β̂ = multiple_regression(X[train,:],y[train])
#         ŷ = insert_ones(X[test,:])*β̂
#         S += dot(y[test]-ŷ,y[test]-ŷ)
#     end
#     return S/n
# end
# ```
# ### 例39
using Joe, Random, Plots
using Plots.PlotMeasures # hide
Plots.reset_defaults() # hide
default( # hide
    titlefont  = font("JuliaMono", default(:titlefontsize),  ), # hide
    guidefont  = font("JuliaMono",  default(:guidefontsize),  ), # hide
    tickfont   = font("JuliaMono", default(:tickfontsize),   ), # hide
    legendfont = font("JuliaMono",  default(:legendfontsize), ), # hide
    left_margin = 30px, # hide
    bottom_margin = 30px # hide
) # hide
using Joe: cv_linear
n = 100; p = 5
Random.seed!(1)
X_39 = randn(n, p)
β = randn(p+1); β[2:3] .= 0
y_39 = insert_ones(X_39)*β + randn(n)
@show cv_linear(X_39[:,[3,4,5]],y_39,10);
@show cv_linear(X_39,y_39,10);
# あれ？変数選択した方が誤差でかいぞ？試行回数を増やして全変数を用いた場合のCV値と
# 変数選択を行った場合のCV値を可視化します。
using Joe # hide
U=Float64[]; V = Float64[]
for _ in 1:100
    global U, V, X_39, β, n # hide
    local y # hide
    y = insert_ones(X_39)*β + randn(n)
    push!(U,cv_linear(X_39[:,[3,4,5]],y,10))
    push!(V,cv_linear(X_39,y,10))
end

p39 = scatter(U,V,xlabel="変数4,5,6を選んだ時の二乗誤差",
                ylabel="全変数を選んだ時の二乗誤差",
                title="変数を多く選びすぎて過学習")
plot!(p39, x->x,xlims=(0.7,1.5),legend=false)
savefig(p39,joinpath(@OUTPUT,"fig3-1.svg")) # hide
# \fig{fig3-1}
# \lineskip
# $y=x$ の直線の上側に大多数の点がプロットされているので、
# 変数選択をした方が誤差が少ないことが分かりました。
# ただ、あまり二乗誤差の違いは大きくないです。
# ### 例40
# kの値を変えたときのCVの予測誤差の変化
n = 100; p = 5;
p40 = plot(ylims=(0.3,1.5),xlabel="k",ylabel="CVの値",
            title="k-foldのkとCVの値の関係",legend=false)
Random.seed!(1)
for _ in 1:10
    global p40 # hide
    local X, β, y, U, V # hide
    X = randn(n,p)
    β = randn(p+1)
    y = insert_ones(X) * β + randn(n)
    U = Int[]
    V = Float64[]
    for k in 2:n
        if n%k==0
            push!(U,k)
            push!(V, cv_linear(X,y,k))
        end
    end
    plot!(p40,U,V)
end
savefig(p40,joinpath(@OUTPUT,"fig3-2.svg")) # hide
# \fig{fig3-2}
# ### 例41 (FInsherのあやめ)
# あやめのデータセットでK近傍法のｋ毎に誤り率を評価します。

using ScikitLearn, StatsBase, Random, Plots
using Joe:knn
@sk_import datasets: load_iris
iris = load_iris()
X_41 = iris["data"]
y_41 = iris["target"]
n = length(y_41);
# データセットをランダムに並べ替えておきます。
Random.seed!(1)
order = StatsBase.sample(1:n,n,replace=false);
X_41 = X_41[order,:];
y_41 = y_41[order];
# 10-foldクロスバリデーションなので、テストデータを15毎に交代させることになります。
errorRate = Vector{Float64}(undef,10)
for k in 1:10
    global errorRate, X_41, y_41 # hide
    S = 0
    for top in 1:15:150
        test = top:top+14
        train = setdiff(1:150,test)
        knn_ans = knn(X_41[train,:],y_41[train],X_41[test,:],k)
        S += sum(y_41[test] .!= knn_ans)
    end
    S /= n
    errorRate[k] = S
end
p41 = plot(errorRate,xlabel="K",ylabel="誤り率",
            legend=false, title="CVによる誤り率の評価")
savefig(p41,joinpath(@OUTPUT,"fig3-3.svg")) # hide
# \fig{fig3-3}
# \lineskip
# $K \ge 3$で大きく誤り率が低下することが分かりました。
# pythonでは実行に10分程度かかったとありますが、いくらなんでも遅すぎでは？(juliaは瞬殺)
# ## 線形回帰の場合の公式(CVの公式)
# クロスバリデーションでの誤差の分析は、訓練データ毎に重回帰分析を行う必要があるので、
# 処理に時間がかかります。$\hat{\beta}_{-S}$をあらわに求めることなく、CVを
# 楽して評価できるのがこの公式です。($S$はテストデータの集合、$-S$は訓練データの集合)
# \lineskip
# $$\sum_{S} || (I - H_S)^{-1}e_S||^2$$
# ただし、$H_S = X_S(X^\mathsf{T}X)^{-1}X_S^\mathsf{T}$、
# $e_S = y_S - X_S\hat{\beta}$ です。プログラムを実装して、ベンチマークをしてみます。
# cv_fast関数は次のように実装しました。
# ```julia
# function cv_fast(X::AbstractArray,y::Vector,K::Int)
#     n = length(y); m = round(Int,n/K)
#     X = insert_ones(X)
#     H = X*((X'X)\X')
#     e = (I(n) - H)*y
#     S = 0
#     @views for j in 1:K
#         test = j*m-m+1:j*m;
#         err = (I(m)-H[test,test]) \ e[test]
#         S += dot(err,err)
#     end
#     return S/n
# end
# ```
using Joe, Random, Plots
using Joe:cv_linear,cv_fast

n = 1000; p = 5
Random.seed!(1)
X = randn(n, p)
β = randn(p+1); β[2:3] .= 0
y = insert_ones(X)*β + randn(n);
# ループのネストを浅くしたいので、K-foldのkのリストを作って、ループすることにします。
K = [i for i in 2:1000 if 1000%i ==0]
t_linear = Vector{Float64}(undef,length(K));
t_fast = similar(t_linear);
# まずは関数の慣らし運転を
@show cv_linear(X,y,10)
@show cv_fast(X,y,10)
# 値が一致することは確認できました。それでは、いざ本番!
for i in eachindex(K)
    global X, y, t_fast, t_linear,K # hide
    start = time_ns();cv_fast(X,y,K[i]);stop = time_ns()
    t_fast[i] = (stop-start)/1e9
    start = time_ns();cv_linear(X,y,K[i]);stop = time_ns()
    t_linear[i] = (stop-start)/1e9
end
p42 = plot(ylims=(0.0,0.5), xlabel = "k" ,ylabel="実行時間",
           label="cv_linear",title="cv_fastとcv_linearの比較",
           legend=:topleft)
plot!(p42,K, t_linear, label="cv_linear")
plot!(p42,K, t_fast, label="cv_fast")
savefig(p42,joinpath(@OUTPUT,"fig3-4.svg")) # hide
# \fig{fig3-4}
# \lineskip
# $n\ge100$ではcv_fastの方が実行時間が短いことが分かります。
# 原著と比べると、マシンの違いがあるとはいえ、juliaは速いということもよくわかります。
# ## ブートストラップ
# 以下のような複合型とbootstrap関数を定義しました。
# 原著のコードよりは可読性が向上したかなぁ。
# ```julia
# using Parameters
# @with_kw mutable struct bt{T<:Number}
#     original::T = 0.0
#     bias::T = 0.0
#     stderr::T = 0.0
# end
# function bootstrap(df::AbstractArray{T}, F::Function, r::Int) where T <:Number
#     m,_ = size(df)
#     original = F(df,1:m)
#     u = Vector{T}(undef, r)
#     for i in 1:r
#         index = sample(1:m,m,replace=true)
#         u[i] = F(df,index)
#     end
#     bias = mean(u) - original; stderr = std(u)
#     result = bt()
#     @pack! result = original, bias, stderr
# end
# ```
# ### 例43
# データは[ここ](https://bitbucket.org/prof-joe/statistical_learning_with_python/src/master/)からダウンロードできます。
using DelimitedFiles, Statistics
using Joe:bootstrap
function func_1(data,index)
    X=data[index,1]; Y= data[index,2]
    (var(Y) - var(X))/(var(X) + var(Y) -2cov(X,Y))
end
Portfolio = readdlm(joinpath("_assets","data","Portfolio.csv"),',', skipstart=1)
bootstrap(Portfolio,func_1,1000)

# ### 例44
# データの読み込み
using Joe,DelimitedFiles
df = readdlm(joinpath("_assets","data","crime.txt"))
X = df[:,3:4]; y =df[:,1]
β = multiple_regression(X,y)

# 回帰係数のそれぞれについて、ブートストラップ法を適用します。ここでは標準偏差だけを出力します。
std_bt = ones(3)
for j in 1:3
    function func_2(data,index)
        X = data[index,3:4];y = data[index,1]
        β = multiple_regression(X,y)
        return β[j]
    end
    std_bt[j] = bootstrap(df,func_2,1000).stderr
end
std_bt
# GLM.jlで回帰係数を推定したときの標準偏差と比較してみます。
using GLM, DataFrames
data = DataFrame(df,:auto)
ols = lm(@formula(x1 ~ x3 + x4), data)
# 標準偏差がほぼ一致することが分かりました。


text:
module MissingDocs
    export f

    "exported"
    f(x) = x

    "unexported"
    g(x) = x
end

using Documenter

for sym in [:none, :exports]
    makedocs(
        root = dirname(@__FILE__),
        source = joinpath("src", string(sym)),
        build = joinpath("build", string(sym)),
        modules = MissingDocs,
        checkdocs = sym,
        sitename = "MissingDocs Checks",
    )
end


text:
abstract BasisFun{T<:Number}
abstract BasisFun1d{T<:Number} <: BasisFun{T}

" Total number of modes of a basis function "
nmodes{T<:BasisFun}(b::T) = b.nmodes

" Total number of modes of a basis function "
length{T<:BasisFun1d}(b::T) = nmodes(b)

" Total number of modes of a basis function "
size{T<:BasisFun1d}(b::T) = (nmodes(b),)

" Can the modes be decoupled into boundary/interior modes? "
isbndryint{T<:BasisFun}(b::T) = false
isbndryint{T<:BasisFun}(::Type{T}) = false

bndidx{T<:BasisFun}(b::T) = bndidx(b.lnum)
nbndry{T<:BasisFun}(b::T) = nbndry(b.lnum)

intidx{T<:BasisFun}(b::T) = intidx(b.lnum)
ninterior{T<:BasisFun}(b::T) = ninterior(b.lnum)

seq2bi!{T}(b::BasisFun, x::AbstractVector{T}, y::AbstractVector{T}) = seq2bi!(b.lnum, x, y)
seq2bi{T}(b::BasisFun, x::AbstractVector{T}) = seq2bi(b.lnum, x)
bi2seq!{T}(b::BasisFun, x::AbstractVector{T}, y::AbstractVector{T}) = bi2seq!(b.lnum, x, y)
bi2seq{T}(b::BasisFun, x::AbstractVector{T}) = bi2seq(b.lnum, x)
bi2seq!{T}(b::BasisFun, xb::AbstractVector{T}, xi::AbstractVector{T}, y::AbstractVector{T}) =
    bi2seq!(b.lnum, xb, xi, y)
bi2seq{T}(b::BasisFun, xb::AbstractVector{T}, xi::AbstractVector{T}) = bi2seq(b.lnum, xb, xi)
seq2b!{T}(b::BasisFun, x::AbstractVector{T}, y::AbstractVector{T}) = seq2b!(b.lnum, x, y)
seq2i!{T}(b::BasisFun, x::AbstractVector{T}, y::AbstractVector{T}) = seq2i!(b.lnum, x, y)
seq2b{T}(b::BasisFun, x::AbstractVector{T}) = seq2b(b.lnum, x)
seq2i{T}(b::BasisFun, x::AbstractVector{T}) = seq2i(b.lnum, x)

"""
Compute the value of mode `i` at all points ξ

The values are returned on array `y`.
"""
function basis1d!{T<:BasisFun1d}(b::T, ξ::AbstractArray, y::AbstractArray, p::Integer)
    for i = 1:length(ξ)
        y[i] = basis1d(b, ξ[i], p)
    end
    y
end

basis1d{T<:BasisFun1d}(b::T, ξ::AbstractArray, p::Integer) = basis1d!(b, ξ, similar(ξ), p)

"""
Karniadakis/Sherwin 1D modal polynomial basis.
"""
immutable ModalC01d{T<:Number} <: BasisFun1d{T}
    nmodes::Int
    lnum::LocalNumSys1d
end
ModalC01d{T<:Number}(n, ::Type{T}=Float64) = ModalC01d{T}(n, LocalNumSys1d([1,2], [3:n;]))

isbndryint(b::ModalC01d) = true
isbndryint(::Type{ModalC01d}) = true

function basis1d{T<:Number}(b::ModalC01d, ξ::T, p::Integer)

    if p == 1
      ϕ = (one(T) - ξ) / 2
    elseif p == 2
      ϕ = (one(T) + ξ) / 2
    else
      ϕ = (one(T) - ξ)*(one(T) + ξ) / 4 * jacobi(ξ, p-3, 1, 1)
    end
    return ϕ
end

"""
Legendre orthogonal basis
"""
immutable Legendre1d{T<:Number} <: BasisFun1d{T}
    nmodes::Int
    lnum::LocalNumSys1d
end
Legendre1d(n) = Legendre1d(n+1, LocalNumSys1d( [1:(n+1);], Int[]) )

basis1d{T<:Number}(b::Legendre1d, ξ::T, p::Integer) = legendre(ξ, p-1)

"""
Lagrange polynomial basis.
"""
immutable Lagrange1d{T<:Number} <: BasisFun1d
    nmodes::Int
    lnum::LocalNumSys1d
    z::Array{T,1}
end

"""
Build a Lagrange polynomial basis at Gauss-Lobatto nodes.
"""
function Lagrange1d{T<:Number}(n, ::Type{T}=Float64)
    z = Jacobi.zglj(n, 0, 0, T)
    Lagrange1d{T}(n, LocalNumSys1d([1,n], [2:n-1;]), z)
end

"""
Build a Lagrange polynomial basis at any set of distinct nodes.
"""
function Lagrange1d{T<:Number}(x::AbstractVector{T})
    n = length(x)
    z = zeros(T, n)

    for i = 1:n
        z[i] = x[i]
    end
    bidx = Int[]
    iidx = Int[]

    for i = 1:n
        if z[i]==-1 || z[i]==1
            append!(bidx, [i])
        else
            append!(iidx, [i])
        end
    end
    Lagrange1d(n, LocalNumSys1d(bidx, iidx), z)
end

isbndryint(b::Lagrange1d) = isbndryint(b.lnum)

basis1d{T<:Number}(b::Lagrange1d{T}, ξ::T, p) = Jacobi.lagrange(p, ξ, b.z)


text:
using LinearAlgebra
using OrderedCollections
using IRTools

include("funres.jl")
include("utils.jl")
include("tape.jl")
include("trace.jl")
include("compile.jl")


text:
# Script to test text layout observable updating on a toggle vector.
# Problem:
# I have a vector of 2 passes, 2 toggles and 2 labels, and upon toggling,
# I want for a function that depends on the internal values of
# `passes` to recompute a function and paint that text.
# Solution 1 (did not work):
# Make each pass an observable, so that the toggles/labels work,
# but also make the vector of observable passes an observable,
# so that painting the text is only done once.
# Solution 2 (morning of day after):
# Use the `@lift` macro, interpolate the toggle observables (`toggles[1].active`, that is)
using GLMakie
    #passes = [
    #          (;checked = true, pass = "pass 1"),
    #          (;checked = true, pass = "pass 2")
    #         ]
    passes = ["pass 1", "pass 2"]
    f = Figure(backgroundcolor = :gray70);
    passbox = f[1,2] = Axis(f[1,2], title = "Passes")
    toggles = [Toggle(f), Toggle(f)]
    labels = [Label(f, lift(x -> x ? "ON $(toggles[i].active.val)" : "OFF $(toggles[i].active.val)", toggles[i].active), textsize = 48f0) for i in 1:2]
    f[1,2] = grid!(hcat(toggles, labels), tellheight = false)
    textbox = f[1,1] = Axis(f[1,1], title = "Compiled code")
    #str = @lift "Toggle 1 is " * repr($(toggles[1].active)) * " and Toggle 2 is " * repr($(toggles[2].active))
    str = @lift "Sum of toggles is " * repr($(toggles[1].active) + $(toggles[2].active))
    #text!(textbox, lift(x ->  string(obs_passes[1][].checked + obs_passes[2][].checked), vector_obs_passes))
    text!(textbox, str)
    #text!(textbox, @lift($(toggles[1].active) + $(toggles[2].active)), vector_obs_passes)
    f


text:
using GAIO
using StaticArrays
using Test

@testset "exported functionality" begin
    partition = BoxPartition(Box(SVector(0.0, 0.0, 0.0, 0.0),
                                     SVector(1.0, 1.0, 1.0, 1.0)))
    @testset "basics" begin
        @test dimension(partition) == 4
    end
    @testset "subdivision" begin
        n = 10
        for k in 1:n
            partition = subdivide(partition, (k%3)+1 )
        end
        @test dimension(partition) == 4
    end
    @testset "size" begin
        partition = BoxPartition(Box(SVector(0.0, 1.0), SVector(1.0, 1.0)), (4,2))
        @test size(partition) == (4, 2)
    end
    @testset "domain with zero radius" begin
        center = SVector(0.0, 0.0)
        radius = SVector(1.0, 0.0)
        @test_throws DomainError Box(center, radius)
    end
end
@testset "internal functionality" begin
    partition = BoxPartition(Box(SVector(0.0, 0.0, 0.0), SVector(1.0, 1.0, 1.0)), (4,2,2))
    @testset "keys all" begin
        @test size(GAIO.keys_all(partition)) == (4*2*2,)
    end
    inside = SVector(0.5, 0.5, 0.5)
    left = SVector(-1.0, -1.0, -1.0)
    right = SVector(1.0, 1.0, 1.0)
    on_boundary_left = SVector(0.0, 0.0, -1.0)
    on_boundary_right = SVector(0.0, 1.0, 0.0)
    outside_left = SVector(0.0, 0.0, -2.0)
    outside_right = SVector(0.0, 2.0, 0.0)
    @testset "point to key" begin
        @test !isnothing(GAIO.point_to_key(partition, inside))
        @test !isnothing(GAIO.point_to_key(partition, left))
        @test isnothing(GAIO.point_to_key(partition, right))
        @test !isnothing(GAIO.point_to_key(partition, on_boundary_left))
        @test isnothing(GAIO.point_to_key(partition, on_boundary_right))
        @test isnothing(GAIO.point_to_key(partition, outside_left))
        @test isnothing(GAIO.point_to_key(partition, outside_right))
        key = GAIO.point_to_key(partition, partition.domain.center)
        @test typeof(key) <: GAIO.keytype(typeof(partition))
    end
    key_inside = GAIO.point_to_key(partition, inside)
    key_left = GAIO.point_to_key(partition, left)
    @testset "key to point" begin
        @test typeof(GAIO.key_to_box(partition, key_inside)) <: typeof(partition.domain)
        @test GAIO.key_to_box(partition, key_inside) != GAIO.key_to_box(partition, key_left)
    end
    @testset "roundtrip" begin
        point = SVector(0.3, 0.3, 0.3)
        key = GAIO.point_to_key(partition, point)
        box = GAIO.key_to_box(partition, key)
        @test point ∈ box
        key_2 = GAIO.point_to_key(partition, box.center)
        @test key == key_2
    end
    @testset "points with wrong dimension" begin
        point_2d = SVector(0.0, 0.0)
        point_4d = SVector(0.0, 0.0, 0.0, 0.0)
        @test_throws Exception GAIO.point_to_key(partition, point_2d)
        @test_throws Exception GAIO.point_to_key(partition, point_4d)
    end
    @testset "non existing keys" begin
        @test_throws Exception GAIO.key_to_box(partition, -1)
        @test_throws Exception GAIO.key_to_box(partition, 2^partiton.depth + 1)
    end
end


text:
# Initialisations
d = Dict() # Dict of Any -> Any
dint = Dict{Int64, Int64}() # Dict of Int64 -> Int64
dstr = Dict("one" => "uno", "two" => "dos", "three" => "tres")
eg1 = Dict([('a', 1), ('c', 3), ('b', 2)])
eg2 = Dict(zip("def", 4:6))
println("eg1: $eg1")
println("eg2: $eg2")

d["abcdef"] = 3
println(d)
dint[7] = 1
println(dint)
# These are not valid:
# dint["abc"] = 3
# dint[11] = "def"

testdel = Dict([(1, 2), (3, 5), (11, 4), (0, 6)])
delete!(testdel, 3)
println("testdel: $(testdel)")
delete!(testdel, 3) # Still valid even if it does not exist
println("testdel: $(testdel)")

println("dstr[\"one\"]: $(dstr["one"])")
println("dint[7]: $(dint[7])")

println("length(dstr): $(length(dstr))")

println("\"one\" ∈ keys(dstr): $("one" ∈ keys(dstr))") # Checks existence of key
println("haskey(dstr, \"one\"): $(haskey(dstr, "one"))") # Also checks existence of key
println("\"dos\" ∈ values(dstr): $("dos" ∈ values(dstr))") # Checks existence of value
println("Pair(\"one\", \"uno\") ∈ dstr: $(Pair("one", "uno") in dstr)") # Checks existence of pair

function printdict(h)
    for c in keys(h)
        println(c, " -> ", h[c])
    end
end

# printdict(dstr)

# Alternatively:
function keyvaluepairs(h)
    for (k, v) in h
        println("$k -> $v")
    end
end

keyvaluepairs(dstr)

function printdictsorted(h)
    for c in sort(collect(keys(h)))
        println(c, " -> ", h[c])
    end
end

# printdictsorted(dstr)

# Reverse lookup
println(findall(isequal("uno"), dstr))


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

using Markdown
using InteractiveUtils

# ╔═╡ 16ddb41a-fc59-11ea-1631-153e3466c75c
using Pkg, DrWatson

# ╔═╡ d65dd2b2-fc58-11ea-2300-4db47ec9a789
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ e875dcfc-fc57-11ea-27e5-c56f1f9d5370
md"## clip-05-01-05s.jl"

# ╔═╡ d65e1360-fc58-11ea-11c0-5928313bb9a0
md"### snippet 5.1"

# ╔═╡ cfa44fec-01c5-11eb-14bf-338eed7e2c9d
md"##### Notice that in below Stan language model we ignore the observed data (the likelihood is commented out). The draws show sampled regression lines implied by the priors. This is an alternative to the formulation chosen in clip-04-01-05s.jl"

# ╔═╡ d65e98dc-fc58-11ea-25e1-9fab97b6125a
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame; delim=';');
	scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
end;

# ╔═╡ d66f515e-fc58-11ea-3fae-cbb82f1a1a6a
stan5_1_alt_priors = "
	data {
	 int < lower = 1 > N; // Sample size
	 vector[N] D; // Outcome
	 vector[N] A; // Predictor
	}

	parameters {
	 real a; // Intercept
	 real bA; // Slope (regression coefficients)
	 real < lower = 0 > sigma;    // Error SD
	}

	model {
	  vector[N] mu;               // mu is a vector
	  a ~ normal(0, 0.2);         // Priors
	  bA ~ normal(0, 0.5);
	  sigma ~ exponential(1);
	  mu = a + bA * A;
	  //D ~ normal(mu , sigma);   // Likelihood
	}
";

# ╔═╡ f4602d4a-fc59-11ea-0d9d-9f58c73c119f
md"### snippet 5.3 - 5.4"

# ╔═╡ d670aefa-fc58-11ea-1c56-4bfb66e1cab2
md"## Define the SampleModel, etc."

# ╔═╡ d67e0602-fc58-11ea-3a27-31d03e1c2318
begin
	m5_1s = SampleModel("MedianAgeMarriage", stan5_1_alt_priors)
	m5_1_data = Dict("N" => size(df, 1), "D" => df.Divorce_s, "A" => df.MedianAgeMarriage_s)
	rc5_1s = stan_sample(m5_1s, data=m5_1_data)
	success(rc5_1s) && (post5_1s_df = read_samples(m5_1s, :dataframe))
end;

# ╔═╡ a4a9351a-01c6-11eb-28d0-71f8fb243719
PRECIS(post5_1s_df)

# ╔═╡ 12fedbca-fc5a-11ea-2d4d-1d5ac93ac4fa
md"### snippet 5.5"

# ╔═╡ 45b2b002-01c6-11eb-3f86-3f9586afcc8b
md"##### Plot priors of the intercept (`:a`) and the slope (`:bA`)."

# ╔═╡ d68ab980-fc58-11ea-342d-31e66a8e7559
if success(rc5_1s)
	xi = -3.0:0.1:3.0
	plot(xlab="Medium age marriage (scaled)", ylab="Divorce rate (scaled)",
		title="Showing 50 regression lines")
	for i in 1:50
		local yi = mean(post5_1s_df[i, :a]) .+ post5_1s_df[i, :bA] .* xi
		plot!(xi, yi, color=:lightgrey, leg=false)
	end
	scatter!(df[:, :MedianAgeMarriage_s], df[!, :Divorce_s], color=:darkblue)
end

# ╔═╡ d69533ba-fc58-11ea-3378-e512a1d55d27
md"## End of clip-05-01-05.jl"

# ╔═╡ Cell order:
# ╟─e875dcfc-fc57-11ea-27e5-c56f1f9d5370
# ╠═16ddb41a-fc59-11ea-1631-153e3466c75c
# ╠═d65dd2b2-fc58-11ea-2300-4db47ec9a789
# ╟─d65e1360-fc58-11ea-11c0-5928313bb9a0
# ╟─cfa44fec-01c5-11eb-14bf-338eed7e2c9d
# ╠═d65e98dc-fc58-11ea-25e1-9fab97b6125a
# ╠═d66f515e-fc58-11ea-3fae-cbb82f1a1a6a
# ╟─f4602d4a-fc59-11ea-0d9d-9f58c73c119f
# ╟─d670aefa-fc58-11ea-1c56-4bfb66e1cab2
# ╠═d67e0602-fc58-11ea-3a27-31d03e1c2318
# ╠═a4a9351a-01c6-11eb-28d0-71f8fb243719
# ╟─12fedbca-fc5a-11ea-2d4d-1d5ac93ac4fa
# ╟─45b2b002-01c6-11eb-3f86-3f9586afcc8b
# ╠═d68ab980-fc58-11ea-342d-31e66a8e7559
# ╟─d69533ba-fc58-11ea-3378-e512a1d55d27


text:
s1 = Statement("a=2", "a")

sdef1 = SDefinition("sdef1", false)
sdef2 = SDefinition("sdef2", true)
sdef3 = SDefinition("sdef3", true, fields=("device" => Ptr{Nothing}, "notype" => nothing))

signature = Signature(first(methods(replace)))
fdef1 = FDefinition("replace", signature, false, Statement[])
fdef2 = FDefinition("replace", signature, false, [Statement("a=5")])
fdef3 = FDefinition("replace", signature, false, Statement.(["i=1", "i+=1", "nothing"]))

cdef1 = CDefinition("c1", 0)
cdef2 = CDefinition("c2", "[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]")

edef1 = EDefinition(name="e1", fields=["a", "b", "c", "d"])
edef2 = EDefinition(name="e1", fields=["a", "b", "c", "d"], with_begin_block=true, enum_macro="@enum")

@testset "Code generation" begin
    @testset "Structure definition" begin
        @test generate(sdef1) == format_text("struct sdef1 end")
        @test generate(sdef2) == format_text("mutable struct sdef2 end")
        @test generate(sdef3) == format_text("""
                         mutable struct sdef3
                             device::Ptr{Nothing}
                             notype::Any
                         end""")
    end

    @testset "Function definition" begin
        @test generate(fdef1) == format_text("function replace(a::Union{Function, Type}, b::Pair; count) end")
        @test generate(fdef2) == format_text("function replace(a::Union{Function, Type}, b::Pair; count) a=5 end")
        @test generate(fdef3) == format_text("""
                         function replace(a::Union{Function, Type}, b::Pair; count)
                             i = 1
                             i += 1
                             nothing
                         end""")
        @test_throws ErrorException generate(FDefinition("replace", signature, true, Statement[])) # short notation for a function without body
    end
    @testset "Statements" begin
        @test generate([s1]) == format_text("a = 2")
    end
    @testset "Constant definition" begin
        @test generate(cdef1) == format_text("const c1 = 0")
        @test generate(cdef2) == format_text("const c2 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]")
    end
    @testset "Enum definition" begin
        @test generate(edef1) == format_text("@enum e1 a b c d")
        @test generate(edef2) == format_text("@enum e1 begin\na\nb\nc\nd\nend")
    end
end


text:
baremodule colors
# FC is fill color and BC is border color

# Brush outlines and fill
const brush_bc = (1.0, 0.5, 0.0, 0.5)
const brush_fc = (1.0, 0.5, 0.0, 0.2)

# The selection tool selection rectangle
const selection_selection_fc = (0.0, 1.0, 0.0, 0.6)
const selection_selection_bc = (0.0, 1.0, 0.0, 0.4)

# Selected items with selection tool
const selection_selected_fc = (1.0, 0.0, 1.0, 0.6)
const selection_selected_bc = (1.0, 0.0, 1.0, 0.4)

# Selection preview with selection tool
const selection_preview_fc = (1.0, 1.0, 0.0, 0.6)
const selection_preview_bc = (1.0, 1.0, 0.0, 0.4)

# Trigger rendering colors
const trigger_fc = (0.6, 0.2, 0.2, 0.6)
const trigger_bc = (0.6, 0.2, 0.2, 0.4)

# Entity/decal placement movement preview
const ghost_preplacement_alpha = 0.5

# Alpha of unselected rooms
const adjacent_room_alpha = 0.5

# Default background fills
const background_room_fill = (0.2, 0.2, 0.2, 1.0)
const background_canvas_fill = (0.15, 0.15, 0.15, 1.0)

# Color coded background fills
const background_room_color_coded_fill = [
    (0.20, 0.20, 0.20, 1.0),
    (0.24, 0.20, 0.20, 1.0),
    (0.20, 0.24, 0.20, 1.0),
    (0.20, 0.21, 0.22, 1.0),
    (0.20, 0.20, 0.25, 1.0),
    (0.24, 0.20, 0.25, 1.0),
    (0.22, 0.20, 0.20, 1.0)
]

# Filler room color
const filler_room_fill = (0.1, 0.4, 0.3, 1.0)

# Unknown tile color
const unknown_tile_color = (1.0, 0.0, 1.0, 1.0)

# Canvas Font color
const canvas_font_color = (0.0, 0.0, 0.0, 1.0)

end


text:
# Records the statistics of two state variables X and Y:
# - covariance cov(X,Y)
# - correlation coefficient
# - marginal variance var(X) and var(Y)
# - joint average E(XY)
# - marginal average E(X) and E(Y)

# Uses the shifted-data method to avoid loss of precision.

mutable struct JointStatsRecorder{T,X,XF,Y,YF}
    x::X
    y::Y
    Kx::Vector{T}
    Ky::Vector{T}
    Σx::Vector{T}
    Σy::Vector{T}
    Σxx::Vector{T}
    Σyy::Vector{T}
    Σxy::Matrix{T}
    n::Int
end

function JointStatsRecorder(x, field_x::Symbol, y, field_y::Symbol)
    data_x = getfield(x, field_x)
    data_y = getfield(y, field_y)
    T = promote_type(eltype(data_x), eltype(data_y))
    X = typeof(x)
    Y = typeof(y)
    XF = field_x
    YF = field_y
    JointStatsRecorder{T,X,XF,Y,YF}(x, y,
        zeros(T, length(data_x)),
        zeros(T, length(data_y)),
        zeros(T, length(data_x)),
        zeros(T, length(data_y)),
        zeros(T, length(data_x)),
        zeros(T, length(data_y)),
        zeros(T, (length(data_x), length(data_y))),
        0)
end

function covariance(rc::JointStatsRecorder)
    return (rc.Σxy .- (rc.Σx.*rc.Σy')/rc.n) / rc.n
end

function variances(rc::JointStatsRecorder)
    varx = (rc.Σxx .- (rc.Σx.^2)/rc.n) / rc.n
    vary = (rc.Σyy .- (rc.Σy.^2)/rc.n) / rc.n
    return (varx, vary)
end

function corrcoef(rc::JointStatsRecorder)
    cov = covariance(rc)
    varx, vary = variances(rc)
    return cov ./ (sqrt(varx).*sqrt(vary'))
end

function marginal_averages(rc::JointStatsRecorder)
    Ex = rc.Σx / rc.n .+ rc.Kx
    Ey = rc.Σy / rc.n .+ rc.Ky
    return (Ex, Ey)
end

function sta(rc::JointStatsRecorder)
    rc.Σxy / rc.n .+ rc.Kx .+ rc.Ky
end

function joint_average(rc::JointStatsRecorder)
    Ex, Ey = marginal_averages(rc)
    return covariance(rc) .+ (Ex .* Ey')
end

function reset!(rc::JointStatsRecorder)
    rc.Kx .= 0
    rc.Σx .= 0
    rc.Σxx .= 0
    rc.Ky .= 0
    rc.Σy .= 0
    rc.Σyy .= 0
    rc.Σxy .= 0
    rc.n = 0
end

@generated function record!{T,X,XF,Y,YF}(rc::JointStatsRecorder{T,X,XF,Y,YF}, step)
    quote
        data_x = rc.x.$XF
        data_y = rc.y.$YF

        if rc.n == 0
            rc.Kx .= data_x
            rc.Ky .= data_y
        end

        rc.n += 1

        @inbounds @simd for i in eachindex(data_x)
            dx = data_x[i] - rc.Kx[i]
            rc.Σx[i] += dx
            rc.Σxx[i] += dx^2
        end
        @inbounds @simd for j in eachindex(data_y)
            dy = data_y[j] - rc.Ky[j]
            rc.Σy[j] += dy
            rc.Σyy[j] += dy^2
        end
        @inbounds for i in eachindex(data_x)
            @simd for j in eachindex(data_y)
                dx = data_x[i] - rc.Kx[i]
                dy = data_y[j] - rc.Ky[j]
                rc.Σxy[i,j] += dx*dy
            end
        end
    end
end


text:
using MixedBags
const MB = MixedBags
using Base.Test

@MixedBag Bar Int64 Float64

@testset "Test Internals" begin
    @test MB.firstindex((x -> x=="b"), ["a", "b", "c"]) == 2
    @test MB.fieldname(Int, eltypes(Bar)) == Symbol("_1")
end

@testset "@MixedBag Bar Int64 Float64" begin
    @test eltypes(Bar) == [Int64, Float64]

    @testset "Empty constructor" begin
        @test eltypes(Bar()) == [Int64, Float64]
        @test Bar()[Int64] == Int64[]
        @test Bar()[Float64] == Float64[]
        acc = 0.
        foreach(Bar()) do x
            acc += x
        end
        @test acc == 0.
    end

    b = Bar([1, 2], [0.3])
    @test b == Bar(1, 2, 0.3)
    @test b != Bar(1,2,3)
    @test push!(b, 3) == b
    @test b == Bar(1,2,3, 0.3)
    @test b[Int64] == [1,2,3]
    @test b[Float64] == [0.3]
    acc = 0.
    foreach(b) do x
        acc += x
    end
    @test acc ≈ 1+2+3+0.3

    @testset "mapreduce" begin
        f = √
        op = +
        v0 = 42
        bcoll = [1,2,3,0.3]
        @test bcoll == collect(b)
        @test mapreduce(f, op, v0, bcoll) ≈ mapreduce(f, op, v0, b)
        @test mapreduce(f, op, bcoll)     ≈ mapreduce(f, op, b)
        @test reduce(op, bcoll)           ≈ reduce(op, b)
        @test reduce(op, v0, bcoll)       ≈ reduce(op, v0, b)

        @test append!(b, [50, 60]) == b
        @test b == Bar(1,2,3,50,60, 0.3)
    end

    # @test Bar(1, 2) == Bar(1, 2)
end

typealias MyInt Int
@MixedBag ComplicatedBag Matrix{String} Matrix{MyInt}

@testset "Bag with compicated types" begin
    cb = ComplicatedBag()
    @test cb == ComplicatedBag()
    @test push!(cb, ["a" "b"; "c" "d"]) == ComplicatedBag(Array{String,2}[String["a" "b"; "c" "d"]],Array{Int64,2}[])

end


text:
@safetestset "ChebyshevImplementation" begin include("test_ChebyshevImplementation.jl") end

@safetestset "SecondPoincareInvariant OOP" begin
    using PoincareInvariants

    @testset "$N Points on Square in $(D)D" for D in [2, 10], N in [10, 123, 4321]
        Ω(z, t, p) = CanonicalSymplecticMatrix(D)
        pinv = SecondPoincareInvariant{Float64}(Ω, D, N, Val(false))

        phasepoints = getpoints(pinv) do x, y
            ntuple(D) do i
                i == 1 && return x
                i == D ÷ 2 + 1 && return y
                return 0
            end
        end

        @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 10
    end
end

@safetestset "Free Particle" begin
    using PoincareInvariants

    function free_particle!(state, δt)
        mid = length(state) ÷ 2
        state[1:mid] += state[mid+1:end] .* δt
    end

    @testset "$N Points on Square in $(D)D" for D in [2, 6], N in [500, 10_000]
        Ω(z, t, p) = CanonicalSymplecticMatrix(D)
        pinv = SecondPoincareInvariant{Float64}(Ω, D, N, Val(false))

        parampoints = getpoints(pinv)

        phasepoints = getpoints(pinv) do x, y
            ntuple(D) do i
                i == 1 && return x
                i == D ÷ 2 + 1 && return y
                return 0
            end
        end

        @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 10

        free_particle!.(eachrow(phasepoints), 10)
        @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 20

        free_particle!.(eachrow(phasepoints), 100)
        @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 100

        free_particle!.(eachrow(phasepoints), 1000)
        @test abs(1 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 1_000
    end

    @testset "$N Points on Quarter Circle in $(D)D " for D in [4, 12], N in [234, 5678]
        Ω(z, t, p) = CanonicalSymplecticMatrix(D)
        pinv = SecondPoincareInvariant{Float64}(Ω, D, N, Val(false))

        parampoints = getpoints(pinv)

        phasepoints = getpoints(pinv) do r, θ
            ntuple(D) do i
                i == 1 && return r * cos(θ * π/2)
                i == D ÷ 2 + 1 && return r * sin(θ * π/2)
                return 0
            end
        end

        @test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 10

        free_particle!.(eachrow(phasepoints), 10)
        @test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 20

        free_particle!.(eachrow(phasepoints), 100)
        @test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 200

        free_particle!.(eachrow(phasepoints), 1000)
        @test abs(π/4 - compute!(pinv, phasepoints, 0, nothing)) / eps() < 2_000
    end
end


text:
# fastest sorting by decreasing magnitude
@inline amaxmin(x::T, y::T) where {T} = ifelse( abs(x) < abs(y), (y,x), (x,y) )

@inline function amaxmin(x::T, y::T, z::T) where {T}
    y, z = amaxmin(y, z)
    x, z = amaxmin(x, z)
    x, y = amaxmin(x, y)
    return x, y, z
end

"""
    two_sum(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_sum(a::T, b::T) where {T<:Real}
    hi = a + b
    v  = hi - a
    lo = (a - (hi - v)) + (b - v)
    return hi, lo
end

"""
    two_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
"""
@inline function two_sum(a::T, b::T, c::T) where {T}
    hi, md, lo = amaxmin(a, b, c)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    md += lo
    return hi, md
end

"""
   three_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `md = err(a+b+c), lo = err(md)`.
"""
function three_sum(a::T, b::T, c::T) where {T}
    hi, md, lo = amaxmin(a, b, c)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    md, lo = two_hilo_sum(md, lo)
    return hi, md, lo
end

"""
    two_diff(a, b)

Computes `s = fl(a-b)` and `e = err(a-b)`.
"""
@inline function two_diff(a::T, b::T) where {T}
    hi = a - b
    v  = hi - a
    lo = (a - (hi - v)) - (b + v)
    return hi, lo
end

"""
    three_diff(a, b, c)

Computes `s = fl(a-b-c)` and `e1 = err(a-b-c), e2 = err(e1)`.
"""
function three_diff(a::T,b::T,c::T) where {T}
    s, t = two_diff(-b, c)
    x, u = two_sum(a, s)
    y, z = two_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    two_prod(a, b)

Computes `hi = fl(a*b)` and `lo = fl(err(a*b))`.
"""
@inline function two_prod(a::T, b::T) where {T}
    hi = a * b
    lo = fma(a, b, -hi)
    hi, lo
end

"""
    two_prod(a, b, c)

Computes `hi = fl(a*b*c)` and `lo = err(hi)`.
"""
function two_prod(a::T, b::T, c::T) where {T}
    t,  k = two_prod(a, b)
    hi, e = two_prod(t, c)
    md, lo = two_fma(c, k, e)
    hi, md = two_hilo_sum(hi, md)
    md += lo
    return hi, md
end

"""
    three_prod(a, b, c)

Computes `hi = fl(a*b*c)` and `md = err(a*b*c), lo = err(md)`.
"""
function three_prod(a::T, b::T, c::T) where {T}
    t,  k = two_prod(a, b)
    hi, e = two_prod(t, c)
    md, lo = two_fma(c, k, e)
    hi, md = two_hilo_sum(hi, md)
    md, lo = two_hilo_sum(md, lo)
    return hi, md, lo
end

#=
   three_fma algorithm from
   Sylvie Boldo and Jean-Michel Muller
   Some Functions Computable with a Fused-mac
=#

"""
   two_fma(a, b, c)

Computes `hi = fl(fma(a,b,c))` and `lo = fl(err(fma(a,b,c)))`.
"""
function two_fma(a::T, b::T, c::T) where {T}
    hi = fma(a, b, c)
    hi0, lo0 = two_prod(a, b)
    hi1, lo1 = two_sum(c, lo0)
    hi2, lo2 = two_sum(hi0, hi1)
    lo = ((hi2 - hi) + lo2) + lo1
    return hi, lo
end

"""
   three_fma(a, b, c)

Computes `hi = fl(fma(a,b,c))` and `md = fl(err(fma(a,b,c))), lo = fl(err(md))`.
"""
function three_fma(a::T, b::T, c::T) where {T}
    hi = fma(a, b, c)
    hi0, lo0 = two_prod(a, b)
    hi1, lo1 = two_sum(c, lo0)
    hi2, lo2 = two_sum(hi0, hi1)
    y = (hi2 - hi) + lo2
    md, lo = two_hilo_sum(y, lo1)
    return hi, md, lo
end

# with arguments sorted by magnitude

"""
    two_hilo_sum(a, b)

*unchecked* requirement `|a| ≥ |b|`

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_hilo_sum(a::T, b::T) where {T}
    hi = a + b
    lo = b - (hi - a)
    return hi, lo
end

"""
    two_hilo_sum(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
"""
function two_hilo_sum(a::T, b::T, c::T) where {T}
    lo, t = two_hilo_sum(b, c)
    hi, lo = two_sum(a, lo)
    lo += t
    hi, lo = two_hilo_sum(hi, lo)
    return hi, lo
end

"""
    three_hilo_sum(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`

Computes `x = fl(a+b+c)` and `y = err(a+b+c), z = err(y)`.
"""
function three_hilo_sum(a::T, b::T, c::T) where {T}
    md, lo = two_hilo_sum(b, c)
    hi, md = two_sum(a, md)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    return hi,md,lo
end

"""
    two_hilo_diff(a, b)

*unchecked* requirement `|a| ≥ |b|`

Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_hilo_diff(a::T, b::T) where {T}
    hi = a - b
    lo = (a - hi) - b
    hi, lo
end

"""
    three_hilo_diff(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`

Computes `x = fl(a-b-c)` and `y = err(a-b-c), z = err(y)`.
"""
function three_hilo_diff(a::T,b::T,c::T) where {T}
    md, lo = two_hilo_sum(-b, -c)
    hi, md = two_sum(a, md)
    md, lo = two_hilo_sum(md, lo)
    hi, md = two_hilo_sum(hi, md)
    return hi,md,lo
end

#=
"Concerning the division [and sqrt], the elementary rounding error is
generally not a floating point number, so it cannot be computed
exactly. Hence we cannot expect to obtain an error
free transformation for the division. ...
This means that the computed approximation is as good as
we can expect in the working precision."
-- http://perso.ens-lyon.fr/nicolas.louvet/LaLo05.pdf

"While the sqrt algorithm [and division] is not strictly an errorfree transformation,
it is known to be reliable and is recommended for general use."
-- https://hal-ens-lyon.archives-ouvertes.fr/ensl-00545591v2/document
=#

@inline function two_div(a::T, b::T) where {T}
     hi = a / b
     lo = fma(-hi, b, a)
     lo /= b
     return hi, lo
end

@inline function two_sqrt(a::T) where {T}
    hi = sqrt(a)
    lo = fma(-hi, hi, a)
    lo /= 2
    lo /= hi
    return hi, lo
end


text:
module TestSearchCv

using Test
using LightGBM

train_matrix = rand(1000,70) # create random dataset
train_labels = rand([0, 1], 1000)
splits = (
    collect(1:200), collect(201:400),
    collect(401:600), collect(601:800),
    collect(801:1000)
)
params = [
    Dict(:num_iterations => num_iterations, :num_leaves => num_leaves)
    for num_iterations in (1, 2), num_leaves in (5, 10)
]

estimator = LightGBM.LGBMClassification(
    objective = "binary",
    num_class = 1,
    is_training_metric = true,
    metric = ["auc"]
)

@testset "test search_cv with X y input -- binary" begin
    # Act
    output = LightGBM.search_cv(estimator, train_matrix, train_labels, splits, params, verbosity = -1);

    # Assert
    @test output isa Array{Tuple{Dict{Symbol,Any},Dict{String,Dict{String,Vector{Float64}}}}}
    @test :num_iterations in keys(output[1][1])
    @test :num_leaves in keys(output[1][1])
    @test length(output[1][2]["validation"]["auc"]) == 5
    @test length(output[1][2]["training"]["auc"]) == 5
end

@testset "test search_cv with Dataset input -- binary" begin
    # Arrange
    train_dataset = LightGBM.LGBM_DatasetCreateFromMat(train_matrix, "")
    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)

    # Act
    output = LightGBM.search_cv(estimator, train_dataset, splits, params, verbosity = -1);

    # Assert
    @test output isa Array{Tuple{Dict{Symbol,Any},Dict{String,Dict{String,Vector{Float64}}}}}
    @test :num_iterations in keys(output[1][1])
    @test :num_leaves in keys(output[1][1])
    @test length(output[1][2]["validation"]["auc"]) == 5
    @test length(output[1][2]["training"]["auc"]) == 5
end

end # module


text:
function validate_available_devices(devices, device::Type{D}) where {D<:PSY.Device}

    if isempty(devices)
        @warn("The data doesn't include devices of type $(device), consider changing the device models")
        return true
    end

    return false

end


text:
using Distributed

## Function to get the number of points inside a circle
function in_circle_points(n)
    points_in = 0 ## Points inside the circle
    for i = 1:n
        x, y = rand(), rand() ## Randomly generate x,y btwn 0 & 1
        ## Count if generated value is within radius of the circle
        ## int the unit square
        points_in += (x*x + y*y) <= 1
    end
    ## Return the number of points in
    return points_in
end

## Run code in parallel across n cores
function parallel_pi(n)
    ## Create the number of workers specified by the
    p = nworkers()
    ## Combine results via addition (+)
    n_in = @distributed (+) for i=1:p
        in_circle_points(n / p)
    end
    return 4 * n_in / n ## Approximation of pi
end

## Run code serially on 1 worker
function serial_pi(n)
    n_in = in_circle_points(n)
    return 4 * n_in / n
end