text:
__precompile__()
module SeisDownload

include("utils.jl")
include("downloadfunc.jl")

using .Utils
using .DownloadFunc

using SeisIO, Dates, Printf, JLD2, FileIO, Distributed

#------------------------------------------------------------------#
#For the time being, we need remove_response function from obspy
#This will be replaced by SeisIO modules in the near future.
#Please activate obspy enviroment before launching Julia.
include("remove_response_obspy.jl")
using .Remove_response_obspy
#------------------------------------------------------------------#

export seisdownload

"""
    seisdownload(InputDict::Dict; MAX_MEM_PER_CPU::Float64=1.0)

    Request seismic data and save into jld2.
# Arguments
- `InputDict`    : dictionary which contains request information
- `MAX_MEM_PER_CPU` : maximum available memory for 1 cpu (unit:GB) (default = 1.0GB)
"""
function seisdownload(InputDict::Dict; MAX_MEM_PER_CPU::Float64=1.0)

    Utils.initlogo()

	DownloadType    = InputDict["DownloadType"]

    if DownloadType == "Noise" || DownloadType == "noise"

		#stationlist
		stationlist     = InputDict["stationinfo"]["stationlist"]
		starttime       = InputDict["starttime"]
		endtime         = InputDict["endtime"]
		DL_time_unit    = InputDict["DL_time_unit"]
		DownloadType    = InputDict["DownloadType"]
		fopath          = InputDict["fopath"]

		if mod((endtime - starttime).value,  DL_time_unit) != 0 || (endtime - starttime).value < DL_time_unit
			error("Total download time cannot be devided by Download Time unit; this may cause unexpected result. abort.")
		end

		# calculate start time list (starttimelist) with each Donwload_time_unit
		starttimelist = Utils.get_starttimelist(starttime, endtime, DL_time_unit)
		# generate DLtimestamplist and ststationlist
		DLtimestamplist = Utils.get_timestamplist(starttimelist)

		InputDict["starttimelist"] = starttimelist
		InputDict["DLtimestamplist"] = DLtimestamplist

		#save info into jld2
		jldopen(fopath, "w") do file
			file["info/DLtimestamplist"] = DLtimestamplist;
			file["info/stationlist"] = stationlist;
			file["info/starttime"]   = string(starttime)
			file["info/endtime"]     = string(endtime)
			file["info/DL_time_unit"]= string(DL_time_unit)
		end

        #Test download to evaluate use of memory and estimate download time.
		max_num_of_processes_per_parallelcycle = testdownload(InputDict, length(starttimelist), MAX_MEM_PER_CPU)

        if max_num_of_processes_per_parallelcycle < 1
            error("Memory allocation is not enought (currently $MAX_MEM_PER_CPU [GB]). Please inclease MAX_MEM_PER_CPU or decrease number of stations")
        end

        if max_num_of_processes_per_parallelcycle >= length(starttimelist)

            S = pmap(x -> seisdownload_NOISE(x, InputDict), 1:length(starttimelist))

            # save data to jld2
            file = jldopen(fopath, "r+")
            unavalilablefile = jldopen(join([fopath[1:end-5], "_unavailablestations.jld2"]), "w+")

            for ii = 1:size(S)[1] #loop at each starttime
                for jj = 1:size(S[1])[1] #loop at each station id

                    requeststr =S[ii][jj].id
                    varname = joinpath(DLtimestamplist[ii], requeststr)
                    #save_SeisData2JLD2(fopath, varname, S[ii][jj])
                    file[varname] = S[ii][jj]

                    if S[ii][jj].misc["dlerror"] == 1
                        unavalilablefile[varname] = S[ii][jj]
                    end
                end
            end
            JLD2.close(file)
            JLD2.close(unavalilablefile)

        else

            #parallelization by time
    	    pitr = 1

    	    # progress bar

    	    while pitr <= length(starttimelist)

    	        startid1 = pitr
    	        startid2 = pitr + max_num_of_processes_per_parallelcycle - 1

    	        if startid1 == length(starttimelist)
					S = pmap(x -> seisdownload_NOISE(x, InputDict), startid1:startid1)

				elseif startid2 <= length(starttimelist)
    	            #use all processors
	                S = pmap(x -> seisdownload_NOISE(x, InputDict), startid1:startid2)

    	        else
    	            #use part of processors
	                startid2 = startid1 + mod(length(starttimelist), nprocs()) - 1
	                S = pmap(x -> seisdownload_NOISE(x, InputDict), startid1:startid2)
    	        end

                # save data to jld2
                file = jldopen(fopath, "r+")
                unavalilablefile = jldopen(join([fopath[1:end-5], "_unavailablestations.jld2"]), "w+")

                for ii = 1:size(S)[1] #loop at each starttime
                    for jj = 1:size(S[1])[1] #loop at each station id

                        requeststr =S[ii][jj].id
                        varname = joinpath(DLtimestamplist[startid1+ii-1], requeststr)
                        #save_SeisData2JLD2(fopath, varname, S[ii][jj])
                        file[varname] = S[ii][jj]

                        if S[ii][jj].misc["dlerror"] == 1
                            unavalilablefile[varname] = S[ii][jj]
                        end
                    end
                end

                JLD2.close(file)
                JLD2.close(unavalilablefile)

    	        pitr += max_num_of_processes_per_parallelcycle

                #println("pitr: $pitr")
            end
        end

    elseif  DownloadType == "Earthquake" || DownloadType == "earthquake"

		method		    = InputDict["method"]
		event		    = InputDict["event"]
		reg			    = InputDict["reg"]
		fopath          = InputDict["fopath"]

		#save info into jld2
		jldopen(fopath, "w") do file
			file["info/method"]  = method;
			file["info/event"]   = event;
			file["info/reg"]     = reg
			file["info/fopath"]  = fopath
		end

		#Test download to evaluate use of memory and estimate download time.
		max_num_of_processes_per_parallelcycle = testdownload(InputDict, length(event), MAX_MEM_PER_CPU)

		if max_num_of_processes_per_parallelcycle < 1
			error("Memory allocation is not enought (currently $MAX_MEM_PER_CPU [GB]). Please inclease MAX_MEM_PER_CPU or decrease number of stations")
		end

		if max_num_of_processes_per_parallelcycle >= length(event)

			S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), 1:length(event))

			# save data to jld2
			file = jldopen(fopath, "r+")
			unavalilablefile = jldopen(join([fopath[1:end-5], "_unavailablestations.jld2"]), "w+")

			for ii = 1:size(S)[1] #loop at each starttime
				varname = joinpath("event",InputDict["event"][ii]["origin"]["time"][1:end-4])
				file[varname] = S[ii]
			end

			JLD2.close(file)
			JLD2.close(unavalilablefile)

		else

			#parallelization by time
			pitr = 1

			while pitr <= length(event)

				startid1 = pitr
				startid2 = pitr + max_num_of_processes_per_parallelcycle - 1

				if startid1 == length(event)
					#use one
					S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), startid1:startid1)

				elseif startid2 <= length(event)
					#use all processors
					S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), startid1:startid2)

				else
					#use part of processors
					startid2 = startid1 + mod(length(event), nprocs()) - 1
					S = pmap(x -> seisdownload_EARTHQUAKE(x, InputDict), startid1:startid2)
				end

				# save data to jld2
				file = jldopen(fopath, "r+")
				unavalilablefile = jldopen(join([fopath[1:end-5], "_unavailablestations.jld2"]), "w+")

				for ii = 1:size(S)[1] #loop at each starttime
					varname = joinpath("event",InputDict["event"][startid1+ii-1]["origin"]["time"][1:end-4])
					file[varname] = S[ii]
				end

				JLD2.close(file)
				JLD2.close(unavalilablefile)

				pitr += max_num_of_processes_per_parallelcycle

				#println("pitr: $pitr")
			end
		end

    else
        println("Download type is not known (chose Noise or Earthquake).")
    end

    println("Downloading and Saving data is successfully done.\njob ended at "*string(now()))
    return 0

end

end


text:
using LogDensities
using Base.Test
using ForwardDiff

struct TestObj{p,T} <: parameters{T,1}
  x::Vector{T}
  prob::ProbabilityVector{p,T}
  pos::PositiveVector{p,T}
  real::RealVector{p,T}
  cov::CovarianceMatrix{p,T}
end

function inv_transform(::Type{TestObj{p,T2} where T2}, x::Vector{T}) where {T <: Real, p}
  Θ = construct(TestObj{p,T}, x)
  ConstrainedParameters.update_Σ!(Θ.cov)
  out = similar(x)
  k = 3p
  for i ∈ 1:p
    out[i] = Θ.prob[i]
    out[i+p] = Θ.pos[i]
    out[i+2p] = Θ.real[i]
    for j ∈ 1:i
      k += 1
      out[k] = Θ.cov.Σ[j,i]
    end
  end
  out
end

p = 2
m1 = Model(TestObj{p})
@testset begin
  @test type_length(typeof(m1.Θ)) == length(m1.Θ) == length(m1.Θ.x) == round(Int, p*(3 + (p+1)/2))

  x = randn(length(m1.Θ))
  copy!(m1.Θ.x, x)
  update!(m1.Θ)

  #The CovarianceMatrix log_jacobian! function drops the constant p*log(2) term.
  @test log_jacobian!(m1.Θ) + p*log(2) ≈ logabsdet(ForwardDiff.jacobian(x -> inv_transform(TestObj{2}, x), x))[1]
end

ext: jl


text:
f = open("input01.txt", "r")
    l = readlines(f)
close(f)

# Part 1
a = collect(zip(l, l[2:end]))
larger = 0
smaller = 0
for l in a
    i = parse(Int, l[1])
    j = parse(Int, l[2])
    if i < j
        global larger += 1
    elseif i > j
        global smaller += 1
end
end
println("Larger : $larger")
println("Smaller : $smaller")

# Part 2
s = 0
larger_2 = 0
smaller_2 = 0
b = collect(zip(l, l[2:end], l[3:end]))
for l in b
    ii = parse(Int, l[1])
    jj = parse(Int, l[2])
    kk = parse(Int, l[3])
    ss = ii + jj + kk
    if s == 0
        global s = ss
    elseif s < ss
        global larger_2 += 1
        global s = ss
    elseif s > ss
        global smaller_2 += 1
        global s = ss
    end
end
println("Larger_2 : $larger_2")
println("Smaller_2 : $smaller_2")

max_stars_repo_name: jung-benjamin/advent-of-code

max_issues_repo_name: jung-benjamin/advent-of-code

max_forks_repo_name: jung-benjamin/advent-of-code


text:
export Constraint, NoCons, L2Cons
export constrain!

@compat abstract type Constraint end

immutable NoCons <: Constraint
  threshold  :: AbstractFloat  # not used, just for consistent API
  every_n_iter :: Int          # also not used
end
NoCons() = NoCons(0.0, 0)

immutable L2Cons <: Constraint
  threshold    :: AbstractFloat
  every_n_iter :: Int
end
L2Cons(threshold) = L2Cons(threshold, 1)

############################################################
# No constraint
############################################################
function constrain!(backend::Backend, cons::NoCons, param)
  # do nothing if no constraints apply
end

############################################################
# L2 norm constraint on the weights
############################################################

function apply_l2_cons!{T <: AbstractFloat}(backend::CPUBackend, blob::CPUBlob{T},
                                            threshold::AbstractFloat, ninputs::Int, nunits::Int)
  param = reshape(blob.data, (ninputs, nunits))
  # we constrain each column vector
  for i = 1:nunits
    # compute norm and scale using blas
    norm = vecnorm(param[:, i])
    if norm > threshold
      scale_factor =  (1. / norm) * threshold
      offset = sizeof(T) * (i-1) * ninputs
      BLAS.scal!(ninputs, convert(T, scale_factor), pointer(param) + offset, 1)
    end
  end
end

# this constraints a given blob along the last dimension that is not of size 1
# it is a bit ugly but was the easiest way to implement it for now
function constrain!(backend :: Backend, cons :: L2Cons, param :: Blob)
  if ndims(param) == 1 # a bias blob
    apply_l2_cons!(backend, param, cons.threshold, length(param), 1)
  else # general weights
    apply_l2_cons!(backend, param, cons.threshold, get_fea_size(param), get_num(param))
  end
end


text:
requiredTypes = [ "AcademicLevel", "Affiliation", "Assignment", "Job", "Rank", "SubJob" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

export CareerStatus

type CareerStatus

    AcademicLevel::AbstractLevel
    Rank::AbstractLevel
    Assignment::AbstractLevel
    Affiliation::AbstractLevel
    Job::AbstractLevel
    SubJob::AbstractLevel

    MinStay::Int
    MaxStay::Int

    Attrition::Float16

    CareerStatus(;academicLevel::AbstractLevel = AcademicLevel(),
                  rank::AbstractLevel = Rank(),
                  assignment::AbstractLevel = Assignment(),
                  affiliation::AbstractLevel = Affiliation(),
                  job::AbstractLevel = Job(),
                  subjob::AbstractLevel = SubJob(),
                  min::Int = 0,
                  max::Int = 0,
                  attrition::Float16 = 0) = new(academicLevel, rank, assignment, affiliation, job, subjob, min, max, attrition)

    function CareerStatus(cs::CareerStatus)
        ncs = new()
        ncs.AcademicLevel = cs.AcademicLevel
        ncs.Rank = cs.Rank
        ncs.Assignment = cs.Assignment
        ncs.Affiliation = cs.Affiliation
        ncs.Job = cs.Job
        ncs.SubJob = cs.SubJob
        ncs.MinStay = cs.MinStay
        ncs.MaxStay = cs.MaxStay
        ncs.Attrition = cs.Attrition
        return ncs
    end

end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
using Random
using GymSpaces: Discrete
abstract type DiscreteEnv <: AbstractEnv end

function categorical_sample(prob_n, seed)
    csprob_n = cumsum(prob_n)
    return argmax(csprob_n .> rand(Float32))
end

mutable struct DiscreteEnvObj
    P
    isd
    lastaction
    nS
    nA

    action_space::Discrete
    observation_space::Discrete

    s
    seed::MersenneTwister
end

function DiscreteEnvObj(nS, nA, P, isd)
    action_space = Discrete(nA)
    observation_space = Discrete(nS)
    seed = MersenneTwister()

    s = categorical_sample(isd, seed)
    DiscreteEnvObj(P, isd, nothing, nS, nA, action_space, observation_space, s, seed)
end

# All DiscreteEnv must have a `discenv_obj` field...
function Base.getproperty(env::DiscreteEnv, sym::Symbol)
    if sym == :P
        return env.discenv_obj.P
    elseif sym == :isd
        return env.discenv_obj.isd
    elseif sym == :lastaction
        return env.discenv_obj.lastaction
    elseif sym == :nS
        return env.discenv_obj.nS
    elseif sym == :nA
        return env.discenv_obj.nA
    elseif sym == :action_space
        return env.discenv_obj.action_space
    elseif sym == :observation_space
        return env.discenv_obj.observation_space
    elseif sym == :s
        return env.discenv_obj.s
    elseif sym == :seed
        return env.discenv_obj.seed
    else
        return Base.getfield(env, sym)
    end
end

seed!(env::DiscreteEnv) = (env.seed = MersenneTwister())
seed!(env::DiscreteEnv, int) = (env.seed = MersenneTwister(int))

function reset!(env::DiscreteEnv)
    env.discenv_obj.s = categorical_sample(env.isd, env.seed)
    env.discenv_obj.lastaction = nothing
    return env.s
end

function step!(env::DiscreteEnv, action)
    @assert action ∈ env.action_space "Invalid action"
    transitions = env.discenv_obj.P[env.s][action]
    i = categorical_sample([t[1] for t ∈ transitions], env.seed)
    p, s, r, d = transitions[i]
    env.discenv_obj.s = s
    env.discenv_obj.lastaction = action
    return s, r, d, Dict(:prob => p)
end


text:
using Test
using Statistics

@test split_into_n(20, 3)  == [6, 7, 7]

num_of_polarimeters = 4
num_of_MPI_proc = 3
total_time = 50
baseline_length_s = 10
fsamp_hz = 10
fknee_hz = [0.01, 0.01]
σ_k = [0.1, 100]
rank = 1

baselines_per_process = [6, 7, 7]

@test split_tod_mpi(total_time, baseline_length_s, baselines_per_process, num_of_MPI_proc) ==  [[datachunk(1, 1, 5, 5), datachunk(2, 1, 1, 1)], [datachunk(2, 2, 5, 4), datachunk(3, 1, 3, 3)], [datachunk(3, 4, 5, 2), datachunk(4, 1, 5, 5)]]

chunks =  [[datachunk(1, 1, 5, 5), datachunk(2, 1, 1, 1)], [datachunk(2, 2, 5, 4), datachunk(3, 1, 3, 3)], [datachunk(3, 4, 5, 2), datachunk(4, 1, 5, 5)]]

@test get_chunk_properties(chunks, baseline_length_s, fsamp_hz, rank) ==  ([2,3], [10.0, 0.0], [49.901, 29.901], [4, 3], [400, 300])


text:
# Define plotting recipe
@recipe function f(s::SimpleSCM; kind = "overall")

    tp = s.treatment_panel

    legend --> :topleft

    if kind == "weights"

        wp = sort(DataFrame(comp = tp.is[Not(treated_ids(tp))],
                            weight = s.w),
                  :weight, rev = true)

        @series begin
            seriestype := :bar
            label --> ""
            xguide --> "Donor"
            yguide --> "Weight"
            linecolor := "white"
            wp[wp.weight .> 0.05, :comp], wp[wp.weight .> 0.05, :weight]
        end

    elseif kind == "overall"

        @series begin
            label --> treated_labels(tp)
            tp.ts, vec(tp.Y[treated_ids(tp), :])
        end

        @series begin
            label --> "Control"
            tp.ts, s.ŷ₁
        end

        @series begin
            label --> "Impact"
            seriestype --> :bar
            seriesalpha --> 0.5
            linecolor --> "white"
            tp.ts[only(length_T₀(tp)):end], s.τ̂
        end

        @series begin
            label --> "Intervention"
            seriestype := :vline
            linestyle := :dash
            seriescolor := "black"
            xguide := "Time"
            yguide := "Outcome"
            [tp.ts[only(length_T₀(tp))]]
        end

    elseif kind == "diffplot"
        @series begin
            label --> tp.treatment[1]
            sort!(unique(tp.data[!, tp.t_var])), [tp.y₁₀; tp.y₁₁] .- s.ŷ₁
        end

        @series begin
            label --> ""
            seriestype := :scatter
            seriescolor := 1
            markerstrokecolor := "white"
            sort!(unique(tp.data[!, tp.t_var])), [tp.y₁₀; tp.y₁₁] .- s.ŷ₁
        end

        @series begin
            label --> "Intervention"
            seriestype := :vline
            linestyle := :dash
            seriescolor := "black"
            xguide := "Time"
            yguide := "Outcome"
            [tp.treatment[2]]
        end
    end

     nothing
end

max_stars_repo_name: nilshg/SynthControl

max_issues_repo_name: nilshg/SynthControl

max_forks_repo_name: nilshg/SynthControl


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ThermalMultiStart <: ThermalGen
        name::String
        available::Bool
        status::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        rating::Float64
        primemover::PrimeMovers.PrimeMover
        fuel::ThermalFuels.ThermalFuel
        activepowerlimits::Min_Max
        reactivepowerlimits::Union{Nothing, Min_Max}
        ramplimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        power_trajectory::Union{Nothing, NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}}
        timelimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        start_time_limits::Union{Nothing, NamedTuple{(:hot, :warm, :cold), Tuple{Float64, Float64, Float64}}}
        start_types::Int
        op_cost::MultiStartCost
        basepower::Float64
        services::Vector{Service}
        time_at_status::Float64
        must_run::Bool
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data Structure for thermal generation technologies.

# Arguments
- `name::String`
- `available::Bool`
- `status::Bool`
- `bus::Bus`
- `activepower::Float64`, validation range: activepowerlimits, action if invalid: warn
- `reactivepower::Float64`, validation range: reactivepowerlimits, action if invalid: warn
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: (0, nothing), action if invalid: error
- `primemover::PrimeMovers.PrimeMover`: PrimeMover Technology according to EIA 923
- `fuel::ThermalFuels.ThermalFuel`: PrimeMover Fuel according to EIA 923
- `activepowerlimits::Min_Max`
- `reactivepowerlimits::Union{Nothing, Min_Max}`
- `ramplimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`, validation range: (0, nothing), action if invalid: error
- `power_trajectory::Union{Nothing, NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}}`: Power trajectory the unit will take during the startup and shutdown ramp process, validation range: (0, nothing), action if invalid: error
- `timelimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: Minimum up and Minimum down time limits in hours, validation range: (0, nothing), action if invalid: error
- `start_time_limits::Union{Nothing, NamedTuple{(:hot, :warm, :cold), Tuple{Float64, Float64, Float64}}}`:  Time limits for startup based on turbine temperature in hours
- `start_types::Int`:  Number of startup based on turbine temperature, validation range: (1, 3), action if invalid: error
- `op_cost::MultiStartCost`
- `basepower::Float64`: Base power of the unit in system base per unit, validation range: (0, nothing), action if invalid: warn
- `services::Vector{Service}`: Services that this device contributes to
- `time_at_status::Float64`
- `must_run::Bool`
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct ThermalMultiStart <: ThermalGen
    name::String
    available::Bool
    status::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "PrimeMover Technology according to EIA 923"
    primemover::PrimeMovers.PrimeMover
    "PrimeMover Fuel according to EIA 923"
    fuel::ThermalFuels.ThermalFuel
    activepowerlimits::Min_Max
    reactivepowerlimits::Union{Nothing, Min_Max}
    ramplimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Power trajectory the unit will take during the startup and shutdown ramp process"
    power_trajectory::Union{Nothing, NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}}
    "Minimum up and Minimum down time limits in hours"
    timelimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    " Time limits for startup based on turbine temperature in hours"
    start_time_limits::Union{Nothing, NamedTuple{(:hot, :warm, :cold), Tuple{Float64, Float64, Float64}}}
    " Number of startup based on turbine temperature"
    start_types::Int
    op_cost::MultiStartCost
    "Base power of the unit in system base per unit"
    basepower::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    time_at_status::Float64
    must_run::Bool
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function ThermalMultiStart(name, available, status, bus, activepower, reactivepower, rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, power_trajectory, timelimits, start_time_limits, start_types, op_cost, basepower, services=Device[], time_at_status=INFINITE_TIME, must_run=false, dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    ThermalMultiStart(name, available, status, bus, activepower, reactivepower, rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, power_trajectory, timelimits, start_time_limits, start_types, op_cost, basepower, services, time_at_status, must_run, dynamic_injector, ext, forecasts, InfrastructureSystemsInternal(), )
end

function ThermalMultiStart(; name, available, status, bus, activepower, reactivepower, rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, power_trajectory, timelimits, start_time_limits, start_types, op_cost, basepower, services=Device[], time_at_status=INFINITE_TIME, must_run=false, dynamic_injector=nothing, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    ThermalMultiStart(name, available, status, bus, activepower, reactivepower, rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, power_trajectory, timelimits, start_time_limits, start_types, op_cost, basepower, services, time_at_status, must_run, dynamic_injector, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function ThermalMultiStart(::Nothing)
    ThermalMultiStart(;
        name="init",
        available=false,
        status=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        rating=0.0,
        primemover=PrimeMovers.OT,
        fuel=ThermalFuels.OTHER,
        activepowerlimits=(min=0.0, max=0.0),
        reactivepowerlimits=nothing,
        ramplimits=nothing,
        power_trajectory=nothing,
        timelimits=nothing,
        start_time_limits=nothing,
        start_types=1,
        op_cost=MultiStartCost(nothing),
        basepower=0.0,
        services=Device[],
        time_at_status=INFINITE_TIME,
        must_run=false,
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::ThermalMultiStart) = value.name
"""Get ThermalMultiStart available."""
get_available(value::ThermalMultiStart) = value.available
"""Get ThermalMultiStart status."""
get_status(value::ThermalMultiStart) = value.status
"""Get ThermalMultiStart bus."""
get_bus(value::ThermalMultiStart) = value.bus
"""Get ThermalMultiStart activepower."""
get_activepower(value::ThermalMultiStart) = value.activepower
"""Get ThermalMultiStart reactivepower."""
get_reactivepower(value::ThermalMultiStart) = value.reactivepower
"""Get ThermalMultiStart rating."""
get_rating(value::ThermalMultiStart) = value.rating
"""Get ThermalMultiStart primemover."""
get_primemover(value::ThermalMultiStart) = value.primemover
"""Get ThermalMultiStart fuel."""
get_fuel(value::ThermalMultiStart) = value.fuel
"""Get ThermalMultiStart activepowerlimits."""
get_activepowerlimits(value::ThermalMultiStart) = value.activepowerlimits
"""Get ThermalMultiStart reactivepowerlimits."""
get_reactivepowerlimits(value::ThermalMultiStart) = value.reactivepowerlimits
"""Get ThermalMultiStart ramplimits."""
get_ramplimits(value::ThermalMultiStart) = value.ramplimits
"""Get ThermalMultiStart power_trajectory."""
get_power_trajectory(value::ThermalMultiStart) = value.power_trajectory
"""Get ThermalMultiStart timelimits."""
get_timelimits(value::ThermalMultiStart) = value.timelimits
"""Get ThermalMultiStart start_time_limits."""
get_start_time_limits(value::ThermalMultiStart) = value.start_time_limits
"""Get ThermalMultiStart start_types."""
get_start_types(value::ThermalMultiStart) = value.start_types
"""Get ThermalMultiStart op_cost."""
get_op_cost(value::ThermalMultiStart) = value.op_cost
"""Get ThermalMultiStart basepower."""
get_basepower(value::ThermalMultiStart) = value.basepower
"""Get ThermalMultiStart services."""
get_services(value::ThermalMultiStart) = value.services
"""Get ThermalMultiStart time_at_status."""
get_time_at_status(value::ThermalMultiStart) = value.time_at_status
"""Get ThermalMultiStart must_run."""
get_must_run(value::ThermalMultiStart) = value.must_run
"""Get ThermalMultiStart dynamic_injector."""
get_dynamic_injector(value::ThermalMultiStart) = value.dynamic_injector
"""Get ThermalMultiStart ext."""
get_ext(value::ThermalMultiStart) = value.ext

InfrastructureSystems.get_forecasts(value::ThermalMultiStart) = value.forecasts
"""Get ThermalMultiStart internal."""
get_internal(value::ThermalMultiStart) = value.internal

InfrastructureSystems.set_name!(value::ThermalMultiStart, val::String) = value.name = val
"""Set ThermalMultiStart available."""
set_available!(value::ThermalMultiStart, val::Bool) = value.available = val
"""Set ThermalMultiStart status."""
set_status!(value::ThermalMultiStart, val::Bool) = value.status = val
"""Set ThermalMultiStart bus."""
set_bus!(value::ThermalMultiStart, val::Bus) = value.bus = val
"""Set ThermalMultiStart activepower."""
set_activepower!(value::ThermalMultiStart, val::Float64) = value.activepower = val
"""Set ThermalMultiStart reactivepower."""
set_reactivepower!(value::ThermalMultiStart, val::Float64) = value.reactivepower = val
"""Set ThermalMultiStart rating."""
set_rating!(value::ThermalMultiStart, val::Float64) = value.rating = val
"""Set ThermalMultiStart primemover."""
set_primemover!(value::ThermalMultiStart, val::PrimeMovers.PrimeMover) = value.primemover = val
"""Set ThermalMultiStart fuel."""
set_fuel!(value::ThermalMultiStart, val::ThermalFuels.ThermalFuel) = value.fuel = val
"""Set ThermalMultiStart activepowerlimits."""
set_activepowerlimits!(value::ThermalMultiStart, val::Min_Max) = value.activepowerlimits = val
"""Set ThermalMultiStart reactivepowerlimits."""
set_reactivepowerlimits!(value::ThermalMultiStart, val::Union{Nothing, Min_Max}) = value.reactivepowerlimits = val
"""Set ThermalMultiStart ramplimits."""
set_ramplimits!(value::ThermalMultiStart, val::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}) = value.ramplimits = val
"""Set ThermalMultiStart power_trajectory."""
set_power_trajectory!(value::ThermalMultiStart, val::Union{Nothing, NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}}) = value.power_trajectory = val
"""Set ThermalMultiStart timelimits."""
set_timelimits!(value::ThermalMultiStart, val::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}) = value.timelimits = val
"""Set ThermalMultiStart start_time_limits."""
set_start_time_limits!(value::ThermalMultiStart, val::Union{Nothing, NamedTuple{(:hot, :warm, :cold), Tuple{Float64, Float64, Float64}}}) = value.start_time_limits = val
"""Set ThermalMultiStart start_types."""
set_start_types!(value::ThermalMultiStart, val::Int) = value.start_types = val
"""Set ThermalMultiStart op_cost."""
set_op_cost!(value::ThermalMultiStart, val::MultiStartCost) = value.op_cost = val
"""Set ThermalMultiStart basepower."""
set_basepower!(value::ThermalMultiStart, val::Float64) = value.basepower = val
"""Set ThermalMultiStart services."""
set_services!(value::ThermalMultiStart, val::Vector{Service}) = value.services = val
"""Set ThermalMultiStart time_at_status."""
set_time_at_status!(value::ThermalMultiStart, val::Float64) = value.time_at_status = val
"""Set ThermalMultiStart must_run."""
set_must_run!(value::ThermalMultiStart, val::Bool) = value.must_run = val
"""Set ThermalMultiStart ext."""
set_ext!(value::ThermalMultiStart, val::Dict{String, Any}) = value.ext = val

InfrastructureSystems.set_forecasts!(value::ThermalMultiStart, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set ThermalMultiStart internal."""
set_internal!(value::ThermalMultiStart, val::InfrastructureSystemsInternal) = value.internal = val


text:
# this is the funtion that appears in the
# FoggySurface and FoggyContour graphics,
# It is the Octave sombrero, tilted a bit
function FoggyMountainObj(θ₁, θ₂)
   	r = sqrt(θ₁^2.0 + θ₂^2.0) + 1e-10
	z = sin(r) / (2*r)
	z = z + θ₁/40.0 - 2*(θ₁/40.0)^2.0 + θ₂/40.0 - 2*(θ₂/40.0)^2.0
	z = -z/10.0 # switch to minimization
end

function FoggyMountainObj(θ)
	θ₁ = θ[1]
	θ₂ = θ[2]
    FoggyMountainObj(θ₁, θ₂)
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module CodeGen

import Base

# using Poptart.Desktop
using ExprTools
using DataStructures
using Comonicon
using Comonicon.Types
using Comonicon.Types: CommandDoc
using Comonicon.CodeGen
using Comonicon.CodeGen: hasparameters

include("poptart/expr.jl")
using .PoptartExpr
include("poptart/context.jl")

export codegen, PoptartCtx

end


text:
"""
```
find_threshold(Entropy(),counts,edges)
```
An algorithm for finding the threshold value for a gray-level histogram using
the entropy of the histogram.

# Output

Returns the point in the `AbstractRange` which corresponds to the threshold bin in the histogram.

# Details

This algorithm uses the entropy of a gray level histogram to produce a threshold
value.

Let ``f_1, f_2, \\ldots, f_I`` be the frequencies in the various bins of the
histogram and ``I`` the number of bins. With ``N = \\sum_{i=1}^{I}f_i``, let
``p_i = \\frac{f_i}{N}`` (``i = 1, \\ldots, I``) denote the probability
distribution of gray levels. From this distribution one derives two additional
distributions. The first defined for discrete values ``1`` to ``s`` and the
other, from ``s+1`` to ``I``. These distributions are

```math
A: \\frac{p_1}{P_s}, \\frac{p_2}{P_s}, \\ldots, \\frac{p_s}{P_s}
\\quad \\text{and} \\quad
B: \\frac{p_{s+1}}{1-P_s}, \\ldots, \\frac{p_n}{1-P_s}
\\quad \\text{where} \\quad
P_s = \\sum_{i=1}^{s}p_i.
```
The entropies associated with each distribution are as follows:

```math
H(A) = \\ln(P_s) + \\frac{H_s}{P_s}
```
```math
H(B) = \\ln(1-P_s) + \\frac{H_n-H_s}{1-P_s}
```
```math
\\quad \\text{where} \\quad
H_s = -\\sum_{i=1}^{s}p_i\\ln{p_i}
\\quad \\text{and} \\quad
H_n = -\\sum_{i=1}^{I}p_i\\ln{p_i}.
```

Combining these two entropy functions we have

```math
\\psi(s) = \\ln(P_s(1-P_s)) + \\frac{H_s}{P_s} + \\frac{H_n-H_s}{1-P_s}.
```
Finding the discrete value ``s`` which maximises the function ``\\psi(s)`` produces
the sought-after threshold value (i.e. the bin which determines the threshold).

See Section 4 of [1] for more details on the derivation of the entropy.

# Options

## Choices for `counts`
You can specify an `AbstractArray` which should be a 1D array of frequencies for
a histogram. You should submit the corresponding `edges` range for the bins of
the histogram. The function will throw an error if it detects that the `edges`
and `counts` have different lengths.

## Choices for `edges`
You can specify an `AbstractRange` which should be the corresponding range for
the bins of the histogram array passed into `counts`.

# Example

```julia

using TestImages, Images

img = testimage("cameraman")
# building a histogram with 256 bins
edges, counts = build_histogram(img, 256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
find_threshold(EntropyThresholding(), counts[1:end], edges)
```

# References
[1] J. N. Kapur, P. K. Sahoo, and A. K. C. Wong, “A new method for gray-level picture thresholding using the entropy of the histogram,” *Computer Vision, Graphics, and Image Processing*, vol. 29, no. 1, p. 140, Jan. 1985.[doi:10.1016/s0734-189x(85)90156-2](https://doi.org/10.1016/s0734-189x%2885%2990156-2)
"""
function find_threshold(algorithm::Entropy, counts::AbstractArray, edges::AbstractRange)
    if length(edges) != length(counts)
        error("the lengths of edges and counts must match")
    end
    nbins = length(counts)
    pdf = counts ./ sum(counts)
    Ψₘₐₓ = 0
    Hₛ = H(pdf)
    Hₙ = -1*Hₛ[nbins]
    Pₛ = cumsum(pdf)
    threshold_bin = 1
    for i in 1:nbins
        Pᵢ = Pₛ[i]
        Hᵢ = -1*Hₛ[i]
        Ψᵢ = log(Pᵢ*(1-Pᵢ)) + (Hᵢ/Pᵢ) + ((Hₙ-Hᵢ)/(1-Pᵢ))
        if Ψᵢ > Ψₘₐₓ
            Ψₘₐₓ = Ψᵢ
            threshold_bin = i
        end
    end
    edges[threshold_bin]
end

function H(pdf::AbstractArray)
    h_dist = zeros(length(pdf))
    for i in 1:length(pdf)
        h_dist[i] = (pdf[i] != 0) ? pdf[i]*log(pdf[i]) : 0.0
    end
    cumsum(h_dist)
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 IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration <: SwaggerModel
    assuredConcurrencyShares::Any # spec type: Union{ Nothing, Int32 } # spec name: assuredConcurrencyShares
    limitResponse::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1LimitResponse } # spec name: limitResponse

    function IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration(;assuredConcurrencyShares=nothing, limitResponse=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration, Symbol("assuredConcurrencyShares"), assuredConcurrencyShares)
        setfield!(o, Symbol("assuredConcurrencyShares"), assuredConcurrencyShares)
        validate_property(IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration, Symbol("limitResponse"), limitResponse)
        setfield!(o, Symbol("limitResponse"), limitResponse)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration

const _property_map_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration = Dict{Symbol,Symbol}(Symbol("assuredConcurrencyShares")=>Symbol("assuredConcurrencyShares"), Symbol("limitResponse")=>Symbol("limitResponse"))
const _property_types_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration = Dict{Symbol,String}(Symbol("assuredConcurrencyShares")=>"Int32", Symbol("limitResponse")=>"IoK8sApiFlowcontrolV1alpha1LimitResponse")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1LimitedPriorityLevelConfiguration)
    true
end

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


text:
__precompile__()
module Nettle
using Compat
import Compat.String

# Load libnettle from our deps.jl
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Nettle not installed properly, run Pkg.build(\"Nettle\"), restart Julia and try again")
end
include(depsjl_path)

include("hash_common.jl")
include("hash.jl")
include("hmac.jl")
include("cipher.jl")

function __init__()
    # Always check your dependencies that live in `deps.jl`
    check_deps()
end

# SnoopCompile acceleration
include("precompile.jl")
_precompile_()

end # module


text:
using Requests
using HttpCommon
using HttpParser

function process_response(stream)
    r = Response()
    rp = Requests.ResponseParser(r,stream)
    while isopen(stream)
        data = readavailable(stream)
        if length(data) > 0
            http_parser_execute(rp.parser, rp.settings, data)
        end
    end
    http_parser_execute(rp.parser,rp.settings,"") #EOF
    r
end

clientside = connect("/tmp/julia.socket")
req = Requests.default_request("GET", "/", "/tmp/julia.socket", "")
dump(req)
write(clientside, Requests.render(req))
dump(process_response(clientside))


text:
using SecretSanta
using Test

@testset "SecretSanta" begin
    @testset "SecretSantaModel" begin
        SecretSanta.run("../test/test.json", test = true)
        @test true
    end
end


text:
#=

This file contains helper functions for creating calculation models

`make_hitran_model` creates a HitranModel with all specified parameters.

`make_interpolation_model` creates a cross-section matrix based on the pressure and
temperature grids, and fits an interpolation.

There are also some save/load convenience functions for the interpolation model.

=#

"""
    make_hitran_model(hitran::HitranTable,
                      broadening::AbstractBroadeningFunction;
                      wing_cutoff::Real=40,
                      vmr::Real=0,
                      CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                      architecture = default_architecture)

Convenience function to make a HitranModel out of the parameters (Matches make_interpolation_model)

"""
function make_hitran_model(hitran::HitranTable,
                           broadening::AbstractBroadeningFunction;
                           wing_cutoff::Integer=40,
                           vmr::Union{Real, Vector}=0,
                           CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                           architecture = default_architecture)

    if architecture isa GPU && !(CEF isa HumlicekWeidemann32SDErrorFunction)
        @warn "Cross-section calculations on GPU may or may not work with this CEF (use HumlicekWeidemann32SDErrorFunction if you encounter issues)"
    end

    return HitranModel(hitran=hitran, broadening=broadening , wing_cutoff=wing_cutoff , vmr=vmr, CEF=CEF, architecture=architecture)
end

"""
    $(FUNCTIONNAME)(hitran::HitranTable,
                    broadening::AbstractBroadeningFunction,
                    wave_grid::AbstractRange{<:Real},
                    p_grid::AbstractRange{<:Real},
                    t_grid::AbstractRange{<:Real};
                    wavelength_flag::Bool=false,
                    wing_cutoff::Real=40,
                    vmr::Real=0,
                    CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                    architecture::AbstractArchitecture=default_architecture )

Using a HitranModel, create an InterpolationModel by interpolating the lineshape function
at the given pressure and temperature grids

"""
function make_interpolation_model(
                                  # Required
                                  hitran::HitranTable,
                                  broadening::AbstractBroadeningFunction,
                                  wave_grid::AbstractRange{<:Real},
                                  p_grid::AbstractRange{<:Real},
                                  t_grid::AbstractRange{<:Real};
                                  # Optionals
                                  wavelength_flag::Bool=false,
                                  wing_cutoff::Integer=40,
                                  vmr::Real=0,
                                  CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                                  architecture::AbstractArchitecture=default_architecture     # Computer `Architecture` on which `Model` is run
                                )

    # Warn user if using incompatible/untested CEF
    if architecture isa Architectures.GPU && !(CEF isa HumlicekWeidemann32SDErrorFunction)
        @warn "Cross-section calculations on GPU may or may not work with this CEF (use HumlicekWeidemann32SDErrorFunction if you encounter issues)"
    end

    # Convert from wavelength to wavenumber if necessary
    ν_grid = wavelength_flag ? reverse(nm_per_m ./ wave_grid) : wave_grid

    # Empty matrix to store the calculated cross-sections
    cs_matrix = zeros(length(p_grid), length(t_grid), length(ν_grid));

    # Calculate all the cross-sections at the pressure and temperature grids
    @showprogress 1 "Computing Cross Sections for Interpolation..." for i in 1:length(p_grid)
        for j in 1:length(t_grid)
            # make_hitran_model(hitran_data, Voigt(), wing_cutoff = 40, CEF=HumlicekWeidemann32SDErrorFunction(), architecture=CPU())
            model = make_hitran_model(hitran, broadening, wing_cutoff=wing_cutoff, CEF=CEF, architecture=architecture)
            cs_matrix[i,j,:] = Array(compute_absorption_cross_section(model, collect(ν_grid), p_grid[i], t_grid[j], wavelength_flag=wavelength_flag))
        end
    end

    # Perform the interpolation
    itp = interpolate(cs_matrix, BSpline(Cubic(Line(OnGrid()))))

    # Get the molecule and isotope numbers from the HitranTable
    mol = hitran.mol[1]
    iso = all(x->x==hitran.iso[1], hitran.iso) ? hitran.iso[1] : -1

    # Return the interpolation model with all the proper parameters
    return InterpolationModel(itp, mol, iso, broadening, ν_grid, p_grid, t_grid, wing_cutoff, vmr, CEF, architecture)
end

""" Convenience function to save an InterpolationModel at a specified filepath (Using JLD2) """
function save_interpolation_model(itp_model::InterpolationModel, filepath::String)
    @save filepath itp_model
end

""" Convenience function to load an InterpolationModel from a specified filepath (Using JLD2) """
function load_interpolation_model(filepath::String)
    itp_model::InterpolationModel = @load filepath itp_model
    return itp_model
end


text:
"""
a, b, c  = stripboundedge(Mesh)

returns a mesh with the inner edges: a
and returns a mesh with the boundary edges: b
and returns a mesh of the faces on the boundary: c
"""

function stripboundedge(Γ::Mesh{3,4,T}) where T
    f = skeleton(Γ,2)
    C = connectivity(f,Γ)

    totfaces = length(f.faces)
    v = SArray{Tuple{3},T,1,3}
    A = zeros(v, 2*totfaces-4*length(Γ.faces))

    l = 1
    for (i,k) in enumerate(f.faces)
        if length(nzrange(C,i)) ≈ 1
            A[l] = k
            l += 1
        end
    end

    #i = 1
    #k = 1
    #while i < totfaces+1
    #    if length(BEAST.nzrange(C,i)) ≈ 1
    #        A[k] = f.faces[i]
    #        k += 1
    #    end
    #    i += 1
    #end

    edges = skeleton(Γ,1)

    Γ2 = Mesh(Γ.vertices,A)
    rem = skeleton(Γ2,1)
    savedges = setdiff(edges.faces,rem.faces)
    B = Mesh(Γ.vertices,savedges)
    return submesh(B,skeleton(Γ,1)), submesh(rem,skeleton(Γ,1)), submesh(Γ2,skeleton(Γ,2))
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
######################################################################
# weierstrasspoints.jl: Addition laws for projective points on Weierstrass curves
######################################################################

export infinity, projective_add, projective_scalar_mul

######################################################################
# Basic methods
######################################################################

"""
Get the point at infinity on an elliptic curve in Weierstrass form.
"""
function infinity(E::AbstractWeierstrass)
    R = base_ring(E)
    return EllipticPoint(Nemo.zero(R), Nemo.one(R), Nemo.zero(R), E)
end

function zero(E::AbstractWeierstrass)
	return infinity(E)
end

######################################################################
# Addition law
######################################################################

"""
Get the opposite of a point on an elliptic curve in Weierstrass form.
"""
function -(P::EllipticPoint)
	E = P.curve
	x, y, z = coordinates(P)
	a1, _, a3, _, _ = a_invariants(E)
    return EllipticPoint(x, - y - a1 * x - a3 * z, z, E)
end

"""
Get the sum of two normalized projective points on the same Weierstrass curve, assuming they are not equal and not inverse of each other.
"""
function _addgeneric(P::EllipticPoint{T}, Q::EllipticPoint{T}) where T<:Nemo.FieldElem
    E = P.curve
	a1, a2, a3, _, _ = a_invariants(E)

	xq, yq, zq = coordinates(Q)
	@assert zq == 1
	xp, yp, zp = coordinates(P)
	@assert zp == 1

	#sanity check
	@assert xp != xq

    denom = xq - xp
	inverse = 1//denom
    lambda = (yq - yp) * inverse
    nu = (yp * xq - xp * yq) * inverse

    Xplus = lambda^2 + a1 * lambda - a2 - xp - xq
    Yplus = -(lambda + a1) * Xplus - nu - a3
    Zplus = Nemo.one(base_ring(P))
    return EllipticPoint(Xplus, Yplus, Zplus, E)
end

"""
Get the sum of two normalized projective points on the same Weierstrass curve, assuming they have equal x-coordinate.
"""
function _addequalx(P::EllipticPoint{T}, Q::EllipticPoint{T}) where T<:Nemo.FieldElem
    E = P.curve
	a1, a2, a3, a4, a6 = a_invariants(E)

	xp, yp, zp = coordinates(P)
	@assert zp == 1
	xq, yq, zq = coordinates(Q)
	@assert zq == 1

	#sanity check
	@assert xp == xq

    denom = yp + yq + a1 * xq + a3
    if iszero(denom)
	    return infinity(E)
    else
		inverse = 1//denom
	    lambda = (3 * xp^2 + 2 * a2 * xp + a4 - a1 * yp) * inverse
	    nu = (- xp^3 + a4 * xp + 2 * a6 - a3 * yp) * inverse
	    Xplus = lambda^2 + a1 * lambda - a2 - xp - xq
        Yplus = -(lambda + a1) * Xplus - nu - a3
        Zplus = one(base_ring(P))
	    return EllipticPoint(Xplus, Yplus, Zplus, E)
    end
end

"""
Get the double of a normalized point.
"""
function _double(P::EllipticPoint)
	if isinfinity(P)
		return infinity(P.curve)
	else
		return _addequalx(P, P)
	end
end

"""
Get the sum of two projective points on the same Weierstrass curve.

"""
function +(P::EllipticPoint{T}, Q::EllipticPoint{T}) where T<:Nemo.FieldElem
    P = normalized(P)
    Q = normalized(Q)
	xp, _, _ = coordinates(P)
	xq, _, _ = coordinates(Q)
    if isinfinity(P)
        return Q
    elseif isinfinity(Q)
        return P
    elseif xp == xq
		return _addequalx(P,Q)
    else
		return _addgeneric(P,Q)
    end
end

function -(P::EllipticPoint{T}, Q::EllipticPoint{T}) where T<:FieldElem
	return P + (-Q)
end

"""
Get a scalar multiple of a point on a Weierstrass curve.
"""

function *(k::Integer, P::EllipticPoint)
	E = P.curve
	P = normalized(P)
	if k == 0
		return infinity(E)
	elseif k<0
		return (-k) * (-P)
	else
		if isinfinity(P)
			return infinity(E)
		else
			return _ladder(k, P)
		end
	end
end

# function *(k::Integer, P::EllipticPoint)
# 	return k * P
# end

function _ladder(m::Integer, P::EllipticPoint)
    p0 = P
    for b in Iterators.drop(Iterators.reverse(digits(Int8, m, base=2)), 1)
        if (b == 0)
            p0 = _double(p0)
        else
            p0 = P + _double(p0)
        end
    end
    return p0
end

######################################################################
# Projective addition law for Short Weierstrass curves
######################################################################

# P = (x1, y1, z1)
# Q = (x2, y2, z2)
# P + Q = (x3, y3, z3)
# if P != Q then
# x3 = (x2 z1 - x1 z2) [ (y2 z1 - y1 z2)^2 z1 z2 - (x2 z1 - x1 z2)^2 (x2 z1 + x1 z2) ]
# y3 = (y2 z1 - y1 z2) [ (x2 z1 - x1 z2)^2 (x2 z1 + 2 x1 z2) - (y2 z1 - y1 z2)^2 z1 z2 ] - (x2 z1 - x1 z2)^3 y1 z2
# z3 = (x2 z1 - x1 z2)^3 z1 z2
# if P = Q then
# x3 = 2 y1 z1 [ (a z1^2 + 3 x1^2)^2 - 8 x1 y1^2 z1 ]
# y3 = (a z1^2 + 3 x1^2 ) [ 12 x1 y1^2 z1 - a^2 (z1^2 + 3 x1^2)^2 ] - 8 y1^4 z1^2
# z3 = (2 y1 z1)^3

#This function is only to be used with distinct points on the same short Weierstrass curve
#xdet is x2 z1 - x1 z2, and ydet is y2 z1 - y1 z2
function _projective_add_neq(P::EllipticPoint{T}, Q::EllipticPoint{T}, xdet::T, ydet::T) where T
	x1, y1, z1 = coordinates(P)
	x2, y2, z2 = coordinates(Q)
	xdet2 = xdet^2
	xdet3 = xdet2 * xdet
	ydet2 = ydet^2
	z1z2 = z1 * z2
	#we could save a few multiplications in what follows
	x3 = xdet * ( ydet2 * z1z2 - xdet2 * (x2 * z1 + x1 * z2) )
	y3 = ydet * ( xdet2 * (x2 * z1 + 2 * x1 * z2) - ydet2 * z1z2 ) - xdet3 * y1 * z2
	z3 = xdet3 * z1z2
	res = Point(x3, y3, z3, base_curve(P))
	#@assert isvalid(res) #sanity check
	return res
end

#This function is only to be used with a point on a short Weierstrass curve
function _projective_dbl(P::EllipticPoint{T}) where T
	x, y, z = coordinates(P)
	_, _, _, a, b = a_invariants(base_curve(P))
	factor = a * z^2 + 3 * x^2
	factor2 = factor^2
	yz = y * z
	xy2z = x * y * yz
	#we could save again a few multiplications in what follows
	xprime = 2 * yz * ( factor2 - 8 * xy2z )
	yprime = factor * ( 12 * xy2z - factor2 ) - 8 * (yz)^2 * y^2
	zprime = 8 * (yz)^3
	res = Point(xprime, yprime, zprime, base_curve(P))
	#@assert isvalid(res) #sanity check
	return res
end

#This function is only to be used with two points on the same short Weierstrass curve
#Compute xdet and ydet, if they are both zero go to _projective_dbl
function projective_add(P::EllipticPoint{T}, Q::EllipticPoint{T}) where T
	x1, y1, z1 = coordinates(P)
	x2, y2, z2 = coordinates(Q)
	xdet = x2 * z1 - x1 * z2
	ydet = y2 * z1 - y1 * z2
	if ((xdet == 0) & (ydet == 0))
		return _projective_dbl(P)
	else
		return _projective_add_neq(P, Q, xdet, ydet)
	end
end

#Here we assume v is a positive integer and P lives on a short Weierstrass curve
function _projective_scalar_mul(P::EllipticPoint, v::Integer)
    P0 = P
    for b in Iterators.drop(Iterators.reverse(digits(Int8, v, base=2)), 1)
        P0 = _projective_dbl(P0)
        if (b == 1)
            P0 = projective_add(P0, P)
        end
    end
    return P0
end

function projective_scalar_mul(P::EllipticPoint, v::Integer)
	if v == 0
		return infinity(base_curve(E))
	elseif v < 0
		return - _projective_scalar_mul(P, -v)
	else
		return _projective_scalar_mul(P, v)
	end
end


text:
# This file imports all of the mlpack functions into the mlpack module.

adaboost = _Internal.adaboost
approx_kfn = _Internal.approx_kfn
bayesian_linear_regression = _Internal.bayesian_linear_regression
cf = _Internal.cf
dbscan = _Internal.dbscan
decision_stump = _Internal.decision_stump
decision_tree = _Internal.decision_tree
det = _Internal.det
emst = _Internal.emst
fastmks = _Internal.fastmks
gmm_train = _Internal.gmm_train
gmm_generate = _Internal.gmm_generate
gmm_probability = _Internal.gmm_probability
hmm_train = _Internal.hmm_train
hmm_loglik = _Internal.hmm_loglik
hmm_viterbi = _Internal.hmm_viterbi
hmm_generate = _Internal.hmm_generate
hoeffding_tree = _Internal.hoeffding_tree
kde = _Internal.kde
kernel_pca = _Internal.kernel_pca
kmeans = _Internal.kmeans
lars = _Internal.lars
linear_regression = _Internal.linear_regression
linear_svm = _Internal.linear_svm
lmnn = _Internal.lmnn
local_coordinate_coding = _Internal.local_coordinate_coding
logistic_regression = _Internal.logistic_regression
lsh = _Internal.lsh
mean_shift = _Internal.mean_shift
nbc = _Internal.nbc
nca = _Internal.nca
knn = _Internal.knn
kfn = _Internal.kfn
nmf = _Internal.nmf
pca = _Internal.pca
perceptron = _Internal.perceptron
preprocess_split = _Internal.preprocess_split
preprocess_binarize = _Internal.preprocess_binarize
preprocess_describe = _Internal.preprocess_describe
preprocess_scale = _Internal.preprocess_scale
preprocess_one_hot_encoding = _Internal.preprocess_one_hot_encoding
image_converter = _Internal.image_converter
radical = _Internal.radical
random_forest = _Internal.random_forest
krann = _Internal.krann
softmax_regression = _Internal.softmax_regression
sparse_coding = _Internal.sparse_coding


text:
########### Defines the parameters & plot for Fig 1 a,c,e,g,i ###########
bi = 0.1; # set the birth rate of infected individuals as a constant 0.1
# import required variables/modules common to all simulations
include("../Model_scripts/fig1params.jl")
include("../Model_scripts/Functions.jl")
# now that we have all the requirements, let's run the simulation
run_simulation();
Np = N';
########### Figure 1a ###########
# plot the number of infected/uninfected hosts
labels = vcat("Infected", ["" for i in 1:1:n_parasites]);
plot_population_numbers(
    hcat(labels...),
    "by=0.1\n \nNumber of infected and uninfected hosts",
    100,
    size(Np)[1],
    "$figure_directory/graph_1a.png"
)
########### Figure 1g ###########
# calculate the weighted average mortality
ui_w_avg = calculate_mortality()
# plot the weighted average mortality rate through time
plot_mortality(1, size(Np)[1], "$figure_directory/graph_1g.png")
########### Figure 1e ###########
# calculate the weighted average β
βy_w_avg = calculate_horizontal_transmission()
# calculate weighted H0
H0_w = calculate_H0()
# calculate weighted V0
V0_w = calculate_V0()
# plot the average weighted V0
plot_vertical_reproductive_ratio(1, size(Np)[1], "$figure_directory/graph_1e.png")
########### Figure 1c ###########
# plot the average weighted R0
plot_reproductive_rate(6, size(Np)[1], "$figure_directory/graph_1c.png")
########### Figure 1i ###########
# calculate & plot the population evenness through time
plot_evenness(1, size(Np)[1], "$figure_directory/graph_1i.png")


text:
module Diderot

using DataStructures

## Core

include("types.jl")
include("interface.jl")
include("implementation.jl")

## Generic implementation

include("restriction.jl")
include("relaxation.jl")

## Specific implementation for problem classes.

include("knapsack.jl")
include("setcover.jl")
include("indexfund.jl")

# add README as docstring to module
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end Diderot

end


text:
#This function implements the algorithm from  https://link.springer.com/content/pdf/10.1007/BF02289694.pdf
# Varibales are generally named simillaryly to thos in the original paper. However
# it felt simples to have a a vector of Is and one of Js, rather than an IJ array.
# All functions used can be found in functions.jl

# The function takes a distance matrix D and returns a K,L array of transformed coordinates
# The algorithm exits once the magnitude of the gradient is less than the maginitude
# from a random arrangement of coordinates multiplied by `local_minimum_criterion`
# or once it has gone through `max_iter` iterations.

function nmds(D,L ::Int64,local_minimum_criterion = 0.02, max_iter ::Int64 = 10000)
    # initialise varables and containers
    local_minimum = false
    S=1
    K = size(D)[1] # Number of observations/points
    X2 = zeros(K,L) # Current gradient
    X3 = zeros(K,L) # Previous gradient
    α = 0.2 # Initial step size
    rand_mag = 0.0 # To be updated as random magnitude of gradient
    stress_vec = CircularBuffer(6) # Current and previous 5 stress calculations
    DISSIM = Vector{Float64}(undef,0) # Original dissimilarities
    Is = Vector{Int64}(undef,0) # i indeces
    Js = Vector{Int64}(undef,0) # j indeces
    @inbounds for i in 2:K
        for j in 1:i-1
            push!(DISSIM,D[i,j])
            push!(Is,i)
            push!(Js,j)
        end
    end
    M = length(Is)
    DIST = Vector{Float64}(undef,M) # euclidean distances between points in X1
    DHAT = Vector{Float64}(undef,M) # fitted distances
    # gradient calculations will be broken into blocks to give vectors of length M
    g2 = Vector{Float64}(undef,M) # solutions to 2nd term on rhs of gradient calculation
    g3 = Array{Float64}(undef,M,L) # solutions to 3rd term on rhs of gradient calculation
    grhs = Vector{Float64}(undef,M) # combined 2nd and third terms for corresponding value of l
                                    # this will be repeatedly overwritten as we iterate over l
    # sort parameters by order of values in DISSIM
    order = sortperm(DISSIM)
    DISSIM = DISSIM[order]
    Is = Is[order]
    Js = Js[order]
    # generate random starting values for X1
    X1 = normalize(randn(K,L),1) # current point values
    iter = 0
    while (local_minimum == false) & (iter <max_iter)
        iter += 1
        get_dist!(X1,DIST,Is,Js,M,L) # calculate euclidean distances between points
        monotone!(DIST,DHAT,M) # fit DHAT values vie monotone regression
        S_star, T_star, S = stress(DIST,DHAT,M) # calculate stress
        push!(stress_vec,S)   # accumulate stress in vector
        grhs2!(DIST, DHAT, S_star,T_star,M,g2) # calculate second term of gradient rhs
        grhs3!(X1,DIST, M,L,g3,Is,Js) # calculate 3rd ter of gradient rhs
        X2 = g(g2,g3,L,M,S,Is,Js,X2,grhs) # put them together and calculate negative gradient
        # if there are 2 or more stress calculations so far, use Kruskals method to get new step size
        α = length(stress_vec) < 2 ? α : α * angle_factor(X2,X3) * relaxation_factor(stress_vec) * good_luck_factor(stress_vec)
        #X1 = normalize(x′(X1,X2,α,K),1) # get new point values
        #X1 = x′(X1,X2,α,K)
        x′!(X1,X2,α,K,L)
        normalize!(X1,K,L)
        X3 .= X2 # update previous gradient

        if iter == 1 # on first iteration
            rand_mag = mag(X2,K) # magnitude of gradient
        elseif  mag(X2,K) < local_minimum_criterion*rand_mag # if local minimum
            local_minimum = true # exit loop
        end
        #println(S)
    end
    return X1, S
end

# run the above algorith n times and return the result with the lowest stress value
function nmds(D,L ::Int64,n ::Int64, local_minimum_criterion = 0.02, max_iter ::Int64 = 10000)
    results = []
    Threads.@threads for i in 1:n
        push!(results,nmds(D,L,local_minimum_criterion,max_iter))
    end
    best = findmin([results[i][2] for i in 1:n])[2]
    return results[best]
end


text:
module AdventOfCode2021

include("day01/day01.jl")
include("day02/day02.jl")
include("day03/day03.jl")

end


text:
#=11)ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля,
на котором могут находиться также внутренние прямоугольные перегородки
(все перегородки изолированы друг от друга, прямоугольники могут вырождаться в отрезки)
РЕЗУЛЬТАТ: Робот - в исходном положении, и в 4-х приграничных клетках, две из которых имеют ту же широту, а две - ту же долготу,
что и Робот, стоят маркеры.=#
#=
Примерный "движения" для решения задачи:
1)спускаться в Ю.З угол (аналогично с задачей про маркеры в углах) но сохранять кол-во шагов в массивах для West и Sud
2.1)Проходить вдоль внешней перегородки известное кол-во шагов и ставить маркеры
2.2)по известным кол-вам шагов из Ю.З можно пойти Ost,Nord потом вернуться обратно к углу и поменять направления
3)вернуться используя массив
=#

function ends_of_cross(r::Robot)
    num_steps=[] #создаём массив для num_steps
    to_West=0
    to_Sud=0
    while isborder(r,Sud)==false || isborder(r,West)==false
        push!(num_steps, get_num_movements!(r, West))
        push!(num_steps, get_num_movements!(r, Sud))
    end
    #робот в ЮЗ углу;из массива в котором индексы и кол-во шагов
    #суммируем шаги проверяя индексы(нечётные-west,чётные-sud)
    for (i,n) in enumerate(num_steps)
        if isodd(i)
            to_West+=n
        else
            to_Sud+=n
        end
    end
   marker_in_line(r,Nord,to_Sud)
   marker_in_line(r,Ost,to_West)
   moves!(r,Sud)
   moves!(r,West)
   marker_in_line(r,Ost,to_West)
   marker_in_line(r,Nord,to_Sud)
   moves!(r,West)
   moves!(r,Sud)

    for (i,n) in enumerate(num_steps)
        side = isodd(i) ? Ost : Nord
        moves!(r,side,n)
    end
end
function marker_in_line(r::Robot,side::HorizonSide,num_steps::Int)#ставим один маркер(num_steps-параметр исходной точки)
    moves!(r,side,num_steps)
    putmarker!(r)
    moves!(r,side)
end
function moves!(r::Robot, side::HorizonSide, num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end
function moves!(r::Robot, side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
    end
end
function get_num_movements!(r::Robot, side::HorizonSide)
    num_steps = 0
    while isborder(r,side)==false
        move!(r,side)
        num_steps += 1
    end
    return num_steps
end

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
using ModelTES, Roots#, Plots
using PyPlot
#pyplot() #make Plots use pyplot backend
using PyCall
@pyimport numpy as np

function highEpixHighBeta()
    Tc=0.107; Tbath = 0.065
    R0 = 1.55e-3; Rl = 0.35e-3; Rn = 10.3e-3; Rpara=0.0
    n=3.3; k = 3.50e-9; C = 0.72e-12
    alpha = 92.0; beta = 1.5; L=300e-9
    model = ModelTES.ShankRIT(alpha, beta, n, Tc, Tbath, k, R0, Rn);
    tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,model)
    bt = BiasedTES(tes_param, R0)
end

pixel = highEpixHighBeta()
params=pixel.p

iv_pt = ModelTES.iv_point(pixel.p, pixel.V);
IV_Vs = pixel.V*collect(0.05:0.05:5);
I,T,R,V = ModelTES.iv_curve(pixel.p, IV_Vs);

function iv_point_from_I(p::TESParams,I)
    k = p.k
    n = p.n
    Tbath = p.Tbath
    Rl = p.Rl
    function thermal(T)
        R = ModelTES.R(I,T,p)
        I^2*R-k*(T^n-Tbath^n)
    end
    T = fzero(thermal, p.Tc)
    R = ModelTES.R(I,T,p)
    V=I*(R+Rl)
    I,T,R,V
end
function iv_curve_from_Is(p::TESParams, Is)
    Is_out = Vector{Float64}(length(Is))
    Ts=Vector{Float64}(length(Is))
    Rs=Vector{Float64}(length(Is))
    Vs=Vector{Float64}(length(Is))
    for i in eachindex(Is)
            I,T,R,V = iv_point_from_T(p,Is[i])
            Is_out[i]=I
            Ts[i]=T
            Rs[i]=R
            Vs[i]=V
    end
    Is_out,Ts,Rs,Vs
end
function iv_curves_from_Is(p::TESParams, Tbaths, Is)
    params = deepcopy(p)
    Is = zeros(length(IV_Vs), length(Tbaths));
    Ts = zeros(length(IV_Vs), length(Tbaths));
    Vs = zeros(length(IV_Vs), length(Tbaths));
    Rs = zeros(length(IV_Vs), length(Tbaths));
    for (i,Tb) in enumerate(Tbaths)
        params.Tbath=Tb
        I,T,R,V = iv_curve_from_Ts(params, Is)
        Is[:,i]=I
        Ts[:,i]=T
        Rs[:,i]=R
        Vs[:,i]=V
    end
    return Is, Ts, Rs, Vs
end

function iv_point_from_T(p::TESParams,T, Iguess=1e-6)
    k = p.k
    n = p.n
    Tbath = p.Tbath
    Rl = p.Rl
    if Tbath!=T
        function thermal(I)
            R = ModelTES.R(I,T,p)
            I^2*R-k*(T^n-Tbath^n)
        end
        I = fzero(thermal, Iguess, [0, 1])
    else
        I=0.0
    end
    R = ModelTES.R(I,T,p)
    V=I*(R+Rl)
    I,T,R,V
end

function iv_curve_from_Ts(p::TESParams, Tmax, n)
    Ts = linspace(p.Tbath+0.005, Tmax,n)
    Is=Vector{Float64}(length(Ts))
    Ts_out=Vector{Float64}(length(Ts))
    Rs=Vector{Float64}(length(Ts))
    Vs=Vector{Float64}(length(Ts))
    for i in eachindex(Ts)
            I,T,R,V = iv_point_from_T(p,Ts[i])
            Is[i]=I
            Ts_out[i]=T
            Rs[i]=R
            Vs[i]=V
    end
    Is,Ts_out,Rs,Vs
end
function iv_curves_from_Ts(p::TESParams, Tbaths, Tmax, n)
    params = deepcopy(p)
    Is = zeros(n, length(Tbaths));
    Ts = zeros(n, length(Tbaths));
    Vs = zeros(n, length(Tbaths));
    Rs = zeros(n, length(Tbaths));
    for (i,Tb) in enumerate(Tbaths)
        params.Tbath=Tb
        I,T,R,V = iv_curve_from_Ts(params, Tmax,n)
        Is[:,i]=I
        Ts[:,i]=T
        Rs[:,i]=R
        Vs[:,i]=V
    end
    return Is, Ts, Rs, Vs
end

#IVs from Ts
Ts = 0.05:0.0005:0.15
Tbaths = 0.05:0.005:0.12
iv_point_from_T(pixel.p,0.106)
Is1, Ts1, Rs1, Vs1 = iv_curve_from_Ts(pixel.p, 0.15, 100)
Is, Ts, Rs, Vs = iv_curves_from_Ts(pixel.p, Tbaths, 0.15, 1000)
Vtess1 = Vs1-Is1*pixel.p.Rl;
Ps1 = Vtess1.*Is1;
Vtess = Vs-Is*pixel.p.Rl
Ps = Vtess.*Is

figure()
cmap = matplotlib[:cm][:get_cmap]("Reds")
xlim(0,0.0001)
ylim(0.095,0.11)
xlabel("current (A)")
ylabel("temperature (K)")
@pyimport scipy.interpolate as spinterpoplate
RIT  = spinterpoplate.interp2d(Is, Ts, Rs)
Inew = linspace(xlim()[1],xlim()[2],500)
Tnew = linspace(ylim()[1],ylim()[2],500)
IInew,TTnew=np.meshgrid(Inew,Tnew)
Rnew = RIT(Inew[:],Tnew[:])
cmap2 = matplotlib[:cm][:get_cmap]("Blues")
scatter(IInew[:],TTnew[:], color=cmap2(2*Rnew[:]/params.Rn))
scatter(Is[:],Ts[:],color=cmap(2*Rs[:]/params.Rn))

# imshow(cmap2(2*Rnew/params.Rn), extent = (xlim()[1], xlim()[2], ylim()[1], ylim()[2]))
#IVs from Is
Is = 0:1e-5:1e-3
iv_point_from_I(pixel.p, 1e-4)
iv_curve_from_Is(pixel.p, Is)

figure()
plot(Ps, Rs)

@pyimport scipy.interpolate as spinterpoplate
function iv_point_from_data_RIT(Rdata,Idata,Tdata, k, n, Tbath, Rl, T, Iguess=1e-6)
    RIT = spinterpoplate.interp2d(Idata, Tdata, Rdata)
    function thermal(I)
       @show R = RIT(I,T)[]
       @show R2 = ModelTES.R(I,T,params)
       @show I^2*R-k*(T^n-Tbath^n)
    end
    I = fzero(thermal, Iguess, [0, 1])
    R=RIT(I,T)[]
    V=I*(R+Rl)
    I,T,R,V
end

iv_point_from_data_RIT(Rs, Is, Ts, params.k, params.n, Tbaths[1], params.Rl, params.Tc )

# I: current in TES
# T: temperature of TES
# V: thevinen equivalen voltage
# Ibias =  V/Rsh
# Vtes = V-I*Rl
Vtes = V-I*pixel.p.Rl;
# plot(Vtes,I)
# xlabel!("V_{tes} (Volts)")
# ylabel!("Ites (Amps)")
# title!("Tb = $(pixel.p.Tbath) K")

# pulserecord = ModelTES.rk8(4000,1e-5,pixel,6000)
#
Tbaths = 0.05:0.005:0.11
params = deepcopy(pixel.p)
Is = zeros(length(IV_Vs), length(Tbaths));
Ts = zeros(length(IV_Vs), length(Tbaths));
Vs = zeros(length(IV_Vs), length(Tbaths));
Rs = zeros(length(IV_Vs), length(Tbaths));
for (i,Tb) in enumerate(Tbaths)
    params.Tbath=Tb
    I,T,R,V = iv_curve_from_Ts(params, 0.15, 100)
    Is[:,i]=I
    Ts[:,i]=T
    Rs[:,i]=R
    Vs[:,i]=V
end
Vtess = Vs-Is*params.Rl;
Ps = Vtess.*Is;
Rsrecalc = Vtess./Is;
Rpercs = Rs/params.Rn;

figure()
plot(Vtess, Is)
xlabel("V_{tes} (Volts)")
ylabel("I_{tes} (Amps)")
legend(loc="best", labels=map(repr, Tbaths))

figure()
plot(Rs,Ts)
xlabel("R_{tes} (Ohms)")
ylabel("T_{tes} (Kelvin)")
ylim(params.Tc-params.RIT.Tw*3, params.Tc+params.RIT.Tw*3)

figure()
plot(Ps,Rs)
xlabel("P_{tes} (Watts)")
ylabel("R_{tes} (Ohms)")
xlim(0,3e-12)
legend(loc="best", labels=map(repr, Tbaths))

# legend(loc="righ", labels=map(repr, Tbaths))
# p = plot(Vtess, Is, leg=false, lab=map(repr, Tbaths'))
# xlabel!("V_{tes} (Volts)")
# ylabel!("I_{tes} (Amps)")
#
# p3 = plot(Vs, Rpercs,  leg=true, lab=map(repr, Tbaths'))
# xlabel!("V_{thevinen} (Volts)")
# ylabel!("R (% Rn)")

# assume R is a function of T only, so R uniquely identifies temperature
# for a given value of R, say 0.4*Rn, we fit P as a function of Tb
"powercurves(Ts, Ps, Ts_chosen)
Ts a 2d array of temperatures from IV curves, where the second index for Tbath, first index is for Vbias
Ps a 2d array of dissipated powers from IV curves, same indicies
Ts_chosen a 1d Vector of temperatures chosen to interpolate power curves on
this may also be used when resistance instead of temperature, relying on the assumption that resistance is only a function of temperature
we know that isn't right so there will be errors assocaited with that"
function powercurves(Ts, Ps, Ts_chosen)
    @show size(Ts)
    @show size(Ps)
    P_interp = zeros(size(Ts,2), length(Ts_chosen))
    @show size(P_interp)
    for (i, T_chosen) in enumerate(Ts_chosen)
        @show i
        PvsTb = [np.interp(T_chosen, Ts[:,i], Ps[:,i]) for i in 1:size(Ps,2)]
        P_interp[:,i]=PvsTb
    end
    P_interp
end

Rpercs_chosen = 0.15:0.25:0.95;
P_interp = powercurves(Rpercs, Ps, Rpercs_chosen)

using LsqFit

Ptofit(Tbaths, T, k, n) =  k*(T^n-Tbaths.^n)
model(x, p) = Ptofit(x,p[1],p[2],p[3])

use_inds = Tbaths.<params.Tc-0.003
Rperc_ind=1
fit = curve_fit(model, Tbaths[use_inds], P_interp[use_inds,Rperc_ind]*1e12, [params.Tc, 3.5e3, 3.3])
Tfit, kfit, nfit = fit.param
Rperc_at_Tfit = Rpercs_chosen[Rperc_ind]

Ts_chosen = [params.Tc-params.RIT.Tw/2, params.Tc, params.Tc+params.RIT.Tw/2]
P_interpT = powercurves(Ts, Ps, Ts_chosen)
T_chosen_ind=1
fitT = curve_fit(model, Tbaths[use_inds], P_interpT[use_inds,T_chosen_ind]*1e12, [params.Tc, 1000, 3])
TfitT, kfitT, nfitT = fitT.param

# p5 = plot(Tbaths, P_interpT*1e12,lab=map(repr, Ts_chosen'),w=2)
# plot!(Tbaths[use_inds], model(Tbaths[use_inds], fit.param), lab="fit",style=:dash,w=2)
# xlabel!("T_{bath} (K)")
# ylabel!("Power required to to get to chosen %Rn (pW)")

j=7 # Tbaths index
TfromTb(Tbath, P, k, n) = (P/k+Tbath^n)^(1/n)
Ps+=1e-30
TrecalcT = [TfromTb(Tbaths[q], Ps[i,q]*1e12, kfitT, nfitT) for i in 1:size(Ps,1), q in 1:size(Ps,2)]
TrecalcR = [TfromTb(Tbaths[q], Ps[i,q]*1e12, kfit, nfit) for i in 1:size(Ps,1), q in 1:size(Ps,2)]

function iv_curve_residual(k,n,I,Vtes,Tbath)
    P = I*Vtes
    Tcalc = TfromTb(Tbath, P, k, n)
    Rcalc = ModelTES.RIT()

Ts_chosen_2 = [params.Tc-params.RIT.Tw, params.Tc, params.Tc+params.RIT.Tw]
P_interp_correct = powercurves(Ts, Ps, Ts_chosen_2)
figure()
plot(Tbaths, P_interp*1e12, lw=2)
plot(Tbaths[use_inds], model(Tbaths[use_inds], fit.param),"r--",lw=2)
plot(Tbaths, P_interp_correct*1e12,"-.",lw=3, label="correct curve")
xlabel("T_{bath} (Kelvin)")
ylabel("Power to reach resistance (pW)")
legend(loc="best", labels=push!(map(repr, Rpercs_chosen), "fit"))

figure()
plot(Tbaths, P_interpTrecalcR*1e12, label="P_interpTrecalcR")
plot(Tbaths, P_interp[:,Rperc_ind]*1e12, label="PinerpR")
legend(loc="best")
xlabel("Tbath")
ylabel("Power to reach %Rn or Temp (pW)")

P_interpTrecalcT = powercurves(TrecalcT, Ps, [TfitT])
figure()
plot(Tbaths, P_interpTrecalcT*1e12, label="P_interpTrecalcR")
plot(Tbaths, P_interpT[:,T_chosen_ind]*1e12, label="PinerpR")
legend(loc="best")
xlabel("Tbath")
ylabel("Power to reach %Rn or Temp (pW)")

# p6=plot(Ps[:,j]*1e12, Ts[:,j]*1e3,lab="Ts")
# plot!(Ps[:,j]*1e12, TrecalcT[:,j]*1e3, lab="TrecalcT")
# plot!(Ps[:,j]*1e12, TrecalcR[:,j]*1e3, lab="TrecalcR")
# xlabel!("power dissipated on TES (pW)")
# ylabel!("temperature of TES (mK)")
# title!("Tbath = $(Tbaths[j]*1e3) mK")
# print("lastline")
function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T})
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end
Is_for_mesh = linspace(0,maximum(Is),90);
Ts_for_mesh = linspace(minimum(Ts), maximum(Ts),100);
Imesh, Tmesh = meshgrid(Is_for_mesh, Ts_for_mesh);
Rdirect = [ModelTES.R(Is_for_mesh[i],Ts_for_mesh[j],params) for i=1:size(Is_for_mesh,1),j=1:size(Ts_for_mesh,1)]

# surface(Imesh,Tmesh,Rdirect ,xlabel="I",ylabel="T",zlabel="R",legend=false)
# surface!(Is, Trecalc, Rs, color=:red)
figure(figsize=(12,8))
surf(Is_for_mesh*1e6,Ts_for_mesh*1e3,1e3*Rdirect')
xlabel("I (uA)")
ylabel("T (mK)")
zlabel("R (mOhm)")
surf(Is*1e6, TrecalcT*1e3, Rs*1e3+0.001, color="cyan")
ax=gca()
ax[:plot](pulserecord.I*1e6, pulserecord.T*1e3, pulserecord.R*1e3, lw=4, color="red")
x1,y1,z1 = pulserecord.I[1]*1e6, pulserecord.T[1]*1e3, pulserecord.R[1]*1e3
x2,y2,z2 = pulserecord.I[end]*1e6, pulserecord.T[end]*1e3, pulserecord.R[end]*1e3
# ax[:quiver](x2,y2,z2,x2-x1,y2-y1,z2-z1)
ax[:plot]([x1,x2],[y1,y2],[z1,z2], color="orange",lw=4)
# ax[:plot]([x1],[y1],[z1], s=[10],color="black")

figure()
Ts_for_thisplot = Ts_chosen
Rs_for_thisplot = [ModelTES.R(Is_for_mesh[i], Ts_for_thisplot[j], params) for i in 1:size(Is_for_mesh,1), j in 1:size(Ts_for_thisplot,1)]
plot(Is_for_mesh, Rs_for_thisplot)
title("R vs T @ fixed I")
xlabel("current (A)")
ylabel("resistance (Ohms)")
legend(labels=map(repr,Ts_for_thisplot), loc="best")
println("lastline")

figure()
Ts_fit_calc = (Ps*1e12/kfit+repmat(Tbaths', size(Ts,1),1).^nfit).^(1/nfit);
surf(Is*1e6, Ts_fit_calc*1e3 ,Rs*1e3)
xlabel("I (uA)")
ylabel("T (mK)")
zlabel("R (mOhm)")

figure()
Ts_fit_calc = (Ps*1e12/kfit+repmat(Tbaths', size(Ts,1),1).^nfit).^(1/nfit);
surf(Is*1e6, Ts*1e3 ,Rs*1e3)
xlabel("I (uA)")
ylabel("T (mK)")
zlabel("R (mOhm)")

figure()
plot(Ts*1e3,  Ts_fit_calc*1e3)
plot()
xlabel("temperature true (mK)")
xlabel("temperature fit_calc (mK)")
legend(loc="best", labels=map(repr, Tbaths),title="Tbath")

figure()
plot(Ts*1e3,  (Ts_fit_calc-Ts)*1e3)
plot([params.Tc-params.RIT.Tw*3, params.Tc+params.RIT.Tw*3]*1e3, [1,1],"k", lw=3)
xlabel("temperature true (mK)")
ylabel("temperature fit_calc-true (mK)")
legend(loc="best", labels=push!(map(repr, Tbaths),"transition"),title="Tbath")


text:
__precompile__()

module AtlasRobot

using RigidBodyDynamics
using RigidBodyDynamics.Contact
using StaticArrays

packagepath() = joinpath(@__DIR__, "..", "deps")
urdfpath() = joinpath(packagepath(), "Atlas", "atlas.urdf")

function default_contact_model()
    SoftContactModel(hunt_crossley_hertz(k = 500e3), ViscoelasticCoulombModel(0.8, 20e3, 100.))
end

function flipsign_if_right(x::Number, side::Symbol)
    side == :left && return x
    side == :right && return -x
    error()
end

function mechanism(::Type{T} = Float64;
        floating = true,
        contactmodel = default_contact_model(),
        remove_fixed_tree_joints = true, add_flat_ground=false) where {T}
    mechanism = RigidBodyDynamics.parse_urdf(urdfpath(); scalar_type=T, floating=floating, remove_fixed_tree_joints=remove_fixed_tree_joints)

    if contactmodel != nothing
        for side in (:left, :right)
            foot = findbody(mechanism, "$(first(string(side)))_foot")
            frame = default_frame(foot)
            z = -0.07645

            # heel
            add_contact_point!(foot, ContactPoint(Point3D(frame, -0.0876, flipsign_if_right(0.066, side), z), contactmodel))
            add_contact_point!(foot, ContactPoint(Point3D(frame, -0.0876, flipsign_if_right(-0.0626, side), z), contactmodel))

            # toe:
            add_contact_point!(foot, ContactPoint(Point3D(frame, 0.1728, flipsign_if_right(0.066, side), z), contactmodel))
            add_contact_point!(foot, ContactPoint(Point3D(frame, 0.1728, flipsign_if_right(-0.0626, side), z), contactmodel))

            # midfoot:
            add_contact_point!(foot, ContactPoint(Point3D(frame, -0.0426, flipsign_if_right(0.066, side), z), contactmodel))
            add_contact_point!(foot, ContactPoint(Point3D(frame, -0.0426, flipsign_if_right(-0.0626, side), z), contactmodel))
        end
    end

    if add_flat_ground
        frame = root_frame(mechanism)
        ground = HalfSpace3D(Point3D(frame, 0., 0., 0.), FreeVector3D(frame, 0., 0., 1.))
        add_environment_primitive!(mechanism, ground)
    end

    remove_fixed_tree_joints && remove_fixed_tree_joints!(mechanism)

    return mechanism
end

function setnominal!(atlasstate::MechanismState)
    mechanism = atlasstate.mechanism
    zero!(atlasstate)
    kneebend = 1.1
    hipbendextra = 0.0
    left_shoulder_x_angle = -1.2
    elbow_y_angle = pi / 2
    left_elbow_x_angle = 0.2
    for sideprefix in ('l', 'r')
        knee = findjoint(mechanism, "$(sideprefix)_leg_kny")
        hippitch = findjoint(mechanism, "$(sideprefix)_leg_hpy")
        anklepitch = findjoint(mechanism, "$(sideprefix)_leg_aky")
        shoulder_x = findjoint(mechanism, "$(sideprefix)_arm_shx")
        elbow_y = findjoint(mechanism, "$(sideprefix)_arm_ely")
        elbow_x = findjoint(mechanism, "$(sideprefix)_arm_elx")
        set_configuration!(atlasstate, knee, kneebend)
        set_configuration!(atlasstate, hippitch, -kneebend / 2 + hipbendextra)
        set_configuration!(atlasstate, anklepitch, -kneebend / 2 - hipbendextra)
        set_configuration!(atlasstate, shoulder_x, sideprefix == 'l' ? left_shoulder_x_angle : -left_shoulder_x_angle)
        set_configuration!(atlasstate, elbow_y, elbow_y_angle)
        set_configuration!(atlasstate, elbow_x, sideprefix == 'l' ? left_elbow_x_angle : -left_elbow_x_angle)
    end
    floatingjoint = first(out_joints(root_body(mechanism), mechanism))
    set_configuration!(atlasstate, floatingjoint, [1; 0; 0; 0; 0; 0; 0.85])
    atlasstate
end

function __init__()
    if !isfile(urdfpath())
        error("Could not find $(urdfpath()). Please run `import Pkg; Pkg.build(\"AtlasRobot\")`.")
    end
end

end # module


text:
const gridfieldstosave = [:nx, :ny, :Lx, :Ly, :X, :Y]

"""
    Output(prob, filename, fieldtuples...)

Define output for the Problem prob with fields and functions that calculate
the output in (fldname, func) tuples.
"""
mutable struct Output
  prob::Problem
  filename::String
  fields::Dict{Symbol,Function}
  init::Bool
end

function Output(prob::Problem, filename::String, fields::Dict{Symbol,Function})
  if filename[end-4:end] == ".jld2" # Initialize output: remove trailing ".jld2" to form 'basefilename'
    basefilename = filename[1:end-5]
  else
    basefilename = filename
  end

  filename = basefilename*".jld2"
  n = 0
  while isfile(filename) # append numbers until unique name is found
    n += 1
    filename = basefilename * "_$n.jld2"
  end

  saveproblem(prob, filename)
  Output(prob, filename, fields, true)
end

function Output(prob::Problem, filename::String, fieldtuples...)
  Output(prob, filename, Dict{Symbol,Function}(
    [ (symfld[1], symfld[2]) for symfld in fieldtuples ]
  ))
end

getindex(out::Output, key) = out.fields[key](out.prob)

"""
    saveoutput(out)

Save current output fields for file in out.filename.
"""
function saveoutput(out::Output)
  groupname = "timeseries"
  jldopen(out.filename, "a+") do file
    file["$groupname/t/$(out.prob.step)"] = out.prob.t
    for fieldname in keys(out.fields)
      file["$groupname/$fieldname/$(out.prob.step)"] = out[fieldname]
    end
  end
  nothing
end

"""
    saveproblem(prob, filename)

Save certain aspects of a problem timestepper, grid, and params. Functions
that are fields in params are not saved.
"""
function saveproblem(prob::AbstractProblem, filename::String)
  jldopen(filename, "a+") do file
    file["timestepper/dt"] = prob.ts.dt   # Timestepper

    for field in gridfieldstosave         # Grid
      file["grid/$field"] = getfield(prob.grid, field)
    end

    for name in fieldnames(typeof(prob.params))   # Params
      field = getfield(prob.params, name)
      if !(typeof(field) <: Function)
        file["params/$name"] = field
      end
    end
  end
  nothing
end

saveproblem(out::Output) = saveproblem(out.prob, out.filename)

"""
    savediagnostic(diag, diagname)

Save diagnostics to file, labeled by the string diagname.
"""
function savediagnostic(diag::AbstractDiagnostic, diagname::String, filename::String)
  jldopen(filename, "a+") do file
    file["diags/$diagname/time"] = diag.time
    file["diags/$diagname/data"] = diag.data
  end
  nothing
end


text:
struct MPIError <: Exception
    code::Cint
end

macro mpichk(expr)
    @assert expr isa Expr && expr.head == :call && expr.args[1] == :ccall
    if Sys.iswindows()
        insert!(expr.args, 3, :stdcall)
    end
    :((errcode = $(esc(expr))) == MPI_SUCCESS || throw(MPIError(errcode)))
end

function error_string(err::MPIError)
    len_ref = Ref{Cint}()
    str_buf = Vector{UInt8}(undef, MPI_MAX_ERROR_STRING)
    # int MPI_Error_string(int errorcode, char *string, int *resultlen)
    @mpichk ccall((:MPI_Error_string, libmpi), Cint,
                  (Cint, Ptr{UInt8}, Ref{Cint}), err.code, str_buf, len_ref)
    resize!(str_buf, len_ref[])
    return String(str_buf)
end

function Base.show(io::IO, err::MPIError)
    print(io, "MPIError(", err.code, "): ", error_string(err))
end


text:
using NLsolve
using ForwardDiff
using LinearAlgebra
using Optim
using Printf
include("../src/algorithms/linesearches.jl")

# ELEMENTAL FUNCTIONS
function gfun!(F, x, u, p)
    # retrieve variables
    VM3 = x[1]
    VA3 = x[2]
    VA2 = x[3]

    VM1 = u[1]
    P2 = u[2]
    VM2 = u[3]

    VA1 = p[1]
    P3 = p[2]
    Q3 = p[3]

    # intermediate quantities
    VA23 = VA2 - VA3
    VA31 = VA3 - VA1
    VA32 = VA3 - VA2

    F[1] = 4.0*VM2*VM2 + VM2*VM3*(-4*cos(VA23) + 10*sin(VA23)) - P2
    F[2] = (8.0*VM3*VM3 + VM3*VM1*(-4*cos(VA31) + 5*sin(VA31))
            + VM3*VM2*(-4*cos(VA32) + 10*sin(VA32)) + P3)
    F[3] = (15.0*VM3*VM3 + VM3*VM1*(-4*sin(VA31) - 5*cos(VA31))
            + VM3*VM2*(-4*sin(VA32) - 10*cos(VA32)) + Q3)
end

function gfun(x, u, p, T=typeof(x))
    # Get Float64 of Vectors{Float64}, that is the first parameter
    F = zeros(T.parameters[1], 3)

    # retrieve variables
    VM3 = x[1]
    VA3 = x[2]
    VA2 = x[3]

    VM1 = u[1]
    P2 = u[2]
    VM2 = u[3]

    VA1 = p[1]
    P3 = p[2]
    Q3 = p[3]

    # intermediate quantities
    VA23 = VA2 - VA3
    VA31 = VA3 - VA1
    VA32 = VA3 - VA2

    F[1] = 4.0*VM2*VM2 + VM2*VM3*(-4*cos(VA23) + 10*sin(VA23)) - P2
    F[2] = (8.0*VM3*VM3 + VM3*VM1*(-4*cos(VA31) + 5*sin(VA31))
            + VM3*VM2*(-4*cos(VA32) + 10*sin(VA32)) + P3)
    F[3] = (15.0*VM3*VM3 + VM3*VM1*(-4*sin(VA31) - 5*cos(VA31))
            + VM3*VM2*(-4*sin(VA32) - 10*cos(VA32)) + Q3)
    return F
end

function cfun(x, u, p)
    VM3 = x[1]
    VA3 = x[2]

    VM1 = u[1]
    P2 = u[2]

    VA1 = p[1]
    VA13 = VA1 - VA3

    # we fix generation weights inside the
    # function to simplify the script and
    # follow the paper closely.
    w1 = 1.0
    w2 = 1.0

    cost = (w1*(4.0*VM1*VM1 + VM1*VM3*(-4*cos(VA13) + 5*sin(VA13))) +
            w2*P2)
    return cost
end

function pslack(x, u, p)
    VM3 = x[1]
    VA3 = x[2]

    VM1 = u[1]
    P2 = u[2]

    VA1 = p[1]
    VA13 = VA1 - VA3

    return (4.0*VM1*VM1 + VM1*VM3*(-4*cos(VA13) + 5*sin(VA13)))
end

# OPF COMPUTATION
function solve_pf(x, u, p, verbose=true)
    fun_pf!(F, x) = gfun!(F, x, u, p)
    x0 = copy(x)
    res = nlsolve(fun_pf!, x0)
    if verbose
        show(res)
        println("")
    end
    return res.zero
end

# initial parameters
x = zeros(3)
u = zeros(3)
p = zeros(3)

# this is an initial guess
x[1] = 1.0 #VM3
x[2] = 0.0 #VA3
x[3] = 0.0 #VA2

# this is given by the problem data, but might be "controlled" via OPF
u[1] = 1.0 #VM1
u[2] = 1.7 #P2
u[3] = 1.0 #VM2

# these parameters are fixed through the computation
p[1] = 0.0 #VA1, slack angle
p[2] = 2.0 #P3
p[3] = 1.0 #Q3

# print initial guesses
println(x)
println(u)

# copy to iteration vectors
xk = copy(x)
uk = copy(u)

function build_callback(x, p)
    xk = copy(x)
    function eval_fg!(f, g, uk)
        ## Inner loop
        xk = solve_pf(xk, uk, p, false)
        gx_x(x) = gfun(x, uk, p, typeof(x))
        gx = x -> ForwardDiff.jacobian(gx_x, x)
        # gradient
        cfun_x(x) = cfun(x, uk, p)
        fx = x ->ForwardDiff.gradient(cfun_x, x)
        # multipliers calculation
        lambda = -inv(gx(xk)')*fx(xk)

        # compute g_u, g_u
        cfun_u(u) = cfun(xk, u, p)
        fu = u ->ForwardDiff.gradient(cfun_u, u)
        gx_u(u) = gfun(xk, u, p, typeof(u))
        gu = u -> ForwardDiff.jacobian(gx_u, u)
        # compute gradient of cost function
        Lu = u -> cfun_u(u) + gx_u(u)'*lambda

        if g != nothing
            grad_Lu = u -> (fu(u) + gu(u)'*lambda)
            g .= grad_Lu(uk)
        end
        if f != nothing
            return Lu(uk)
        end
    end

    return eval_fg!
end

u0 = copy(u)
eval_fg! = build_callback(x, p)
options = Optim.Options(iterations=250, g_tol=1e-5, show_trace=true,
                        allow_f_increases=false)
res = Optim.optimize(Optim.only_fg!(eval_fg!),
                     u0, BFGS(), options)