text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Initialize node_indices of MPI interface container
@inline function init_mpi_interface_node_indices!(mpi_interfaces::P4estMPIInterfaceContainer{2},
                                                  faces, local_side, orientation,
                                                  mpi_interface_id)
  # Align interface in positive coordinate direction of primary element.
  # For orientation == 1, the secondary element needs to be indexed backwards
  # relative to the interface.
  if local_side == 1 || orientation == 0
    # Forward indexing
    i = :i_forward
  else
    # Backward indexing
    i = :i_backward
  end

  if faces[local_side] == 0
    # Index face in negative x-direction
    mpi_interfaces.node_indices[mpi_interface_id] = (:begin, i)
  elseif faces[local_side] == 1
    # Index face in positive x-direction
    mpi_interfaces.node_indices[mpi_interface_id] = (:end, i)
  elseif faces[local_side] == 2
    # Index face in negative y-direction
    mpi_interfaces.node_indices[mpi_interface_id] = (i, :begin)
  else # faces[local_side] == 3
    # Index face in positive y-direction
    mpi_interfaces.node_indices[mpi_interface_id] = (i, :end)
  end

  return mpi_interfaces
end

end # muladd

ext: jl


text:
"""
    uipathorchestratorschedulreadjustment(scheduleplan::DataFrame,robotn::Int,run_unit_time::Int,jobn::Int,timen::Int,schedulcolumn::Int=6)

# 処理概要
スケジュール作成、処理結果の整合性確認および整形

# 引数
* `scheduleplan`:スケジュール作成のために必要な諸元
* `robotn`:同時最大稼働ロボット数
* `run_unit_time`:実行単位時間、コマ割時間
* `jobn`:ジョブ数
* `timen`:時間コマ数
* `schedulcolumn`:スケジュールON,OFF開始列を指定（デフォルト：６）

# 結果（戻り値）
スケジュール調整結果を出力
"""
function uipathorchestratorschedulreadjustment(scheduleplan::DataFrame,robotn::Int,run_unit_time::Int,jobn::Int,timen::Int;schedulcolumn::Int=6)
## 数理モデル
 m1 = JuMP.Model(Ipopt.Optimizer)

  ## 変数
 JuMP.@variable(m1, 0<=s[1:jobn,1:timen] <=1 ) #ジョブ毎のコマ割りごとに実行有無を表示

  ## 定数
  runtime=zeros(Int,jobn,1) #ジョブ毎の実行時間

  ## ジョブ実行時間
  ###割り当てるコマ数を算出
  for i in 1:jobn
    runtime[i]=scheduleplan[i,:runtime]/run_unit_time
  end

  ## 目的関数
  JuMP.@objective(m1, Min, sum(scheduleplan[1:end,:runtime])/run_unit_time-sum(s[1:jobn,1:timen] ))

  ## 制約条件
  ### ジョブ実行時間予約指定
  for i in 1:jobn
    flag1 =true
    if scheduleplan[i,:Specifiedtime] == 1
      for j in schedulcolumn:schedulcolumn+timen-1
        if typeof(scheduleplan[i,j]) !=Missing && flag1
          index=j-schedulcolumn+1
          JuMP.@constraint(m1, s[i,index:index+runtime[i]-1] .== 1.0 )
           flag1=false
        end
      end
    elseif scheduleplan[i,:Specifiedtime] == 2
      df = DateFormat("HH:MM");
      StartTime=DateTime(scheduleplan[i,:JobStartTime],df)
      EndTime=DateTime(scheduleplan[i,:JobEndTime],df)
      for j in schedulcolumn:schedulcolumn+timen-1
        if DateTime(string(names(scheduleplan)[j]),df) >= StartTime && DateTime(string(names(scheduleplan)[j]),df) < EndTime
        else
          JuMP.@constraint(m1, s[i,j-schedulcolumn+1] == 0.0 )
        end
      end
    elseif scheduleplan[i,:Specifiedtime] == 0 || ismissing(scheduleplan[i,:Specifiedtime])
    else
      println(i,"行目の指定パターンを確認してください")
    end
  end

  ### ジョブ実行時間制限
  for i in 1 :jobn
    JuMP.@constraint(m1,sum(s[i,1:timen])==Float64(runtime[i]))
  end

  ### ロボット同時実行制限(ロボット数を超過させない)
  for i in 1 :timen
    JuMP.@constraint(m1,sum(s[1:jobn,i])<=Float64(robotn) )
  end

  ### ジョブ連続実行制限
  JuMP.register(m1, :ContinuousOperation, size(s)[2], ContinuousOperation; autodiff = true)

  for i in 1:jobn
    x=s[i,1:timen]
    @constraint(m1,sum(ContinuousOperation(x...))== Float64(runtime[i]-1) )
  end

  ## ソルバーの実行
  status = JuMP.optimize!(m1)

  ## スケジュール案表示
  plan=zeros(Int,jobn,timen)
  plan=map(Int,map(round,JuMP.value.(s)))

  # 割当時間1コマの不具合対策
  if(sum(plan) != sum(runtime))
    schedulesubplan=scheduleplan[:,schedulcolumn:end ]
    schedulesubplan[:,:].= hcat(plan,zeros(Int,jobn,1))
    schedulesubplan=hcat(scheduleplan[:,1:schedulcolumn-1 ] , schedulesubplan)
    for i in 1:jobn
      if(sum(schedulesubplan[i,schedulcolumn:end]) == runtime[i])
        schedulesubplan[i,:Specifiedtime] =1
      end
    end

    plan = uipathorchestratorschedulreadjustmentsub1(schedulesubplan,robotn,run_unit_time,jobn,timen,schedulcolumn=schedulcolumn)
  end

  return plan,runtime
end

"""
    adjustedresultcheck(plan::Array,runtime::Array,scheduleplan::DataFrame)

# 処理概要
スケジュール調整した結果が正しく結果かどうかチェック。作成されたスケジュールの妥当性チェック

# 引数
* `plan`:uipathorchestratorschedulreadjustmentの実行結果を指定
* `runtime`:readprerequisiteで読み込んだ結果を指定
* `scheduleplan`:スケジュール作成のために必要な諸元

# 結果（戻り値）
スケジュール調整結果を出力
"""
function adjustedresultcheck(plan::Array,runtime::Array,scheduleplan::DataFrame,robotn::Int,jobn::Int,timen::Int;schedulcolumn::Int=6,checkreturn::Bool=false)
adjustedresultcheckflag=true
adjustedresultcheckmastarflag=Array{Bool}(undef,jobn,3) # #1:ジョブごとに実行時間確保されているかチェック,#2:ロボット数超過,#3:処理時間が連続しているか
adjustedresultcheckmastarflag .= true

  #ジョブごとに実行時間確保されているかチェック
  for i in 1:jobn
    adjustedresultcheckflag1=true
    for j in 1:timen
      if(plan[i,j]==1 && adjustedresultcheckflag1)
        adjustedresultcheckflag1=false
        if(sum(plan[i,j:(j+runtime[i]-1)])!=runtime[i])
          adjustedresultcheckmastarflag[i,1]=false
        end
      end
    end
  end

  # ロボット数超過していないか確認
  for i in 1:timen
    if sum(plan[:,i]) > robotn
      adjustedresultcheckmastarflag[:,2] .=false
    end
  end

  # 処理時間が連続しているか確認
  for i in 1:jobn
    x=plan[i,:]
    if(sum(ContinuousOperation(x...))!=runtime[i]-1)
      adjustedresultcheckmastarflag[i,3]=false
    end
  end

  if sum(adjustedresultcheckmastarflag) != size(adjustedresultcheckmastarflag,1)*size(adjustedresultcheckmastarflag,2)
    plan=zeros(Int,jobn,timen)
    adjustedresultcheckflag = false
  end

  result=scheduleplan[:,schedulcolumn:end ]
  result[:,:].= hcat(plan,zeros(Int,jobn,1))
  result=hcat(scheduleplan[:,1:schedulcolumn-1 ] , result)

  #処理開始時間および終了時間を取得
  times=names(scheduleplan)

  for i in 1 :  size(scheduleplan,1)
    startflag=true
    endflag=true

    for j in schedulcolumn : size(scheduleplan,2)
      if(result[i,j] == 1 && startflag)
        startflag = false
        result[i,:JobStartTime]=string(times[j])
      end

      if result[i,j] == 0 && !startflag && endflag
        endflag = false
        result[i,:JobEndTime]=string(times[j])
      end
    end
  end

  result[:,end] .=missing

  # 調整が失敗した場合明示的にわかるように、JobStartTime,JobEndTimeを空白にする
  if !adjustedresultcheckflag
    result[:,:JobStartTime] .= missing
    result[:,:JobEndTime].= missing
  end

  if checkreturn
    return result,adjustedresultcheckflag
  else
    return result
  end

end

# スケジュール作成サブ１
function uipathorchestratorschedulreadjustmentsub1(schedulesubplan::DataFrame,robotn::Int,run_unit_time::Int,jobn::Int,timen::Int;schedulcolumn::Int=6)
  ## 数理モデル
  msub1 = JuMP.Model(Cbc.Optimizer)

   ## 変数
  JuMP.@variable(msub1, s[1:jobn,1:timen] ,Bin ) #ジョブ毎のコマ割りごとに実行有無を表示

   ## 定数
   runtime=zeros(Int,jobn,1) #ジョブ毎の実行時間

   ## ジョブ実行時間
   ###割り当てるコマ数を算出
   for i in 1:jobn
     runtime[i]=schedulesubplan[i,:runtime]/run_unit_time
   end

   ## 目的関数
   JuMP.@objective(msub1, Min, sum(schedulesubplan[1:end,:runtime])/run_unit_time-sum(s[1:jobn,1:timen] ))

   ## 制約条件
   ### ジョブ実行時間予約指定
   for i in 1:jobn
     flag =true
     if( convert(Bool,schedulesubplan[i,:Specifiedtime]) )
       for j in schedulcolumn:schedulcolumn+timen-1
         if(schedulesubplan[i,j] != 0 && flag)
           index=j-schedulcolumn+1
           JuMP.@constraint(msub1, s[i,index:index+runtime[i]-1] .== 1.0 )
            flag=false
         end
       end
     end
   end

   ### ジョブ実行時間制限
   for i in 1 :jobn
     JuMP.@constraint(msub1,sum(s[i,1:timen])==Float64(runtime[i]))
   end

   ### ロボット同時実行制限(ロボット数を超過させない)
   for i in 1 :timen
     JuMP.@constraint(msub1,sum(s[1:jobn,i])<=Float64(robotn) )
   end

   ## ソルバーの実行
   status = JuMP.optimize!(msub1)

   ## スケジュール案表示
   subplan=zeros(Int,jobn,timen)
   try
    subplan=map(Int,map(round,JuMP.value.(s)))
   catch
    subplan=zeros(Int,jobn,timen)
   end

 return subplan
end


text:
workspace()
using JukaiDepParser

p = Parser()
JukaiDepParser.train(p)


text:
# Importing modules inside functions is not allowed in Julia and must be done at top-level
using reactor: max_temperature

"""
    check_reactor_temperature(temperature_celsius)

Checks whether temperature is above max_temperature
and returns a status.
"""
function check_reactor_temperature(temperature_celsius)
    if temperature_celsius > max_temperature
        status = 1
    else
        status = 0
    end
end


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

mutable struct DedicatedHostAvailableCapacity <: SwaggerModel
    allocatableVMs::Any # spec type: Union{ Nothing, Vector{DedicatedHostAllocatableVM} } # spec name: allocatableVMs

    function DedicatedHostAvailableCapacity(;allocatableVMs=nothing)
        o = new()
        validate_property(DedicatedHostAvailableCapacity, Symbol("allocatableVMs"), allocatableVMs)
        setfield!(o, Symbol("allocatableVMs"), allocatableVMs)
        o
    end
end # type DedicatedHostAvailableCapacity

const _property_map_DedicatedHostAvailableCapacity = Dict{Symbol,Symbol}(Symbol("allocatableVMs")=>Symbol("allocatableVMs"))
const _property_types_DedicatedHostAvailableCapacity = Dict{Symbol,String}(Symbol("allocatableVMs")=>"Vector{DedicatedHostAllocatableVM}")
Base.propertynames(::Type{ DedicatedHostAvailableCapacity }) = collect(keys(_property_map_DedicatedHostAvailableCapacity))
Swagger.property_type(::Type{ DedicatedHostAvailableCapacity }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DedicatedHostAvailableCapacity[name]))}
Swagger.field_name(::Type{ DedicatedHostAvailableCapacity }, property_name::Symbol) =  _property_map_DedicatedHostAvailableCapacity[property_name]

function check_required(o::DedicatedHostAvailableCapacity)
    true
end

function validate_property(::Type{ DedicatedHostAvailableCapacity }, name::Symbol, val)
end


text:
module ThreadedIterables

include("ensureThreaded.jl")

include("foreach.jl");
include("map.jl");
include("mapreduce.jl");
include("mapfoldl.jl");
include("mapfoldr.jl");
include("reduce.jl");

include("threadedMacro.jl");

end


text:
include("type_FoamPreprocess.jl")
using HDF5, JLD
type FOAM_Data <: Data_Type_Selector
  prefix :: ASCIIString
  var :: ASCIIString
end

OceanT_FOAM = FOAM_Data("/ocean/history.ocean.", "T");

function makereading(dim,p :: Foam_Data_Path, day :: Int64,
                     Data_by :: FOAM_Data = OceanT_FOAM,
                     layer_selector :: Layer_1st = SurfaceLayer)
  data_path = string(p.dir_input,Data_by.prefix, @sprintf("%07.7i",day),".nc");
  (start, count) = select_layer(dim, layer_selector);
  res = ncread(data_path, Data_by.var, start= start, count = count);
  ncclose(data_path);
  return(res);
end

function makereading_all(p :: Foam_Data_Path,
                         Data_by :: FOAM_Data = OceanT_FOAM,
                         layer_selector :: Layer_1st = SurfaceLayer,
                         months_selector :: Annual_Months = Months_each_Year
                         ; check_starter :: Bool = true, set_NaN :: Bool = true)
  first_day = p.start_on_month
  last_day = p.last_on_month - 360
  if (check_starter)
    first_day = (first_day+330); first_day = first_day - (first_day%360);
  end
  years_days = [first_day:360:last_day]

  data_path = string(p.dir_input, Data_by.prefix, @sprintf("%07.7i",first_day),".nc");
  dim = size(ncread(data_path,Data_by.var));
  sst = fill(NaN, (dim[1], dim[2],1, size(months_selector.seasonidx,2),length(years_days)) )
  sst_temp = fill(NaN, (dim[1],dim[2],1,12))
  for (index, value) in enumerate(years_days)
    for (imon,vmon) in enumerate(0:30:330)
      sst_temp[:,:,:,imon] = makereading(dim,p,vmon+value , Data_by,layer_selector)
    end
    sst[:,:,:,:,index] = select_months(sst_temp, months_selector);
  end
  if (set_NaN)
    sst[sst.> 1000] = NaN
    sst[sst.< -1000] = NaN
  end
  return (sst)
end
function make_rw(p :: Foam_Data_Path,
                         Data_by :: FOAM_Data = OceanT_FOAM,
                         layer_selector :: Layer_1st = SurfaceLayer,
                         months_selector :: Annual_Months = Months_each_Year
                         ; check_starter :: Bool = true, set_NaN :: Bool = true)
  first_day = p.start_on_month
  last_day = p.last_on_month - 360
  if (check_starter)
    first_day = (first_day+330); first_day = first_day - (first_day%360);
  end
  years_days = [first_day:360:last_day]
  sst = makereading_all(p,Data_by, layer_selector,
                        months_selector,check_starter = check_starter, set_NaN= set_NaN);
  sst = float64(squeeze(sst,3))

  data_path = string(p.dir_input, Data_by.prefix, @sprintf("%07.7i",first_day),".nc");
  lon = float64(ncread(data_path,"lon"));
  lat = float64(ncread(data_path,"lat"));
  outputfile = string(p.filename,".jld");
  if (isfile(outputfile))
    rm(outputfile)
  end
  save(outputfile,"lon", lon, "lat",lat,"time",years_days,"sst", sst);
end

#=   latdim = NetCDF.NcDim("lat",lat,lat_att) =#
#=   londim = NetCDF.NcDim("lon",lon,lon_att) =#
#=   timedim = NetCDF.NcDim("times", years_days,@Compat.Dict("Units"=>"Day")) =#
#=   permonth = size(months_selector.seasonidx,1) =#
#=   monthdim = NetCDF.NcDim("month", [1:size(months_selector.seasonidx,2)], =#
#=                                     Dict{ASCIIString,ASCIIString}("Units"=>"$permonth months")) =#
#=   sst_att = Dict{ASCIIString,ASCIIString}("Units"=>"C"); =#
#=   sstVar = NetCDF.NcVar("sst", [londim,latdim,monthdim,years_days],sst_att, Float64) =#
#=   sstAVar = NetCDF.NcVar("sst_annual", [londim,latdim,years_days], sst_att, Float64) =#

#=   nc = NetCDF.create(p.filename, [sstVar,sstVar]); =#
#=   NetCDF.putVar(nc,"sst", sst); =#
#=   NetCDF.putVar(nc,"sst_annual", squeeze(mean(sst,3),1)); =#
#=   NetCDF.close(nc); =#


text:
__precompile__()

"""
This module simulates gene flow by gene dropping.
"""
module MendelGeneDropping
#
# Required OpenMendel packages and modules.
#
using MendelBase
# namely: DataStructures, GeneralUtilities
#
# Required external modules.
#
using CSV
using DataFrames

export GeneDropping

"""
This is the wrapper function for the Gene Dropping analysis option.
"""
function GeneDropping(control_file = ""; args...)
  #
  # Print the logo. Store the initial directory.
  #
  print(" \n \n")
  println("     Welcome to OpenMendel's")
  println("  Gene Dropping analysis option")
  print(" \n \n")
  println("Reading the data.\n")
  initial_directory = pwd()
  #
  # The user specifies the analysis to perform via a set of keywords.
  # Start the keywords at their default values.
  #
  keyword = set_keyword_defaults!(Dict{AbstractString, Any}())
  #
  # Keywords unique to this analysis should be first defined here
  # by setting their default values using the format:
  # keyword["some_keyword_name"] = default_value
  #
    # Allowed values are: "Unordered", "Ordered", "Sourced", or "Population".
  keyword["gene_drop_output"] = "Unordered"
  keyword["interleaved"] = true
  keyword["keep_founder_genotypes"] = false
  keyword["missing_data_pattern"] = "ExistingData" # Not yet implemented.
  keyword["missing_rate"] = 0.0
  keyword["repetitions"] = 1
  #
  # Process the run-time user-specified keywords that will control the analysis.
  # This will also initialize the random number generator.
  #
  process_keywords!(keyword, control_file, args)
  #
  # Check that the correct analysis option was specified.
  #
  lc_analysis_option = lowercase(keyword["analysis_option"])
  if (lc_analysis_option != "" &&
      lc_analysis_option != "genedropping")
     throw(ArgumentError(
       "An incorrect analysis option was specified.\n \n"))
  end
  keyword["analysis_option"] = "GeneDropping"
  #
  # Read the genetic data from the external files named in the keywords.
  #
  (pedigree, person, nuclear_family, locus, snpdata,
    locus_frame, phenotype_frame, person_frame, snp_definition_frame) =
    read_external_data_files(keyword)
  #
  # Check if SNP data were read.
  #
  if snpdata.snps != 0
    println(" \n\nERROR: This analysis does not use data from SNP files!\n")
  else
  #
  # Execute the specified analysis.
  #
    println(" \nAnalyzing the data.\n")
    execution_error = false
    new_person_frame = genedropping_option(pedigree, person,
      locus, locus_frame, phenotype_frame, person_frame, keyword)
    show(new_person_frame)
    if execution_error
      println(" \n \nERROR: Mendel terminated prematurely!\n")
    else
      println(" \n \nMendel's analysis is finished.\n")
    end
  end
  #
  # Finish up by closing, and thus flushing, any output files.
  # Return to the initial directory.
  #
  close(keyword["output_unit"])
  cd(initial_directory)
  return nothing
end # function GeneDropping

"""
This function conducts gene dropping at the model loci.
The results are placed in a dataframe.
The process is controlled by the keywords gene_drop_output
(with values Unordered, Ordered, Sourced, and Population);
missing_rate; new_pedigree_file; repetitions; seed; and
and the boolean keyword 'interleaved'.
"""
function genedropping_option(pedigree::Pedigree, person::Person,
  locus::Locus, locus_frame::DataFrame,
  phenotype_frame::DataFrame, person_frame::DataFrame,
  keyword::Dict{AbstractString, Any})

  (rows, columns) = size(person_frame)
  new_person_frame = deepcopy(person_frame[1:1, 1:columns])
  allowmissing!(new_person_frame)
  gene_drop_output = keyword["gene_drop_output"]
  #
  # Determine whether to use the ordered or unordered allele separator in output.
  # In either case use the first character in the defined strings.
  #
  if gene_drop_output == "Unordered"
    separator_char = keyword["allele_separator"][1]
  elseif gene_drop_output == "Ordered"
    separator_char = keyword["ordered_allele_separator"][1]
  elseif gene_drop_output == "Sourced"
    separator_char = keyword["ordered_allele_separator"][1]
  elseif gene_drop_output == "Population"
    separator_char = keyword["ordered_allele_separator"][1]
  else
    throw(ArgumentError(
      "Illegal value for gene_drop_output. OpenMendel terminated!\n \n"))
  end
  separator = ascii(string(separator_char))
  #
  # Loop over all pedigree/repeat combinations.
  #
  if keyword["interleaved"]
    for repetition = 1:keyword["repetitions"]
      for ped = 1:pedigree.pedigrees
        #
        # Copy the relevant part of the pedigree data frame.
        #
        pedstart = pedigree.pedstart[ped]
        twinfinish = pedigree.twin_finish[ped]
        extent = twinfinish - pedstart + 1
        frame = deepcopy(person_frame[pedstart:twinfinish, 1:columns])
        allowmissing!(frame)
        #
        # Perform gene dropping.
        #
        (sampled_genotype, source) =
          simulate_genotypes(pedigree, person, locus, keyword, ped)
        #
        # Name the pedigree replicate.
        #
        frame[1:extent, :Pedigree] .= pedigree.name[ped] * string(repetition)
        #
        # Output can be sampled alleles or sources.
        #
        if gene_drop_output == "Unordered" || gene_drop_output == "Ordered"
          converted_genotype = convert_sampled_genotype(locus,
            sampled_genotype, separator, gene_drop_output)
        else
          converted_genotype = convert_sampled_genotype(locus,
            source, separator, gene_drop_output)
        end
        #
        # Load the simulated data.
        #
        for loc = 1:locus.loci
          l = locus.locus_field_in_person_frame[loc]
          frame[1:extent, l] = converted_genotype[:, loc]
        end
        #
        # Sort the new pedigree so that it is the same order as the input file.
        # NB: Thus the new pedigree frame will not agree
        # with the order of the individuals in the data structures.
        # Append the simulated data to the new pedigree frame.
        #
        sort!(frame, :EntryOrder)
        append!(new_person_frame, frame)
      end
    end
  else
    for ped = 1:pedigree.pedigrees
      #
      # Copy the relevant part of the pedigree data frame.
      #
      pedstart = pedigree.pedstart[ped]
      pedfinish = pedigree.pedfinish[ped]
      extent = pedfinish - pedstart + 1
      #
      # Perform gene dropping.
      #
      for repetition = 1:keyword["repetitions"]
        frame = deepcopy(person_frame[pedstart:pedfinish, 1:columns])
        allowmissing!(frame)
        (sampled_genotype, source) =
          simulate_genotypes(pedigree, person, locus, keyword, ped)
        #
        # Name the pedigree replicate.
        #
        frame[1:extent, :Pedigree] = pedigree.name[ped] * string(repetition)
        #
        # Output can be sampled alleles or sources.
        #
        if gene_drop_output == "Unordered" || gene_drop_output == "Ordered"
          converted_genotype = convert_sampled_genotype(locus,
            sampled_genotype, separator, gene_drop_output)
        else
          converted_genotype = convert_sampled_genotype(locus,
            source, separator, gene_drop_output)
        end
        #
        # Load the simulated data.
        #
        for l = 1:locus.model_loci
          l = locus.locus_field_in_person_frame[loc]
          frame[1:extent, l] = converted_genotype[:, loc]
        end
        #
        # Sort the new pedigree so that it is the same order as the input file.
        # NB: Thus the new pedigree frame will not agree
        # with the order of the individuals in the data structures.
        # Append the simulated data to the new pedigree frame.
        #
        sort!(frame, :EntryOrder)
        append!(new_person_frame, frame)
      end
    end
  end
  #
  # If an output file was named for the simulations, then output the data.
  # First, remove the first row, which is a stub of the input pedigree.
  # Second, remove the EntryOrder or Inverse_Perm columns.
  #
##  deleterows!(new_person_frame, 1)
  new_person_frame = new_person_frame[2:end, :]
  new_pedigree_file = string(keyword["new_pedigree_file"])
  if new_pedigree_file != ""
    names_list = names(new_person_frame)
    deleteat!(names_list, findall((in)([:EntryOrder]), names_list))
    deleteat!(names_list, findall((in)([:Inverse_Perm]), names_list))
    CSV.write(new_pedigree_file, new_person_frame[:, names_list];
      writeheader = true, delim = keyword["output_field_separator"],
      missingstring = keyword["output_missing_value"])
  end
  return new_person_frame
end # function gene_dropping_option

"""
This function carries out gene dropping for a single pedigree
over the model loci. These are the loci common to the Pedigree
and Locus frames.
"""
function simulate_genotypes(pedigree::Pedigree, person::Person,
  locus::Locus, keyword::Dict{AbstractString, Any}, ped::Int)

  theta = locus.theta
  #
  # Assign a source label to each founder gene.
  #
  source =
    founder_source(pedigree, person, locus, ped, keyword["gene_drop_output"])
  #
  # Compute pedigree attributes and offsets.
  #
  populations = person.populations
  founders = pedigree.founders[ped]
  ped_size = pedigree.individuals[ped]
  descendants = ped_size - founders
  model_loci = locus.model_loci
  q = pedigree.pedstart[ped] - 1 # offset for founders
  r = q + founders # offset for descendants
  #
  # Allocate arrays.
  #
  sampled_genotype = zeros(Int, 2, ped_size, model_loci)
  gene = ones(Int, 2, ped_size)
  meiosis = ones(Int, 2, descendants)
  old_meiosis = ones(Int, 2, descendants)
  #
  # Loop over all model loci.
  #
  for l = 1:model_loci
    loc = locus.model_locus[l]
    xlinked = locus.xlinked[loc]
    #
    # Choose a compatible ordered genotype for each founder.
    # Each allele frequency is a convex combination of the
    # allele frequencies in the component populations.
    # The frequency variable is a column vector.
    #
    for i = 1:founders
      frequency =
        vec(transpose(transpose(person.admixture[i + q, :]) *
	              locus.frequency[loc]))
      male = person.male[i + q]
      #
      # Generate a founder genotype at the current locus.
      #
      if keyword["keep_founder_genotypes"]
        genotype_set = person.genotype[i + q, l]
        gene[:, i] = choose_genotype(frequency, genotype_set, xlinked, male)
      else
        gene[:, i] = random_genotype(frequency, xlinked, male)
      end
    end
    #
    # Propagate these choices to the descendants.
    #
    for i = 1:descendants
      for j = 1:2
        if j == 1
          parent = person.mother[i + r] - q
        else
          parent = person.father[i + r] - q
        end
        #
        # Determine the parental source of each gamete and store
        # a value of 1 or 2 in the meiosis indicator.
        #
        uniform = rand()
        if l == 1
          meiosis[j, i] = round(Integer, uniform) + 1
        else
          #
          # Take into account recombination after the first locus.
          #
          if uniform < theta[j, l - 1]
            meiosis[j, i] = 3 - old_meiosis[j, i]
          else
            meiosis[j, i] = old_meiosis[j, i]
          end
        end
        #
        # Record the allele passed and the source.
        #
        gene[j, i + founders] = gene[meiosis[j, i], parent]
        source[j, i + founders, l] = source[meiosis[j, i], parent, l]
      end
      if xlinked && person.male[i + r]
        source[2, i + founders, l] = source[1, i + founders, l]
      end
    end
    #
    # Save the sampled genotypes.
    #
    for j = 1:ped_size
      for i = 1:2
        sampled_genotype[i, j, l] = gene[i, j]
      end
    end
    #
    # Save the meiosis indicators.
    #
    if l < locus.model_loci
      old_meiosis = deepcopy(meiosis)
    end
  end
  #
  # Extend the simulation to co-twins.
  #
  for i = pedigree.pedfinish[ped] + 1:pedigree.twin_finish[ped]
    m = person.primary_twin[i]
    sampled_genotype[:, i - q, :] = sampled_genotype[:, m - q, :]
    source[:, i - q, :] = source[:, m - q, :]
  end
  #
  # Delete genotypes that should be missing.
  #
  for l = 1:locus.model_loci
    for i = 1:ped_size
      if rand() < keyword["missing_rate"]
        sampled_genotype[:, i, l] = [0 0]
      end
    end
  end
  return (sampled_genotype, source)
end # function simulate_genotypes

"""
This function computes source labels for founder genes in a single pedigree.
Gene sources can be labeled by founders or ancestral populations.
"""
function founder_source(pedigree::Pedigree, person::Person,
  locus::Locus, ped::Int, gene_drop_output::AbstractString)

  q = pedigree.pedstart[ped] - 1
  founders = pedigree.founders[ped]
  ped_size = pedigree.individuals[ped]
  model_loci = locus.model_loci
  source = zeros(Int, 2, ped_size, model_loci)
  #
  # Assign a source number to each founder gene.
  #
  for i = 1:founders
    if gene_drop_output == "Population"
      for l = 1:model_loci
        source[1, i, l] = random_category(vec(person.admixture[i + q, :]))
        source[2, i, l] = random_category(vec(person.admixture[i + q, :]))
      end
    else
      source[1, i, :] .= 2i - 1
      source[2, i, :] .= 2i
    end
  end
  #
  # At an x-linked locus, male founders have only one source.
  #
  for l = 1:model_loci
    loc = locus.model_locus[l]
    if locus.xlinked[loc]
      for i = 1:founders
        if person.male[i + q]
          source[2, i, l] = source[1, i, l]
        end
      end
    end
  end
  return source
end # function founder_source

"""
This function randomly chooses an ordered genotype from
a set of ordered genotypes. If the set is empty,
it returns a completely random genotype.
"""
function choose_genotype(frequency::Vector{Float64},
  genotype_set::Set{Tuple{Int, Int}}, xlinked::Bool, male::Bool)

  if length(genotype_set) == 0
    return random_genotype(frequency, xlinked, male)
  end
  #
  # Normalize the probabilities of the possible genotypes.
  #
  prob = zeros(length(genotype_set))
  i = 0
  for (gm, gp) in genotype_set
    i = i + 1
    if xlinked && male
      if gm == gp
        prob[i] = frequency[gm]
      end
    else
      prob[i] = frequency[gm] * frequency[gp]
    end
  end
  #
  # Check whether any compatible genotypes exist.
  #
  total = sum(prob)
  if total <= 0.0
    return random_genotype(frequency, xlinked, male)
  end
  #
  # Normalize the probabilities and return a sampled genotype.
  #
  prob = prob / total
  j = random_category(prob)
  i = 0
  for (gm, gp) in genotype_set
    i = i + 1
    if i == j
      return [gm gp]
    end
  end
end # function choose_genotype

"""
This function samples a random ordered genotype from
the universe of possible ordered genotypes.
"""
function random_genotype(frequency::Vector{Float64}, xlinked::Bool, male::Bool)

  if xlinked && male
    i = random_category(frequency)
    return [i i]
  else
    i = random_category(frequency)
    j = random_category(frequency)
    return [i j]
  end
end # function random_genotype

"""
This function converts sampled numerical genotypes into ordinary genotypes.
"""
function convert_sampled_genotype(locus::Locus, sampled_genotype::Array{Int, 3},
  separator::AbstractString, gene_drop_output::AbstractString)

  (m, n) = (size(sampled_genotype, 2), size(sampled_genotype, 3))
  converted_genotype = Array{AbstractString}(undef, m, n)
  if gene_drop_output == "Unordered" || gene_drop_output == "Ordered"
    #
    # Convert allele numbers to allele names and concatenate.
    #
    for l = 1:locus.model_loci
      loc = locus.model_locus[l]
      for j = 1:m
        a = locus.allele_name[loc][sampled_genotype[1, j, l]]
        b = locus.allele_name[loc][sampled_genotype[2, j, l]]
        converted_genotype[j, l] = a * separator * b
      end
    end
  else
    #
    # Convert alleles numbers to integer strings and concatenate.
    #
    for l = 1:locus.model_loci
      loc = locus.model_locus[l]
      for j = 1:m
        a = string(sampled_genotype[1, j, l])
        b = string(sampled_genotype[2, j, l])
        converted_genotype[j, l] = a * separator * b
      end
    end
  end
  return converted_genotype
end # function convert_sampled_genotype
#
# Method to obtain path to this package's data files
# so they can be used in the documentation and testing routines.
# For example, datadir("Control file.txt") will return
# "/path/to/package/data/Control file.txt"
#
datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)

end # module MendelGeneDropping


text:
# cartpole_train_sample.jl
using Gym
using Knet
using ArgParse

include(joinpath(@__DIR__, "dqn", "dqn.jl"))

main() = main(ARGS)
main(args::String) = main(split(args))
function main(args::Vector{<:AbstractString})
    s = ArgParseSettings()
    # s.description = "(c) Ozan Arkan Can, 2018. An implementation of the deep q network."
    @add_arg_table s begin
        ("--frames"; arg_type=Int; default=30001; help="number of frames")
        ("--lr"; arg_type=Float64; default=0.001; help="learning rate")
        ("--gamma"; arg_type=Float64; default=0.99; help="discount factor")
        ("--hiddens"; arg_type=Int; nargs='+'; default=[32]; help="number of units in the hiddens for the mlp")
        ("--env_id"; default="CartPole-v0")
        ("--render"; action=:store_true)
        ("--memory"; arg_type=Int; default=10000; help="memory size")
        ("--bs"; arg_type=Int; default=32; help="batch size")
        ("--stack"; arg_type=Int; default=4; help="length of the frame history")
        # ("--save"; help="model name")
        ("--load"; help="load model name")
        ("--save_interval"; arg_type=Int; default=10000; help="save interval (per frames)")
        ("--no_save"; action=:store_true; help="to prevent to save the model")
        ("--model_prefix"; arg_type=String; help="model name prefix (default: YYYYmmddHHMMSS)")
        ("--play"; action=:store_true; help="only play")
        ("--printinfo"; action=:store_true; help="print the training messages")
        ("--tupdate"; arg_type=Int; default=500; help="update frequency for the target network")
    end
    if "--help" ∈ args || "-h" ∈ args
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    o = parse_args(args, s)
    o["atype"] = Array{Float32}
    # o["play"] = false

    if o["model_prefix"] === nothing
        o["model_prefix"] = Dates.format(Dates.now(), Dates.dateformat"YYYYmmddHHMMSS")
    end

    # train
    env = GymEnv(o["env_id"])
    # seed!(env, 12345)
    # cartpoleenv = env.gymenv
    # cartpoleenv[:_max_episode_steps] = typemax(Int)
    # cartpoleenv[:_max_episode_seconds] = typemax(Int)

    if o["load"] === nothing
        INPUT = env.observation_space.shape[1] * o["stack"]
        OUTPUT = env.action_space.n
        w = DQN.init_weights(INPUT, o["hiddens"], OUTPUT, o["atype"])
        o["n_hiddens"] = length(o["hiddens"])
    else
        w = DQN.load_model(o["load"], o["atype"])
        o["n_hiddens"] = length(w) ÷ 2 - 1
    end

    opts = Dict(k => Rmsprop(lr=o["lr"]) for k in keys(w))
    buffer = DQN.ReplayBuffer(o["memory"])
    exploration = DQN.PiecewiseSchedule([(0, 1.0),
        (round(Int, o["frames"]/5), 0.1)])
    # exploration = DQN.PiecewiseSchedule([(0, 1.0), (min(o["frames"] ÷ 3, 10000), 0.1)])

    rewards, frames = DQN.dqn_learn(w, opts, env, buffer, exploration, o)

    close!(env)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

max_stars_repo_name: JuliaTagBot/JuliaGymDemo

max_issues_repo_name: JuliaTagBot/JuliaGymDemo

max_forks_repo_name: JuliaTagBot/JuliaGymDemo


text:
function test_sigmoid_neuron(backend::Backend, T, eps)
  println("-- Testing Sigmoid neuron on $(typeof(backend)){$T}...")

  data = rand(T, 3,4,5,6) - 0.5
  data_blob = make_blob(backend, data)
  neuron = Neurons.Sigmoid()

  println("    > Forward")
  forward(backend, neuron, data_blob)
  expected_data = 1 ./ (1 + exp(-data))
  got_data = zeros(T, size(data))
  copy!(got_data, data_blob)

  @test all(-eps .< got_data - expected_data .< eps)

  println("    > Backward")
  grad = rand(T, size(data))
  grad_blob = make_blob(backend, grad)
  backward(backend, neuron, data_blob, grad_blob)

  expected_grad = grad .* (expected_data .* (1-expected_data))
  got_grad = zeros(T, size(expected_grad))
  copy!(got_grad, grad_blob)

  @test all(-eps .< got_grad - expected_grad .< eps)
end
function test_sigmoid_neuron(backend::Backend)
  test_sigmoid_neuron(backend, Float32, 1e-3)
  test_sigmoid_neuron(backend, Float64, 1e-9)
end

if test_cpu
  test_sigmoid_neuron(backend_cpu)
end
if test_gpu
  test_sigmoid_neuron(backend_gpu)
end


text:
@static if VERSION < v"1.1"
    mutable struct Event
        lock::Base.Threads.Mutex
        q::Vector{Task}
        set::Bool
        # TODO: use a Condition with its paired lock
        Event() = new(Base.Threads.Mutex(), Task[], false)
    end

    function Base.wait(e::Event)
        e.set && return
        lock(e.lock)
        while !e.set
            ct = current_task()
            push!(e.q, ct)
            unlock(e.lock)
            try
                wait()
            catch
                filter!(x->x!==ct, e.q)
                rethrow()
            end
            lock(e.lock)
        end
        unlock(e.lock)
        return nothing
    end

    function Base.notify(e::Event)
        lock(e.lock)
        if !e.set
            e.set = true
            for t in e.q
                schedule(t)
            end
            empty!(e.q)
        end
        unlock(e.lock)
        return nothing
    end
elseif VERSION < v"1.2"
    using Base.Threads: Event
else
    using Base: Event
end


text:
# This file contains functions used in the high level service definitions

_merge(a::AbstractDict, b::AbstractDict) = merge(a, b)
_merge(a::AbstractString, b::AbstractString) = b
Base.merge(_merge, a::AbstractString, b::AbstractString) = _merge(a, b)


text:
function test_linop()
  @testset "linop" begin
    (nrow, ncol) = (10, 6)
    ϵ = eps(Float64)
    rtol = sqrt(ϵ)
    A1 = rand(nrow, ncol)
    opA1 = LinearOperator(A1)
    v = rand(ncol)
    res = zeros(nrow)

    # test mul3
    mul!(res, opA1, v) # build
    opA1 = LinearOperator(A1)
    allocs = @allocated mul!(res, opA1, v)
    @test allocs == 0
    res_true = zeros(nrow)
    mul!(res_true, A1, v)
    @test norm(res-res_true) ≤ sqrt(ϵ)

    # test mul5
    α, β = 2.0, -3.0
    res = rand(nrow)
    res_true .= res
    opA1 = LinearOperator(A1)
    mul!(res, opA1, v, α, β)
    mul!(res, A1, v, α, β)
    res .= res_true
    allocs2 = @allocated mul!(res, opA1, v, α, β)
    allocs = @allocated mul!(res_true, A1, v, α, β)
    @test allocs2 == 0
    @test allocs == 0
    @test norm(res-res_true) ≤ sqrt(ϵ)

    # test prod
    opA1 = LinearOperator(A1)
    res = opA1 * v
    res_true = A1 * v
    @test norm(res-res_true) ≤ sqrt(ϵ)

    # test +op
    opA2 = +opA1
    α, β = 2.0, -3.0
    res = rand(nrow)
    res_true .= res
    mul!(res, opA2, v, α, β)
    mul!(res_true, A1, v, α, β)
    @test norm(res-res_true) ≤ sqrt(ϵ)

    # test -op
    opA2 = -opA1
    α, β = 2.0, -3.0
    res = rand(nrow)
    res_true .= res
    opA1 = LinearOperator(A1)
    mul!(res, opA2, v, α, β)
    mul!(res_true, -A1, v, α, β)
    @test norm(res-res_true) ≤ sqrt(ϵ)
    mul!(res, opA2, v, α, β)

    # test LinearOperator with preallocated vectors
    popA1 = LinearOperator(A1)
    mul!(res, popA1, v, α, β) # compile
    popA1.Mv .= rand(nrow)
    res_true .= popA1.Mv
    mul!(popA1.Mv, popA1, v, α, β)
    mul!(res_true, A1, v, α, β)
    @test norm(popA1.Mv-res_true) ≤ sqrt(ϵ)

    (nrow2, ncol2) = (ncol, 7)
    A2 = rand(nrow2, ncol2)
    v2 = rand(ncol2)
    popA2 = LinearOperator(A2)
    popA3 = popA1 * popA2
    popA3.Mv .= rand(nrow)
    res_true .= popA3.Mv
    mul!(popA3.Mv, popA3, v2, α, β)
    mul!(res_true, A1*A2, v2, α, β)
    @test norm(popA3.Mv-res_true) ≤ sqrt(ϵ)

    # Matrix * LinearOperator
    opA3 = opA1 * A2
    opA3.Mv .= rand(nrow)
    res_true .= opA3.Mv
    mul!(opA3.Mv, opA3, v2, α, β)
    mul!(res_true, A1*A2, v2, α, β)
    @test norm(opA3.Mv-res_true) ≤ sqrt(ϵ)

    # test opEye
    # opE = opEye(Float64, nrow)
    # vec, res, res_true = ones(nrow), ones(nrow), ones(nrow)
    # mul!(res, opE, vec, α, β)
    # mul!(res_true, I, vec, α, β)
    # @test norm(res-res_true) ≤ sqrt(ϵ)
  end
end

test_linop()


text:
using MiniTB
using Test
using TBComponents
using CSV
using DataFrames
using StatsBase
using Utils
using BenchmarkTools

include("../src/web/web_model_libs.jl" )

print_test = false

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120
BenchmarkTools.DEFAULT_PARAMETERS.samples = 2

function basic_run( params, num_households, num_people, base_results)

    global print_test

    num_repeats = 100
    json_out=missing
    results=missing
    params.it_rate[1] = 0.66
    results = doonerun( params, num_households, num_people, num_repeats )
    summary_output = summarise_results!( results=results, base_results=base_results )

    if print_test
        print( "   deciles = $( summary_output.deciles)\n\n" )

        print( "   poverty_line = $(summary_output.poverty_line)\n\n" )

        print( "   inequality = $(summary_output.inequality)\n\n" )

        print( "   poverty = $(summary_output.poverty)\n\n" )

        print( "   gainlose_by_sex = $(summary_output.gainlose_by_sex)\n\n" )
        print( "   gainlose_by_thing = $(summary_output.gainlose_by_thing)\n\n" )

        print( "   metr_histogram= $(summary_output.metr_histogram)\n\n")
        println( "SUMMARY OUTPUT")
        println( summary_output )
        println( "as JSON")
        println( JSON.json( summary_output ))
    end

end # summary_output.timing summary_output.blockt = JSON.json( summary_output )

example_names, num_households, num_people = load_data( load_examples = true, load_main = true, start_year = 2015 )

params = deepcopy(DEFAULT_PARAMS)

base_results = create_base_results( num_households, num_people )

t = @benchmark basic_run( params, num_households, num_people, base_results ) seconds=120 samples=10 evals=1

print(t)


text:
using PyPlot, Revise, Random, ProgressMeter
#includet("2NNlogistic-param-predict.jl")
using3D()

##

"""
Gradient ascent to train a 2-layers ReLU neural net for the soft-min loss
INPUT: X (training input), Y (training output), m (nb neurons), both: training both layers or just the output
OUTPUT: Ws (training trajectory)
"""
function logistic_reg(X, Y, stepsize, niter)
    (n,d) = size(X)
    a     = zeros(d+1)
    as    = zeros(d+1, niter) # store optimization path
    betas = zeros(niter)
    margins = zeros(niter)
    Z = Y .* cat(X, 1*ones(n), dims=2)/10
    for iter = 1:niter
        as[:,iter] = a
        perf = Z*a
        margin = minimum(perf)
        p = exp.( margin .- Z*a)
        p = p./ sum(p)
        grad = Z'*p
        betas[iter] = maximum([1,sqrt(sum(a.^2))])
        margins[iter] = margin/sqrt(sum(a.^2))
        a = a + betas[iter] * stepsize * grad /(sqrt(iter+1))
    end
    as, betas, margins
end

##

"""
Gradient ascent to train a 2-layers ReLU neural net for the soft-min loss
INPUT: X (training input), Y (training output), m (nb neurons), both: training both layers or just the output
OUTPUT: Ws (training trajectory)
"""
function logistic_reg_lin(X, Y, stepsize, niter)
    (n,d) = size(X)
    a     = zeros(d)
    as    = zeros(d, niter) # store optimization path
    betas = zeros(niter)
    margins = zeros(niter)
    Z = Y .* X

    for iter = 1:niter
        as[:,iter] = a
        perf = Z*a
        margin = minimum(perf)
        p = exp.( margin .- Z*a)
        p = p./ sum(p)
        grad = Z'*p
        betas[iter] = maximum([1,sqrt(sum(a.^2))])
        margins[iter] = margin/sqrt(sum(a.^2))
        a = a + betas[iter] * stepsize * grad /(sqrt(iter+1))
    end
    as, betas, margins
end

# Random.seed!(4)
# nx = 3
# nxx = 50
# n = 2*(nx + nxx)
# XX = rand(nx,2)
# A1 = (XX[:,1] .+ max.(XX[:,2] .- XX[:,1], 0) .+0.0)/2.5
# A2 = XX[:,2] .- max.(XX[:,2] .- XX[:,1], 0) .- 0.5
# AA = [0; -0.5]' .+ 0.03*randn(nxx,2)
# XX = rand(nx,2)
# B1 = (XX[:,1] .+ min.(XX[:,2] .- XX[:,1], 0) .-1.0)/2.5
# B2 = XX[:,2] .- min.(XX[:,2] .- XX[:,1], 0) .-0.5
# BB = [-0.5; -0.5]' .+ 0.03*randn(nxx,2)
# X = [A1 A2;AA; B1 B2;BB]
# Y = cat(ones(nx+nxx), -ones(nx+nxx), dims=1)
# plot(X[Y.==1,1], X[Y.==1,2],"oC0",label=L"Y=+1")
# plot(X[Y.==-1,1], X[Y.==-1,2],"oC1",label=L"Y=-1")
# legend()

Random.seed!(4)
nx = 3
nxx = 50
n = 2*(nx + nxx)
XX = rand(nx,2)
A1 = (XX[:,1] .+ max.(XX[:,2] .- XX[:,1], 0) .+0.0)/2.5
A2 = XX[:,2] .- max.(XX[:,2] .- XX[:,1], 0) .- 0.5
AA = [0; -0.5]' .+ 0.03*randn(nxx,2)
XX = rand(nx,2)
B1 = (XX[:,1] .+ min.(XX[:,2] .- XX[:,1], 0) .-1.0)/2.5
B2 = XX[:,2] .- min.(XX[:,2] .- XX[:,1], 0) .-0.5
BB = [-0.5; -0.5]' .+ 0.03*randn(nxx,2)
X = [A1 A2;AA; B1 B2;BB]
Y = cat(ones(nx+nxx), -ones(nx+nxx), dims=1)
plot(X[Y.==1,1], X[Y.==1,2],"oC0",label=L"Y=+1")
plot(X[Y.==-1,1], X[Y.==-1,2],"oC1",label=L"Y=-1")
legend()

stepsize = 0.05
niter = 100000
as, betas, margins = logistic_reg(X, Y, stepsize, niter)
abars = as ./ betas'
plot(margins[end] .- margins)

# Create a sphere
fig = figure(figsize=[10,5])
ax1 = subplot(121, projection="3d")

res = 50
phi = range(0.0, 2*π, length=res)#.* ones(res)'

#theta =  ones(res) .* range(0.0, 2π, length=res)'
#x = sin.(phi) .* cos.(theta)
#y = sin.(phi) .* sin.(theta)
#z = cos.(phi)

plot3D(sin.(phi), cos.(phi),zeros(res),".k",ms=1)
plot3D(sin.(phi), zeros(res),cos.(phi),".k",ms=1)
plot3D(zeros(res),sin.(phi), cos.(phi),".k",ms=1)
#plot_surface(x, y, z,  rstride=1, cstride=1, color="C0", alpha=0.1, linewidth=0,shade=true)
ax1.set_xlim3d(-1.2, 1.2)
ax1.set_ylim3d(-1.2, 1.2)
ax1.set_zlim3d(-1.2, 1.2)
ax1.view_init(15,-80)
plot3D([0],[0],[0],"ok")
#plot3D(abars[1,:], abars[2,:],abars[3,:],"C0",lw=3)
N = max.(atan.(betas)./atan(1),1)
plot3D(abars[1,:] .* N, abars[2,:].* N,abars[3,:],"C1",lw=3)

ax2 = subplot(122)
plot(X[Y.==1,1], X[Y.==1,2],"oC0",label=L"Y=+1")
plot(X[Y.==-1,1], X[Y.==-1,2],"oC1",label=L"Y=-1")
ts = -2:0.1:2
iter=1000
    plot(-(abars[2,iter].*ts .+ abars[3,iter])/abars[1,iter], ts,"k")
iter=niter
    plot(-(abars[2,iter].*ts .+ abars[3,iter])/abars[1,iter], ts,"r")

legend()

stepsize = 0.5
niter = 500000
as, betas, margins = logistic_reg(X, Y, stepsize, niter)

nframes=200
resolution = 0.005
# define the sequence of time steps ts to be plotted
a = (niter-1)/(nframes-1)^4
ts = setdiff(Int.(floor.(a*(0:nframes-1).^4)) .+ 1,1)
as = as[:,ts]
N =tanh.(0.5*sqrt.(sum(as.^2, dims=1)))
alogs = as .* N./ sqrt.(sum(as.^2, dims=1))

@showprogress 1 "Plotting images..." for k = 1:length(ts)
       ioff() # turns off interactive plotting
    fig = figure(figsize=[7,4])
    ax1 = subplot(121, projection="3d")
    ax1.set_position([0,0.1,0.5,0.8])

    if k<21
        indt = 1:k
    else
        indt = 1:k #(k-20):k
    end

    plot3D(alogs[1,indt], alogs[2,indt],alogs[3,indt],"k",lw=1) # tail
    plot3D([alogs[1,k]], [alogs[2,k]],[alogs[3,k]],"o",color="C3", markersize=3)

    res = 50
    phi = range(0.0, 2π, length=res)#.* ones(res)'
    r = 0.5
    plot3D(r*sin.(phi), r*cos.(phi),zeros(res),".k",ms=0.3)
    plot3D(r*sin.(phi), zeros(res),r*cos.(phi),".k",ms=0.3)
    plot3D(zeros(res),r*sin.(phi), r*cos.(phi),".k",ms=0.3)

    ax1.set_xlim3d(-1, 1)
    ax1.set_ylim3d(-1, 1)
    ax1.set_zlim3d(-1, 1)
    ax1.set_xticks([-1/2, 0, 1/2])
    ax1.set_yticks([-1/2, 0, 1/2])
    ax1.set_zticks([-1/2, 0, 1/2])
    ax1.view_init(25-20*sinpi(k/length(ts)),-45*cospi(k/length(ts)))
   # ax1.view_init(25-20*sinpi(0.2),45-45*cospi(0.2))

    ax2 = subplot(122)
    ax2.set_position([0.45,0.25,0.5,0.5])

    f(x1,x2,k) = sum( as[:,k] .* [x1;x2;1.0]) # prediction function

    xs = -0.8:resolution:0.8
    tab = [f(xs[i],xs[j],k) for i=1:length(xs), j=1:length(xs)]
    pcolormesh(xs', xs, tanh.(tab'), cmap="coolwarm", shading="gouraud", vmin=-1.0, vmax=1.0, edgecolor="face")

    xs = -0.8:resolution:0.8
    tab = [f(xs[i],xs[j],k) for i=1:length(xs), j=1:length(xs)]
    contour(xs', xs, tanh.(tab'), levels =0, colors="k", antialiased = true, linewidths=2)

    # plot training set
    X1 = X[(Y .== 1),:]
    X2 = X[(Y .== -1),:]
    plot(X1[:,1],X1[:,2],"+k")
    plot(X2[:,1],X2[:,2],"_k")
    axis("equal");axis("off");
    ax2.set_xticks([-1/2, 0, 1/2])
    ax2.set_yticks([-1/2, 0, 1/2])

    savefig("dynamics_$(k).png",bbox_inches="tight", dpi=300)
    close(fig)
end

ion()

Random.seed!(7)
X = rand(4,2)
X = [-3 -1; -1 1; 1 -2; 2 1]
Y = cat(-ones(2), ones(2), dims=1)
as, betas, margins = logistic_reg_lin(X, Y, 0.1, 1000)
abars = as ./ betas'
a = abars[:,end]

fig = figure(figsize=[4,4*6/7])
plot(X[Y.==1,1], X[Y.==1,2],"oC3",label=L"y=+1")
plot(X[Y.==-1,1], X[Y.==-1,2],"oC0",label=L"y=-1")
ts = -5:0.01:5
fill_betweenx(ts, -a[2].*ts /a[1],4,alpha=0.1,color="C3")
fill_betweenx(ts, -4,-a[2].*ts /a[1],alpha=0.1,color="C0")
arrow(0,0,0.9*a[1],0.9*a[2],lw=2,head_width=0.15,color="k",overhang=0)
text(0.6,-0.2,L"a")
plot([0],[0],"ok")
#text(-0.95,0.1,L"(0,0)")
ts=0.0:0.01:2.0
plot(1.2*cospi.(ts),1.2*sinpi.(ts),":k",label="unit circle")
f = 2.2
plot([X[1,1],X[1,1]+f*a[1]],[X[1,2], X[1,2] + f*a[2]],"--C0")
f = 1.3
plot([X[2,1],X[2,1]+f*a[1]],[X[2,2], X[2,2] + f*a[2]],"--C0")
f = -1.8
plot([X[3,1],X[3,1]+f*a[1]],[X[3,2], X[3,2] + f*a[2]],"--C3")
f = -1.4
plot([X[4,1],X[4,1]+f*a[1]],[X[4,2], X[4,2] + f*a[2]],"--C3")
text(X[1,1]+0.2,X[1,2],L"x_1")
text(X[2,1]-0.5,X[2,2]-0.1,L"x_2")
text(X[3,1]+0.2,X[3,2],L"x_3")
text(X[4,1]+0.2,X[4,2],L"x_4")
axis([-3.5,3.5,-2.5,3.5])
arrow(1.2,1.7,0.8*a[1],0.8*a[2],lw=1,head_width=0.15,color="k",overhang=0)
f =0.2
arrow(1.2+f*a[1],1.7+f*a[2],-f*a[1],-f*a[2],lw=1,head_width=0.15,color="k",overhang=0)
text(1.4,1.8,L"margin of $a$")
ts = -5:0.01:5
plot(-a[2].*ts /a[1], ts,"k",lw=3)
axis("off")
legend(loc=2)
savefig("max_margin.png",dpi=200, bbox_inches="tight")

?arrow

SS(u,v) = - log(exp(log(1+exp(-u))/2 + log(1+exp(-v))/2) - 1)
S(u,v) = - log(exp(-u)/2 +exp(-v)/2)
u = -20:0.05:20
v = -20:0.05:20

plot(u,SS.(u,10))

plot(u[1:end-1],SS.(u,10)[2:end].- SS.(u,10)[1:end-1])

surf(u,v,SS.(u,v'))

show()

max_stars_repo_name: brechetp/ImplicitBias

max_issues_repo_name: brechetp/ImplicitBias

max_forks_repo_name: brechetp/ImplicitBias


text:
module Performance

using StandardizedMatrices, BenchmarkTools, StatsBase

n, p = 1000, 1000
x = randn(n, p)
x2 = zscore(x, 1)
z = StandardizedMatrix(x)
β = randn(p)
isapprox(x2 * β, z * β) ?
	info("Dense Matrix-Vector multiplication correct") :
	warn("Dense Matrix-Vector multiplication incorrect")
x = sprandn(n, p, .01)
x2 = zscore(x, 1)
z = StandardizedMatrix(x)
isapprox(x2 * β, z * β) ?
	info("Sparse Matrix-Vector multiplication correct") :
	warn("Sparse Matrix-Vector multiplication incorrect")

info("Dense A_mul_B! timing ratio")
b1 = @benchmark(
	A_mul_B!(storage, x, β),
	setup = (storage = zeros(n); x = StandardizedMatrix(randn(n, p)); β = randn(p))
)
b2 = @benchmark(
	A_mul_B!(storage, x, β),
	setup = (storage = zeros(n); x = randn(n, p); β = randn(p))
)
@show ratio(minimum(b1), minimum(b2))

info("Sparse A_mul_B! timing ratio")
b1 = @benchmark(
	A_mul_B!(storage, x, β),
	setup = (storage = zeros(n); x = StandardizedMatrix(sprandn(n, p, .05)); β = randn(p))
)
b2 = @benchmark(
	A_mul_B!(storage, x, β),
	setup = (storage = zeros(n); x = zscore(sprandn(n, p, .05), 1); β = randn(p))
)
@show ratio(minimum(b1), minimum(b2))

end


text:
function readinput(filename)
    map(group -> filter(!isempty, split(group, "\n")), split(read(filename, String), "\n\n"))
end

P1, P2 = 0, 0

groups = readinput("d06.in")

for group in groups
    global P1 += length(Set(join(group)))
    commonanswers = Set(first(group))
    for answers in group
        intersect!(commonanswers, answers)
    end
    global P2 += length(commonanswers)
end

println("P1: ", P1)
println("P2: ", P2)

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
abstract type FusionMethod end

"""
    Mean <: FusionMethod

Use the mean to combine annotations.

# Constructor

    function Mean(; f = x -> x, g = x -> x)

Construct a Mean method. f is the function applied when taking the mean. For example:

```julia
julia> using Statistics

julia> f = x -> √x; mean(f, [1, 2, 3])
1.3820881233139908
```

g is the function applied after taking the mean. For example:

```julia
julia> using Statistics

julia> g = x -> round(x); round(g([1, 1, 3]))
2.0
```
"""
struct Mean <: FusionMethod
    f::Function # used in mean(f, x)
    g::Function # used after: g(mean(x))

    function Mean(; f = x -> x, g = x -> x)
        new(f, g)
    end
end

"""
    fuse(annotations::DataFrame, index::Symbol, method::Mean[; g::Function=x -> x])

Fuse annotations by applying the mean over annotations. This assumes an annotations matrix where each row contains
several annotations for an item or session.

# Parameters

 - g is a function to be applied _after_ the mean is computed. Can be used for rounding or casting.

# Examples
```julia
julia> annotations = DataFrame(items = string.(1:3), annotator1 = 1:3, annotator2 = 2:4, annotator3 = 3:5)
3×4 DataFrame
 Row │ items   annotator1  annotator2  annotator3
     │ String  Int64       Int64       Int64
─────┼────────────────────────────────────────────
   1 │ 1                1           2           3
   2 │ 2                2           3           4
   3 │ 3                3           4           5

julia> fuse(annotations, :items, Mean)
3×5 DataFrame
 Row │ items   annotator1  annotator2  annotator3  mean
     │ String  Int64       Int64       Int64       Float64
─────┼─────────────────────────────────────────────────────
   1 │ 1                1           2           3      2.0
   2 │ 2                2           3           4      3.0
   3 │ 3                3           4           5      4.0

julia> fuse(annotations, :items, Mean(); g = x -> round(Int, x))
3×5 DataFrame
 Row │ items   annotator1  annotator2  annotator3  mean
     │ String  Int64       Int64       Int64       Int64
─────┼───────────────────────────────────────────────────
   1 │ 1                1           2           3      2
   2 │ 2                2           3           4      3
   3 │ 3                3           4           5      4
```

"""
function fuse!(annotations::DataFrame, index::Symbol, method::Mean)
    annotations.mean = fuse(annotations, index, method)
    return annotations
end

function fuse(annotations::DataFrame, index::Symbol, method::Mean)
    return [method.g(mean(method.f, skipmissing(row[Not(index)]))) for row in eachrow(annotations)]
end

function fuse(annotations::AbstractMatrix, method::Mean)
    return [method.g(mean(method.f, skipmissing(row))) for row in eachrow(annotations)]
end

"""
    Median <: FusionMethod

Use the median to combine annotations.

# Constructor

    function Median(; g = x -> x)

Construct a Median method. g is the function applied after taking the mean. For example:

```julia
julia> using Statistics

julia> g = x -> round(x); g(median([1, 1, 3]))
1.0
```
"""
struct Median <: FusionMethod
    g::Function # Used after taking the median

    function Median(; g = x -> x)
        new(g)
    end
end

"""
    fuse(annotations::DataFrame, index::Symbol, method::Median)

Fuse annotations by applying the median over annotations. This assumes an annotations matrix where each row contains
several annotations for an item or session.

# Parameters

 - g is a function to be applied _after_ the mean is computed. Can be used for rounding or casting.

# Examples
```julia
julia> annotations = DataFrame(items = string.(1:3), annotator1 = 1:3, annotator2 = 2:4, annotator3 = 3:5)
3×4 DataFrame
 Row │ items   annotator1  annotator2  annotator3
     │ String  Int64       Int64       Int64
─────┼────────────────────────────────────────────
   1 │ 1                1           2           3
   2 │ 2                2           3           4
   3 │ 3                3           4           5

julia> fuse(annotations, :items, Median())
3×5 DataFrame
 Row │ items   annotator1  annotator2  annotator3  mean
     │ String  Int64       Int64       Int64       Float64
─────┼─────────────────────────────────────────────────────
   1 │ 1                1           2           3      2.0
   2 │ 2                2           3           4      3.0
   3 │ 3                3           4           5      4.0

julia> fuse(annotations, :items, Median(); g = x -> round(Int, x))
3×5 DataFrame
 Row │ items   annotator1  annotator2  annotator3  mean
     │ String  Int64       Int64       Int64       Int64
─────┼───────────────────────────────────────────────────
   1 │ 1                1           2           3      2
   2 │ 2                2           3           4      3
   3 │ 3                3           4           5      4
```

"""
function fuse!(annotations::DataFrame, index::Symbol, method::Median; g::Function=x -> x)
    annotations.mean = [g(median(skipmissing(row[Not(index)]))) for row in eachrow(annotations)]
    return annotations
end

function fuse(annotations::DataFrame, index::Symbol, method::Median; g::Function=x -> x)
    return [g(median(skipmissing(row[Not(index)]))) for row in eachrow(annotations)]
end

function fuse(annotations::AbstractMatrix, method::Median; g::Function=x -> x)
    return [g(median(skipmissing(row))) for row in eachrow(annotations)]
end

"""
    TE([loss::T=tSTE(α=30), ntriplets::Symbol = :all, verbose::Bool=true, print_every::Int=50]) where T <: TripletEmbeddings.AbstractLoss

Create a Triplet Embedding fusion method.

# Arguments

 - loss: Loss to be used (between STE and tSTE). Default is `tSTE(α=30)` (equivalent to `STE()` but more numerically stable)
 - ntriplets: Number of triplets to sample between `[:all, :auto]`. Use :auto if you run into memory issues, it randomly samples 20nlog(n) triplets.
 - scaling: Scaling of the embeddings. Options: `[:distribution, :procrustes]`. Use `:procrustes` for continuous scales, `:distribution` for discrete scales.

# Summary

struct TE <: AnnotationFusion.FusionMethod

# Fields

    loss        :: TripletEmbeddings.AbstractLoss
    ntriplets   :: Symbol
    verbose     :: Bool
    print_every :: Int64

# Supertype Hierarchy

TE <: AnnotationFusion.FusionMethod <: Any
"""
struct TE <: FusionMethod
    loss::T where T <: TripletEmbeddings.AbstractLoss
    ntriplets::Symbol
    constant::Int
    scaling::Symbol
    verbose::Bool
    print_every::Int

    function TE(; loss::T=tSTE(α=30), ntriplets::Symbol = :all, constant::Int = 40, scaling::Symbol=:procrustes, verbose::Bool=true, print_every::Int=50) where T <: TripletEmbeddings.AbstractLoss
        ntriplets in [:all, :auto] || throw(ArgumentError("ntriplets must be one of [:auto, :all]"))
        constant ≥ 1 || throw(ArgumentError("Constant in Cnlog(n) must be ≥ 1."))
        scaling in [:distribution, :procrustes] || throw(ArgumentError("scaling must be one of [:distribution, :procrustes]"))
        new(loss, ntriplets, constant, scaling, verbose, print_every)
    end
end

function fuse(annotations::DataFrame, index::Symbol, method::TE)
    n = size(annotations, 1)

    ntriplets = if method.ntriplets == :auto
        # If auto, uses the min between all triplets and method.constant * nlog(n)
        floor(Int, min(size(annotations, 1) * binomial(size(annotations, 1) - 1, 2), method.constant * size(annotations, 1) * log(size(annotations, 1))))
    elseif method.ntriplets == :all
        n * binomial(n - 1, 2)
    end

    triplets = Triplets(annotations, index, ntriplets)

    μ = fuse(annotations, index, Mean())
    X = Embedding(μ) # set the initial condition as the mean, since they are somewhat close

    misclassifications = TripletEmbeddings.fit!(method.loss, triplets, X; verbose=method.verbose, print_every=method.print_every)

    if method.scaling == :procrustes
        X, tr = procrustes(X, Matrix(μ'))
        return TripletEmbeddings.Vector(X)
    elseif method.scaling == :distribution
        return rankings_to_ratings(annotations, index, μ, TripletEmbeddings.Vector(X))
    end
end

function fuse!(annotations::DataFrame, index::Symbol, method::TE)
    annotations.TE = fuse(annotations, index, method)
end

struct Copeland <: FusionMethod
    imputation::Function
    scaling::Symbol

    """
        function Copeland([imputation::Function = row -> mean(skipmissing(row)), scaling::Symbol=:distribution])

    Create a Copeland's method struct.

    # Arguments

     - imputation: Function to impute values (per row). Defaults to the mean of the row values.
     - scaling: Either by reusing the ratings' distribution of the annotations matrix (`:distribution`) or `:procrustes`. Defaults to `:distribution`.
    """
    function Copeland(; imputation::Function = row -> mean(skipmissing(row)), scaling::Symbol=:distribution)
        scaling in [:distribution, :procrustes] || throw(ArgumentError("scaling mus be one of [:distribution, :procrustes]"))

        new(imputation, scaling)
    end
end

function copeland(annotations::DataFrame)
    points = zeros(size(annotations,1))
    nannotations, nannotators = size(annotations)

    # For each annotator, we compare all of their annotations with each other
    for i in 1:nannotations, j = i + 1:nannotations
        votesᵢ = 0 # Votes for each option
        votesⱼ = 0

        for a in 1:nannotators
            ratingᵢ = annotations[i, a]
            ratingⱼ = annotations[j, a]

             if !ismissing(ratingᵢ) || !ismissing(ratingⱼ)
                 if !ismissing(ratingᵢ) && ismissing(ratingⱼ)
                     ratingⱼ = 0
                 elseif ismissing(ratingᵢ) && ismissing(ratingⱼ)
                     ratingᵢ = 0
                 end

                if !ismissing(ratingᵢ) && !ismissing(ratingⱼ)
                    if ratingᵢ > ratingⱼ
                        votesᵢ = votesᵢ + 1
                    elseif ratingᵢ < ratingⱼ
                        votesⱼ = votesⱼ + 1
                    end
                end
             end
        end

        if votesᵢ > votesⱼ
            points[i] = points[i] + 1
        elseif votesᵢ < votesⱼ
            points[j] = points[j] + 1
        else
            points[i] = points[i] + 1/2
            points[j] = points[j] + 1/2
        end
    end

    return points
end

function fuse(annotations::DataFrame, index::Symbol, method::Copeland)
    μ = fuse(annotations, index, Mean())
    imputed = fillmissing(annotations, index, method.imputation)
    points = copeland(imputed[:, Not(index)])

    if method.scaling == :procrustes
        points, tr = procrustes(Matrix(points'), Matrix(μ'))
        return Vector(dropdims(points', dims=2))
    elseif method.scaling == :distribution
        return rankings_to_ratings(annotations, index, μ, points)
    end
end

function fuse!(annotations::DataFrame, index::Symbol, method::Copeland)
    annotations.copeland = fuse(annotations, index, method)
end

function fuse(annotations::AbstractMatrix, method::Copeland)
    μ = fuse(annotations, Mean())
    imputed = fillmissing(annotations, method.imputation)

    if method.scaling == :procrustes
        points, tr = procrustes(Matrix(points'), Matrix(μ'))
        return Vector(dropdims(points', dims=2))
    elseif method.scaling == :distribution
        return rankings_to_ratings(annotations, index, μ, points)
    end
end

function rankings_to_ratings(annotations::DataFrame, index::Symbol, means::Vector{T}, points::Vector{T}) where T <: Real
    size(annotations, 1) == length(points) || throw(ArgumentError("Number of rows in annotations must be the same as length of points"))

    results = sort(DataFrame([annotations[!,index],  points], [index, :points]), :points)

    vals = collect(skipmissing(vec(Matrix(annotations[:, Not(index)]))))
    scale = sort(unique(vals))

    counts = fit(Histogram, vals, nbins=length(scale)).weights
    ratings = zeros(Float64, length(points))
    args = ceil.(Int, cumsum(counts) ./ sum(counts) .* length(ratings))

    ratings[1:args[1]] .= scale[1]
    for i in 2:length(scale)
        ratings[args[i-1] + 1:args[i]] .= scale[i]
    end

    results.ratings = ratings
    return sort(results, index).ratings
end


text:
using Caesar, RoME, IncrementalInference
using Gadfly, KernelDensityEstimatePlotting, Cairo, Fontconfig
using KernelDensityEstimate
using DocStringExtensions
using DelimitedFiles, JLD, DistributedFactorGraphs

include(joinpath(@__DIR__,"slamUtils.jl"))
include(joinpath(@__DIR__,"plotSASUtils.jl"))
include(joinpath(@__DIR__,"expDataUtils.jl"))

function main(expID::String, datastart::Int, dataend::Int, fgap::Int, gps_gap::Int; saswindow::Int=8, trialID::Int=1, debug::Bool=false, dataloc::String="lj")
    ## Default parameters
    println("Start Script\n")

    sas_gap = saswindow+fgap;
    tstart = 1_000_000;
    element = 2

    totalposes = dataend-datastart;
    datawindow = collect(datastart:dataend);
    allpaths = getPaths(expID,"dynsas", currloc = dataloc, trialID = trialID);
    scriptHeader = joinpath(allpaths[3],expID*"_fgap$(fgap)_gpsgap$(gps_gap)_swind$(saswindow)_f$(datastart)t$(dataend)_")
    savefgHeader = joinpath(allpaths[3]*"_fg",expID*"_fgap$(fgap)_gpsgap$(gps_gap)_swind$(saswindow)_f$(datastart)t$(dataend)_")

    poses = Dict{Int,Array}();
    nav = Dict{Int,Array}();
    errorinds = [];
    allsasframes = Dict{Int,Array}();

    igt = loadGT(datawindow,expID,allpaths);
    posData = loadNav(datawindow,allpaths)
    dposData = deepcopy(posData)
    i=1
    while i+gps_gap < size(dposData,1)
        cumulativeDrift!(@view(dposData[i:i+gps_gap-1,:]),[0.0;0],[0.2,0.2])
         i+=gps_gap
        if i+gps_gap > size(dposData,1)
            cumulativeDrift!(@view(dposData[i-gps_gap:end,:]),[0.0;0],[0.2,0.2])
        end
    end

    # posPl = Gadfly.plot(layer(x=posData[:,1],y=posData[:,2], Geom.path, Theme(default_color=colorant"green")), layer(x=dposData[:,1],y=dposData[:,2],Geom.path)); posPl |> PDF(scriptHeader*"posGT.pdf")

    chirpFile = joinpath("/media","data1","data","kayaks","chirp250.txt");
    cfgFile = joinpath("/media","data1","data","kayaks","SAS2D.yaml")
    cfgd=loadConfigFile(cfgFile)

    fg = initfg();
    getSolverParams(fg).maxincidence = 1000

    beacon = :l1
    addVariable!(fg, beacon, Point2 )
    tree = BayesTree();

    pose_counter = 1
    sas_counter = 1
    sas_gap_counter = 0;

    while pose_counter < totalposes-2
        current_symbol = Symbol("x$pose_counter")
        addVariable!(fg, current_symbol, DynPoint2(ut=tstart))
        println("Adding Variable x$pose_counter")
        sas_gap_counter += 1
        tstart += 1_000_000

        #Interpolate velocity from GPS ground truth every gps_gap poses
       if pose_counter == 1 || mod(pose_counter,gps_gap) == 0
           xdotp = posData[pose_counter+1,1] - posData[pose_counter,1];
           ydotp = posData[pose_counter+1,2] - posData[pose_counter,2];
            dpμ = [posData[pose_counter,1];posData[pose_counter,2];xdotp;ydotp];
            dpσ = Matrix(Diagonal([0.1;0.1;0.2;0.2].^2))

            pp = DynPoint2VelocityPrior(MvNormal(dpμ,dpσ))
            addFactor!(fg, [current_symbol;], pp, autoinit=true)
            # initManual!(fg,current_symbol,kde!(rand(MvNormal(dpμ,dpσ),100)))
       end

       #odo factors
       if pose_counter > 1
           xdotp = dposData[pose_counter,1] - dposData[pose_counter-1,1];
           ydotp = dposData[pose_counter,2] - dposData[pose_counter-1,2];
           xdotf = dposData[pose_counter+1,1] - dposData[pose_counter,1];
           ydotf = dposData[pose_counter+1,2] - dposData[pose_counter,2];
           # dpμ = [xdotp;ydotp;xdotf-xdotp;ydotf-ydotp];
           dpμ = [xdotp;ydotp;0;0];
           dpσ = Matrix(Diagonal([0.2;0.2;0.2;0.2].^2))
           vp = VelPoint2VelPoint2(MvNormal(dpμ,dpσ))
           addFactor!(fg, [Symbol("x$(pose_counter-1)");current_symbol], vp, autoinit=true)
           # initManual!(fg,current_symbol,kde!(rand(MvNormal(dpμ,dpσ),100)))
       end

       if sas_gap_counter >= sas_gap
           println("Trying to add SAS at pose: $pose_counter")
           navframes = collect(pose_counter-saswindow+1:1:pose_counter)
           navtemp = posData[navframes,:];
           navchecked, errorind = sanitycheck_nav(navtemp);
           push!(errorinds,errorind);

           if navchecked
               poses[sas_counter] = [Symbol("x$i") for i in pose_counter-saswindow+1:pose_counter]
               currentstart = datastart+pose_counter-saswindow
               println("Adding SAS at dataframe: $currentstart")
               sasframes = collect(currentstart:1:currentstart+saswindow-1);
               allsasframes[sas_counter] = sasframes;
               waveformData = importdata_waveforms(sasframes,element,datadir=allpaths[1]);
               sas2d = prepareSAS2DFactor(saswindow, waveformData, rangemodel=:Correlator,cfgd=cfgd, chirpFile=chirpFile)
               addFactor!(fg, [beacon;poses[sas_counter]], sas2d, autoinit=false)

               getSolverParams(fg).drawtree = false
               getSolverParams(fg).showtree = false
               getSolverParams(fg).limititers = 500

               lstpose = Symbol("x$(pose_counter)");
               tmpInit = approxConv(fg,Symbol("x$(pose_counter-1)x$(pose_counter)f1"),lstpose);
               XXkde = manikde!(tmpInit,getManifolds(fg,lstpose));
               setValKDE!(fg,lstpose,XXkde);

               if sas_counter > 1
                   tree, smt, hist = solveTree!(fg,tree)
               else
                   tree, smt, hist = solveTree!(fg)
               end

               writeGraphPdf(fg,viewerapp="", engine="neato", filepath=scriptHeader*"fg.pdf")
               drawTree(tree, filepath=scriptHeader*"bt.pdf")

               plotSASDefault(fg,expID, posData,igt,dposData,savedir=scriptHeader*"$sas_counter.pdf")

               if debug
                   plk = [];
                   push!(plk,plotBeaconGT(igt));
                   plotBeaconContours!(plk,fg);
                   for mysym in poses[sas_counter]
                       xData = getKDEMax(getBelief(getVariable(fg,mysym)))
                       push!(plk, plotPoint(xData,colorIn=colorant"orange"))
                   end
                   L1p = approxConv(fg,ls(fg,:l1)[end],:l1);
                   push!(plk,layer(x=L1p[1,:],y=L1p[2,:],Geom.histogram2d(xbincount=300, ybincount=300)))
                   Gadfly.plot(plk...) |> PDF(scriptHeader*"debug$sas_counter.pdf");
               end

               if expID == "dock"
                   l1fit = getKDEMean(getBelief(getVariable(fg,:l1)))
                   meanerror = sqrt((igt[1]-l1fit[1])^2+(igt[2]-l1fit[2])^2)

                   l1max = getKDEMax(getBelief(getVariable(fg,:l1)))
                   maxerror = sqrt((igt[1]-l1max[1])^2+(igt[2]-l1max[2])^2)

                   mykde = getBelief(getVariable(fg,:l1))
                   mynorm = kde!(rand(fit(MvNormal,getVal(fg,:l1)),200))
                   kld = min(abs(KernelDensityEstimate.kld(mynorm,mykde)),abs(KernelDensityEstimate.kld(mykde,mynorm)))
               elseif expID == "drift"
                   l1fit = getKDEMean(getBelief(getVariable(fg,:l1)))
                   meanerror = sqrt((mean(igt[:,1])-l1fit[1])^2+(mean(igt[:,2])-l1fit[2])^2)

                   l1max = getKDEMax(getBelief(getVariable(fg,:l1)))
                   maxerror = sqrt((mean(igt[:,1])-l1max[1])^2+(mean(igt[:,2])-l1max[2])^2)

                   mykde = getBelief(getVariable(fg,:l1))
                   mynorm = kde!(rand(fit(MvNormal,getVal(fg,:l1)),200))
                   kld = min(abs(KernelDensityEstimate.kld(mynorm,mykde)),abs(KernelDensityEstimate.kld(mykde,mynorm)))
               end

               ev = 0;
               evi = 0;
               es = 0;
               for tmpi = 1:pose_counter
                   rv = getVal(fg,Symbol("x$tmpi"));
                   dxt = (rv[1,:].-posData[tmpi,1]).^2;
                   dyt = (rv[2,:].-posData[tmpi,2]).^2;
                   ev += sum(sqrt.(dxt+dyt));
                   dxt2 = (dposData[tmpi,1].-posData[tmpi,1]).^2;
                   dyt2 = (dposData[tmpi,2].-posData[tmpi,2]).^2;
                   evi += sum(sqrt.(dxt2+dyt2));

                   if tmpi > 1
                       dv = (rv[3,:].-(posData[tmpi,1]-posData[tmpi-1,1])).^2;
                       dw = (rv[4,:].-(posData[tmpi,2]-posData[tmpi-1,2])).^2;
                       es += sum(sqrt.(dv+dw));
                   end
               end
               ev = ev./pose_counter;
               evi = evi./pose_counter;
               es = es./(pose_counter-1);

               # jldname2 = scriptHeader * "solve_$(sas_counter).jld"
               # JLD.save(jldname2,"beacon",getVal(fg,:l1),"posData",posData,"dposData", dposData,"gps_gap", gps_gap, "poses",poses,"sasframes", allsasframes, "l1fit",l1fit, "meanerror",meanerror,"l1max",l1max,"maxerror",maxerror,"kld",kld,"ev",ev,"evi",ev, "es", es)

               # saveDFG(fg,savefgHeader * "fg$(sas_counter)")

               writedlm(scriptHeader*"stats$(sas_counter).txt", [pose_counter meanerror maxerror kld evi ev es], ",")

               sas_counter +=1
               sas_gap_counter = 0
           end
       end
       pose_counter+=1

       println("Next Variable is $pose_counter")
       println("Next SAS is $sas_counter")

   end

   println("This Solve was saved under: "*scriptHeader)
   return fg
end

# main("drift", 1550,1580,2,49)
#
# fg = main("drift", 1550,1580,5,49)
#
# main("drift", 1550,1600,5,49)
#
# fg = main("drift", 1550,1800,5,49)
#

# fg = main("dock", 410,700,22,120)

# fg = main("dock", 410,700,40,120)


text:
module RocketCountOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: count()" begin

    println("Testing: operator count()")

    run_proxyshowcheck("Count", count())

    run_testset([
        (
            source = from(1:42) |> count(),
            values = @ts([ 42, c ])
        ),
        (
            source = from(1:42) |> async(0) |> count(),
            values = @ts([ 42, c ])
        ),
        (
            source = completed(Int) |> count(),
            values = @ts([ 0, c ])
        ),
        (
            source = never(Int) |> count(),
            values = @ts()
        ),
        (
            source = faulted("e") |> count(),
            values = @ts(e("e"))
        )
    ])

end

end


text:
abstract type IsingMeasurement <: AbstractMeasurement end

"""
    IsingEnergyMeasurement(mc::MC, model::IsingModel, rate=1)

Measures observables related to the energy of the IsingModel. This includes
- `E`: the total energy
- `E2`: the squared energy
- `e`: the energy per spin
- `C`: the specific heat
"""
struct IsingEnergyMeasurement <: IsingMeasurement
    invN::Float64
    E::Observable
    E2::Observable
    e::Observable
    C::Observable

    IsingEnergyMeasurement(mc::MC, model::IsingModel) = new(
        1.0 / model.l.sites,
        Observable(Float64, "Total energy"),
        Observable(Float64, "Total energy squared"),
        Observable(Float64, "Energy per site"),
        Observable(Float64, "Specific heat")
    )
end

function measure!(m::IsingEnergyMeasurement, mc::MC, model::IsingModel, i::Int64)
    push!(m.E, model.energy[])
    push!(m.E2, model.energy[]^2)
    push!(m.e, model.energy[] * m.invN)
    nothing
end

function finish!(m::IsingEnergyMeasurement, mc::MC, model::IsingModel)
    E = mean(m.E)
    E2 = mean(m.E2)
    push!(m.C, mc.p.beta^2 * m.invN * (E2 - E^2))
    nothing
end

"""
    IsingMagnetizationMeasurement(mc::MC, model::IsingModel, rate=1)

Measures observables related to the magnetization of the IsingModel. This
includes
- `M`: the total magnetization
- `M2`: the squared magnetization
- `m`: the magnetization per spin
- `chi`: the magentic susceptibility
"""
struct IsingMagnetizationMeasurement <: IsingMeasurement
    invN::Float64
    M::Observable
    M2::Observable
    m::Observable
    chi::Observable

    IsingMagnetizationMeasurement(mc::MC, model::IsingModel) = new(
        1.0 / model.l.sites,
        Observable(Float64, "Total magnetization"),
        Observable(Float64, "Total magnetization squared"),
        Observable(Float64, "Magnetization per site"),
        Observable(Float64, "Magnetic susceptibility")
    )
end

function measure!(m::IsingMagnetizationMeasurement, mc::MC, model::IsingModel, i::Int64)
    M = abs(sum(mc.conf))
    push!(m.M, M)
    push!(m.M2, M^2)
    push!(m.m, M * m.invN)
    nothing
end

function finish!(m::IsingMagnetizationMeasurement, mc::MC, model::IsingModel)
    M = mean(m.M)
    M2 = mean(m.M2)
    push!(m.chi, mc.p.beta * m.invN * (M2 - M^2))
    nothing
end

function default_measurements(mc::MC, model::IsingModel)
    Dict(
        :conf => ConfigurationMeasurement(mc, model),
        :Magn => IsingMagnetizationMeasurement(mc, model),
        :Energy => IsingEnergyMeasurement(mc, model)
    )
end


text:
using SampledDomains
using Test

@testset "SampledDomains.jl" begin
    a= SampledDomains.CartesianDomain2D(1:3, -5:.5:-3.5)
    f(x,y) = x*y
    b = SampledDomains.SampledDomain(f, a)
    @test b.vals ==[
        -5.0  -10.0  -15.0
        -4.5   -9.0  -13.5
        -4.0   -8.0  -12.0
        -3.5   -7.0  -10.5
    ]
    @test b.dom === a
end


text:
module cthmisc
	using DataFrames
	export merge_columns!,indexer
	include("dataframes.jl")

	using Distributions
	export inverse_normal_transform, mean_imputation!, group_mean_imputation!
	include("statistics.jl")

	export  joincmds,mkdir_ifnot,parse_danish_float
	include("other.jl")
end # module


text:
# ---
# title: 11. Container With Most Water
# id: problem11
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/container-with-most-water/description/>
# hidden: true
# ---
#
# Given `n` non-negative integers `a1, a2, ..., an` , where each represents a
# point at coordinate `(i, ai)`. `n` vertical lines are drawn such that the two
# endpoints of the line `i` is at `(i, ai)` and `(i, 0)`. Find two lines, which,
# together with the x-axis forms a container, such that the container contains
# the most water.
#
# **Notice**  that you may not slant the container.
#
#
#
# **Example 1:**
#
# ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/07/17/question_11.jpg)
#
#
#
#     Input: height = [1,8,6,2,5,4,8,3,7]
#     Output: 49
#     Explanation: The above vertical lines are represented by array [1,8,6,2,5,4,8,3,7]. In this case, the max area of water (blue section) the container can contain is 49.
#
#
# **Example 2:**
#
#
#
#     Input: height = [1,1]
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: height = [4,3,2,1,4]
#     Output: 16
#
#
# **Example 4:**
#
#
#
#     Input: height = [1,2,1]
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `n = height.length`
#   * `2 <= n <= 3 * 104`
#   * `0 <= height[i] <= 3 * 104`
#
#
## @lc code=start
using LeetCode

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


text:
using Compat
using BinDeps

const CXXWRAP_JL_VERSION = v"0.5.0+"

function prompt_cmake()
  try
    run(pipeline(`cmake --version`, stdout=DevNull, stderr=DevNull))
    return
  catch
    error("command \"cmake\" not found in path, please install CMake and make sure it can be found")
  end
end

@static if is_windows()
  # prefer building if requested
  if haskey(ENV, "BUILD_ON_WINDOWS") && ENV["BUILD_ON_WINDOWS"] == "1"
    prompt_cmake()
    saved_defaults = deepcopy(BinDeps.defaults)
    empty!(BinDeps.defaults)
    append!(BinDeps.defaults, [BuildProcess])
  end
else
  prompt_cmake()
end

@BinDeps.setup

build_type = get(ENV, "JLCXX_BUILD_TYPE", "Release")
build_examples = get(ENV, "JLCXX_BUILD_EXAMPLES", "ON")

# List of the libraries that will be built
lib_labels = [:jlcxx]
if build_examples == "ON"
  lib_labels = vcat(lib_labels, [:jlcxx_containers, :except, :extended, :functions, :hello, :inheritance, :parametric, :types])
end
deps = BinDeps.LibraryDependency[]
for l in lib_labels
  @eval $l = $(library_dependency(string(l), aliases=["lib"*string(l)]))
  push!(deps, eval(:($l)))
end

depsdir = BinDeps.depsdir(first(deps))
prefix = joinpath(depsdir, "usr")

@static if is_windows()
    bindir = joinpath(prefix, "bin")
else
    bindir = joinpath(prefix, "lib")
end

jlcxx_srcdir = joinpath(depsdir, "src", "jlcxx")
jlcxx_builddir = joinpath(depsdir, "builds", "jlcxx")
lib_prefix = @static is_windows() ? "" : "lib"
lib_suffix = @static is_windows() ? "dll" : (@static is_apple() ? "dylib" : "so")
julia_base_dir = splitdir(JULIA_HOME)[1]
julia_executable = split(string(Base.julia_cmd()))[1][2:end]

makeopts = ["--", "-j", "$(Sys.ARCH == :armv7l ? 2 : Sys.CPU_CORES+2)"]

# Set generator if on windows
genopt = "Unix Makefiles"
@static if is_windows()
  if isdir(prefix)
    rm(prefix, recursive=true)
  end
  if get(ENV, "MSYSTEM", "") == ""
    makeopts = "--"
    if Sys.WORD_SIZE == 64
      genopt = "Visual Studio 14 2015 Win64"
    else
      genopt = "Visual Studio 14 2015"
    end
  else
    lib_prefix = "lib" #Makefiles on windows do keep the lib prefix
  end
end

deps_paths = AbstractString[]
for l in lib_labels
  push!(deps_paths, joinpath(bindir, "$(lib_prefix)$(string(l)).$lib_suffix"))
end

cxx_steps = @build_steps begin
  `cmake -G "$genopt" -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_BUILD_TYPE="$build_type" -DCMAKE_PROGRAM_PATH=$JULIA_HOME -DJLCXX_BUILD_EXAMPLES=$build_examples -DCMAKE_INSTALL_LIBDIR=lib $jlcxx_srcdir`
  `cmake --build . --config $build_type --target install $makeopts`
end

# If built, always run cmake, in case the code changed
if isdir(jlcxx_builddir)
  BinDeps.run(@build_steps begin
    ChangeDirectory(jlcxx_builddir)
    cxx_steps
  end)
end

provides(BuildProcess,
  (@build_steps begin
    CreateDirectory(jlcxx_builddir)
    @build_steps begin
      ChangeDirectory(jlcxx_builddir)
      FileRule(deps_paths, cxx_steps)
    end
  end), deps)

@static if is_windows()
  shortversion = "$(VERSION.major).$(VERSION.minor)"
  zipfilename = "CxxWrap-julia-$(shortversion)-win$(Sys.WORD_SIZE).zip"
  archname = Sys.WORD_SIZE == 64 ? "x64" : "x86"
  pkgverstring = string(CXXWRAP_JL_VERSION)
  if endswith(pkgverstring,"+")
    bin_uri = URI("https://ci.appveyor.com/api/projects/barche/cxxwrap-jl/artifacts/$(zipfilename)?job=Environment%3a+JULIAVERSION%3djulialang%2fbin%2fwinnt%2f$(archname)%2f$(shortversion)%2fjulia-$(shortversion)-latest-win$(Sys.WORD_SIZE).exe%2c+BUILD_ON_WINDOWS%3d1")
  else
    bin_uri = URI("https://github.com/JuliaInterop/CxxWrap.jl/releases/download/v$(pkgverstring)/CxxWrapv$(pkgverstring)-julia-$(VERSION.major).$(VERSION.minor)-win$(Sys.WORD_SIZE).zip")
  end
  provides(Binaries, Dict(bin_uri => deps), os = :Windows)
end

if build_examples == "ON"
  @BinDeps.install Dict([(:jlcxx, :_l_jlcxx),
                         (:jlcxx_containers, :_l_containers),
                         (:except, :_l_except),
                         (:extended, :_l_extended),
                         (:functions, :_l_functions),
                         (:hello, :_l_hello),
                         (:inheritance, :_l_inheritance),
                         (:parametric, :_l_parametric),
                         (:types, :_l_types)])
else
  @BinDeps.install Dict([(:jlcxx, :_l_jlcxx)])
end

@static if is_windows()
  if haskey(ENV, "BUILD_ON_WINDOWS") && ENV["BUILD_ON_WINDOWS"] == "1"
    empty!(BinDeps.defaults)
    append!(BinDeps.defaults, saved_defaults)
    if get(ENV, "MSYSTEM", "") == "MINGW32"
      run(`cp -f $(joinpath("/mingw32", "bin", "libwinpthread-1.dll")) $(bindir)`)
      run(`cp -f $(joinpath("/mingw32", "bin", "libgcc_s_dw2-1.dll")) $(bindir)`)
    else
      redist_dlls = ["concrt140.dll", "msvcp140.dll", "vccorlib140.dll", "vcruntime140.dll"]
      redistbasepath = "C:\\Program Files (x86)\\Microsoft Visual Studio 14.0\\VC\\redist\\$(Sys.WORD_SIZE == 64 ? "x64" : "x86")\\Microsoft.VC140.CRT"
      for dll in redist_dlls
        cp(joinpath(redistbasepath, dll), joinpath(bindir, dll), remove_destination=true)
      end
    end
  end
end


text:
module AbInitioTopo

export compute_Z2, compute_Z4

include("Calculators.jl")
include("Parsers.jl")

end


text:
module PsychrometricCharts

# package code goes here

using Psychro

function calcpoints(Pa=101325.0, Tlim=[0.0, 50.0], nT=100,
                    Wlim=[0.0, 30.0], nW=100)

    Tmax = maximum(Tlim) + 273.15
    Tmin = minimum(Tlim) + 273.15
    dT = (Tmax - Tmin) / (nT-1)
    T = Tmin:dT:Tmax

    Wmax = maximum(Wlim)/1000
    Wmin = minimum(Wlim)/1000
    dW = (Wmax - Wmin) / (nW-1)

    h = zeros(nW, nT)
    s = zeros(nW, nT)
    r = zeros(nW, nT)
    Td = zeros(nW, nT)
    w = zeros(nW, nT)
    q = zeros(nW, nT)
    ρ = zeros(nW, nT)
    Tb = zeros(nW, nT)
    xv = zeros(nW, nT)
    q = zeros(nW, nT)
    Z = zeros(nW, nT)

    wsat = humrat.(MoistAir, T, RelHum, 1.0, Pa)

    for i = 1:nT
        w1 = linspace(Wmin, max(wsat[i], Wmax), nW)
        for k = 1:nW
            h[i,k] = enthalpy(MoistAir, T[i], HumRat, w1[k], Pa)
            s[i,k] = entropy(MoistAir, T[i], HumRat, w1[k], Pa)
            r[i,k] = relhum(MoistAir, T[i], HumRat, w1[k], Pa)

        end
    end
end

function pretty(x, n=8, dx=[1.0, 2.0, 5.0])
    xmin = minimum(x)
    xmax = maximum(x)
    xr = xmax-xmin
    dx0 = xr / n
    logdx = [log10(y) - floor(log10(y))  for y in dx]
    pwr = floor(log10(dx0))
    l0 = log10(dx0) - pwr

    val, idx = findmin(abs.(abs(l0) .- logdx))
    #return val, idx, dx[idx], (l0, pwr)
    dx1 = dx[idx]*10^pwr

    xlim1 = round(xmin/dx1, RoundDown)*dx1
    xlim2 = round(xmax/dx1, RoundUp)*dx1

    return xlim1:dx1:xlim2

end

function drawlines(humfun, HumType, n=10, Tc, W, Pa=101325.0)

    Tmin = minimum(Tc) + 273.15
    Tmax = maximum(Tc) + 273.15
    Wmin = minimum(W)
    Wmax = maximum(W)

    wsat1 = humrat(MoistAir, Tmin, RelHum, 1.0, Pa)
    wsat2 = humrat(MoistAir, Tmax, RelHum, 1.0, Pa)

    wmin = min(Wmin, wsat1)
    wmax = min(wsat2, Wmax)

    p1 = humfun(MoistAir, Tmin, HumRat, wmin, Pa)
    p2 = humfun(MoistAir, Tmin, HumRat, wmax, Pa)
    p3 = humfun(MoistAir, Tmax, HumRat, wmin, Pa)
    p4 = humfun(MoistAir, Tmax, HumRat, wmax, Pa)

    pmin = min(p1, p2, p3, p4)
    pmax = max(p1, p2, p3, p4)

end
function relhumlines(Tc, r, Pa)
    nT = size(T,1)
    nr = size(r,1)

    w = zeros(nT, nr)
    h = zeros(nT, nr)

    T0 = 273.15

    for i in 1:nr
        for k in 1:nT
            w[k,i] = humrat(MoistAir, Tc[k] + T0, RelHum, r[i], Pa)
            h[k,i] = enthalpy(MoistAir, T[k] + T0, RelHum, r[i], Pa)
        end
    end

    return w, h
end

end # module


text:
struct Agent{T,E}
    id::Int
    values::Vector{T}
    vote::Ref{Int}
    energy::Ref{E}
end

function initNetwork(data, N)
    #Generate Barabasi graph with N nodes, 3 conntections each, 10 initial nodes
    network = barabasi_albert(N, 10, 3, #=seed=1,=# is_directed=true)
    #global Network = erdos_renyi(N, 4, is_directed=true)

    #Initialize agents
    nodes = Agent{Int,Int}[]
    for i in 1:N
        push!(nodes, Agent(i, rand(0:10, 5), Ref(rand(1:10)), Ref(0)))
    end

    #Initialize energy
    for i in 1:N
        dE!(nodes, i, network)
    end
    computeEnergy!(data, nodes)

    #Count initial preferences
    trackPreference!(data, nodes)

    return network, nodes
end

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
"""
    assert(x::Bool)

Assert that the condition generating `x` holds true.
Returns `x`

"""
assert(x::Bool) = x

"""
    prove(x::Bool)

```julia
assert(x < 10)
y = x - 10
prove(y < 0)
```

Will turn into the corresponding SMT2 program:

```
(declare-const x Int)
(declare-const y Int)
(assert (< x 10))
(assert (= y (- x 10)))
(assert (not (< y 0)))
(check-sat)
```

If the model is not satisfiable, given the pre-condition defined by `assert`
the post-condition given by `prove` holds.

For more information see `check`.
"""
prove(x::Bool) = x

nodub(f, args...) = Expr(:call, Expr(:nooverdub, f), args...)

# Define a Cassette pass to insert asserts after branches
function insertasserts(ctx, ref)
    ir = ref.code_info
    locations = Int[]
    Cassette.insert_statements!(ir.code, ir.codelocs,
        (stmt, i) -> Base.Meta.isexpr(stmt, :gotoifnot) ? 2 : nothing,
        (stmt, i) -> [
            Expr(:call, ConcolicFuzzer.assert, stmt.args[1]),
            stmt
        ])

    # Taint foreigncalls
    # TODO: llvmcall, check if stmts.args[2] supported
    # This works but is way to aggressive.
    #=
    Cassette.insert_statements!(ir.code, ir.codelocs,
        (stmt, i) -> begin
            stmt = Base.Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
            if Base.Meta.isexpr(stmt, :foreigncall) && supported(stmt.args[2])
                12
            else
                nothing
            end
        end,
        (stmt, i) -> begin
            if Base.Meta.isexpr(stmt, :(=))
                tagstmt = Expr(:(=), stmt.args[1])
                stmt = stmt.args[2]
            else
                tagstmt = nothing
            end
            # setup slot for if-else
            slot = Core.SlotNumber(length(ir.slotflags) + 1)
            push!(ir.slotflags, 0x10) # TODO check flag value
            push!(ir.slotnames, gensym(:temp))
            # obtain location as stacktrace
            type = stmt.args[2]

            # Helper to more easily create result
            result = Any[]
            insert!(stmt) = (id = i+length(result); push!(result, stmt); Core.SSAValue(id))

            #=
            stack = StackTraces.stacktrace()
            loc = StackTraces.remove_frames!(stack, :execute)
            val, sym = record!(ctx, loc, type)
            if val === nothing
                val = stmt
            end
            return tag(val, ctx, sym)
            =#
            stack  = insert!(nodub(StackTraces.stacktrace))
            loc    = insert!(nodub(StackTraces.remove_frames!, stack, QuoteNode(:execute)))
            record = insert!(nodub(ConcolicFuzzer.record!, Expr(:contextslot), loc, type))
            val    = insert!(nodub(Base.getindex, record, 1))
            sym    = insert!(nodub(Base.getindex, record, 2))
            cond   = insert!(nodub(Core.:(===), val, nothing))
            ifstmt = insert!(Expr(:gotoifnot, cond, #=fill in later=#))
            bb1    = insert!(Expr(:(=), slot, stmt))
            elstmt = insert!(Expr(:goto, #=fill in later=#))
            bb2    = insert!(Expr(:(=), slot, val))
            tag    = insert!(nodub(Cassette.tag, slot , Expr(:contextslot), sym))

            # fixup if and else
            push!(result[ifstmt.id - i + 1].args, bb2.id)
            result[elstmt.id - i + 1] = Core.GotoNode(tag.id)

            if tagstmt === nothing
                push!(result, tag)
            else
                push!(tagstmt.args, tag)
                push!(result, tagstmt)
            end
            return result
        end)
    =#
    ir.ssavaluetypes = length(ir.code)
    return ir
end

const InsertAssertsPass = Cassette.@pass insertasserts

# We can't use primitives since otherwise our tracer won't trace assert and prove statements
# Which is the entire point.
function Cassette.overdub(ctx::TraceCtx, f::F, x::Tagged) where F <: Union{Core.typeof(assert), Core.typeof(prove)}
    call = Callsite(f, record((x, ), ctx))
    retval = untag(x, ctx)
    call.retval = Some(retval)
    push!(ctx.metadata, call)
    return retval
end


text:
include("histogram_test.jl")

# vars = 75.0:1.0:85.0                 # angle
vars = [-179.0:5.0:-1.0; 1.0:5.0:179.0] # angle
# vars = 1.0:5.0:200                   # mass
n_frames = length(vars)
anim = @animate for (i, var) in enumerate(vars)
    print("Frame: $i / $n_frames complete")
    fvar = Float64(var)
    energy_histogram(obs_angle = fvar)
end
gif(anim, "histgif.gif", fps = 20)