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

mutable struct VirtualMachineProperties <: SwaggerModel
    hardwareProfile::Any # spec type: Union{ Nothing, HardwareProfile } # spec name: hardwareProfile
    storageProfile::Any # spec type: Union{ Nothing, StorageProfile } # spec name: storageProfile
    additionalCapabilities::Any # spec type: Union{ Nothing, AdditionalCapabilities } # spec name: additionalCapabilities
    osProfile::Any # spec type: Union{ Nothing, OSProfile } # spec name: osProfile
    networkProfile::Any # spec type: Union{ Nothing, NetworkProfile } # spec name: networkProfile
    securityProfile::Any # spec type: Union{ Nothing, SecurityProfile } # spec name: securityProfile
    diagnosticsProfile::Any # spec type: Union{ Nothing, DiagnosticsProfile } # spec name: diagnosticsProfile
    availabilitySet::Any # spec type: Union{ Nothing, SubResource } # spec name: availabilitySet
    virtualMachineScaleSet::Any # spec type: Union{ Nothing, SubResource } # spec name: virtualMachineScaleSet
    proximityPlacementGroup::Any # spec type: Union{ Nothing, SubResource } # spec name: proximityPlacementGroup
    priority::Any # spec type: Union{ Nothing, Priority } # spec name: priority
    evictionPolicy::Any # spec type: Union{ Nothing, EvictionPolicy } # spec name: evictionPolicy
    billingProfile::Any # spec type: Union{ Nothing, BillingProfile } # spec name: billingProfile
    host::Any # spec type: Union{ Nothing, SubResource } # spec name: host
    hostGroup::Any # spec type: Union{ Nothing, SubResource } # spec name: hostGroup
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    instanceView::Any # spec type: Union{ Nothing, VirtualMachineInstanceView } # spec name: instanceView
    licenseType::Any # spec type: Union{ Nothing, String } # spec name: licenseType
    vmId::Any # spec type: Union{ Nothing, String } # spec name: vmId
    extensionsTimeBudget::Any # spec type: Union{ Nothing, String } # spec name: extensionsTimeBudget

    function VirtualMachineProperties(;hardwareProfile=nothing, storageProfile=nothing, additionalCapabilities=nothing, osProfile=nothing, networkProfile=nothing, securityProfile=nothing, diagnosticsProfile=nothing, availabilitySet=nothing, virtualMachineScaleSet=nothing, proximityPlacementGroup=nothing, priority=nothing, evictionPolicy=nothing, billingProfile=nothing, host=nothing, hostGroup=nothing, provisioningState=nothing, instanceView=nothing, licenseType=nothing, vmId=nothing, extensionsTimeBudget=nothing)
        o = new()
        validate_property(VirtualMachineProperties, Symbol("hardwareProfile"), hardwareProfile)
        setfield!(o, Symbol("hardwareProfile"), hardwareProfile)
        validate_property(VirtualMachineProperties, Symbol("storageProfile"), storageProfile)
        setfield!(o, Symbol("storageProfile"), storageProfile)
        validate_property(VirtualMachineProperties, Symbol("additionalCapabilities"), additionalCapabilities)
        setfield!(o, Symbol("additionalCapabilities"), additionalCapabilities)
        validate_property(VirtualMachineProperties, Symbol("osProfile"), osProfile)
        setfield!(o, Symbol("osProfile"), osProfile)
        validate_property(VirtualMachineProperties, Symbol("networkProfile"), networkProfile)
        setfield!(o, Symbol("networkProfile"), networkProfile)
        validate_property(VirtualMachineProperties, Symbol("securityProfile"), securityProfile)
        setfield!(o, Symbol("securityProfile"), securityProfile)
        validate_property(VirtualMachineProperties, Symbol("diagnosticsProfile"), diagnosticsProfile)
        setfield!(o, Symbol("diagnosticsProfile"), diagnosticsProfile)
        validate_property(VirtualMachineProperties, Symbol("availabilitySet"), availabilitySet)
        setfield!(o, Symbol("availabilitySet"), availabilitySet)
        validate_property(VirtualMachineProperties, Symbol("virtualMachineScaleSet"), virtualMachineScaleSet)
        setfield!(o, Symbol("virtualMachineScaleSet"), virtualMachineScaleSet)
        validate_property(VirtualMachineProperties, Symbol("proximityPlacementGroup"), proximityPlacementGroup)
        setfield!(o, Symbol("proximityPlacementGroup"), proximityPlacementGroup)
        validate_property(VirtualMachineProperties, Symbol("priority"), priority)
        setfield!(o, Symbol("priority"), priority)
        validate_property(VirtualMachineProperties, Symbol("evictionPolicy"), evictionPolicy)
        setfield!(o, Symbol("evictionPolicy"), evictionPolicy)
        validate_property(VirtualMachineProperties, Symbol("billingProfile"), billingProfile)
        setfield!(o, Symbol("billingProfile"), billingProfile)
        validate_property(VirtualMachineProperties, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(VirtualMachineProperties, Symbol("hostGroup"), hostGroup)
        setfield!(o, Symbol("hostGroup"), hostGroup)
        validate_property(VirtualMachineProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(VirtualMachineProperties, Symbol("instanceView"), instanceView)
        setfield!(o, Symbol("instanceView"), instanceView)
        validate_property(VirtualMachineProperties, Symbol("licenseType"), licenseType)
        setfield!(o, Symbol("licenseType"), licenseType)
        validate_property(VirtualMachineProperties, Symbol("vmId"), vmId)
        setfield!(o, Symbol("vmId"), vmId)
        validate_property(VirtualMachineProperties, Symbol("extensionsTimeBudget"), extensionsTimeBudget)
        setfield!(o, Symbol("extensionsTimeBudget"), extensionsTimeBudget)
        o
    end
end # type VirtualMachineProperties

const _property_map_VirtualMachineProperties = Dict{Symbol,Symbol}(Symbol("hardwareProfile")=>Symbol("hardwareProfile"), Symbol("storageProfile")=>Symbol("storageProfile"), Symbol("additionalCapabilities")=>Symbol("additionalCapabilities"), Symbol("osProfile")=>Symbol("osProfile"), Symbol("networkProfile")=>Symbol("networkProfile"), Symbol("securityProfile")=>Symbol("securityProfile"), Symbol("diagnosticsProfile")=>Symbol("diagnosticsProfile"), Symbol("availabilitySet")=>Symbol("availabilitySet"), Symbol("virtualMachineScaleSet")=>Symbol("virtualMachineScaleSet"), Symbol("proximityPlacementGroup")=>Symbol("proximityPlacementGroup"), Symbol("priority")=>Symbol("priority"), Symbol("evictionPolicy")=>Symbol("evictionPolicy"), Symbol("billingProfile")=>Symbol("billingProfile"), Symbol("host")=>Symbol("host"), Symbol("hostGroup")=>Symbol("hostGroup"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("instanceView")=>Symbol("instanceView"), Symbol("licenseType")=>Symbol("licenseType"), Symbol("vmId")=>Symbol("vmId"), Symbol("extensionsTimeBudget")=>Symbol("extensionsTimeBudget"))
const _property_types_VirtualMachineProperties = Dict{Symbol,String}(Symbol("hardwareProfile")=>"HardwareProfile", Symbol("storageProfile")=>"StorageProfile", Symbol("additionalCapabilities")=>"AdditionalCapabilities", Symbol("osProfile")=>"OSProfile", Symbol("networkProfile")=>"NetworkProfile", Symbol("securityProfile")=>"SecurityProfile", Symbol("diagnosticsProfile")=>"DiagnosticsProfile", Symbol("availabilitySet")=>"SubResource", Symbol("virtualMachineScaleSet")=>"SubResource", Symbol("proximityPlacementGroup")=>"SubResource", Symbol("priority")=>"Priority", Symbol("evictionPolicy")=>"EvictionPolicy", Symbol("billingProfile")=>"BillingProfile", Symbol("host")=>"SubResource", Symbol("hostGroup")=>"SubResource", Symbol("provisioningState")=>"String", Symbol("instanceView")=>"VirtualMachineInstanceView", Symbol("licenseType")=>"String", Symbol("vmId")=>"String", Symbol("extensionsTimeBudget")=>"String")
Base.propertynames(::Type{ VirtualMachineProperties }) = collect(keys(_property_map_VirtualMachineProperties))
Swagger.property_type(::Type{ VirtualMachineProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineProperties }, property_name::Symbol) =  _property_map_VirtualMachineProperties[property_name]

function check_required(o::VirtualMachineProperties)
    true
end

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


text:
is_empty(bitset::RSparseBitSet) = bitset.last_ptr == 0

function clear_mask(bitset::RSparseBitSet)
    indices = bitset.indices
    mask = bitset.mask
    @inbounds for i in 1:(bitset.last_ptr)
        idx = indices[i]
        mask[idx] = UInt64(0)
    end
end

function full_mask(bitset::RSparseBitSet)
    indices = bitset.indices
    mask = bitset.mask
    @inbounds for i in 1:length(indices)
        idx = indices[i]
        mask[idx] = typemax(UInt64)
    end
end

function invert_mask(bitset::RSparseBitSet)
    indices = bitset.indices
    mask = bitset.mask
    @inbounds for i in 1:(bitset.last_ptr)
        idx = indices[i]
        mask[idx] = ~mask[idx]
    end
end

function add_to_mask(bitset::RSparseBitSet, add)
    mask = bitset.mask
    indices = bitset.indices
    @inbounds for i in 1:(bitset.last_ptr)
        idx = indices[i]
        mask[idx] |= add[idx]
    end
end

function intersect_mask_with_mask(bitset::RSparseBitSet, intersect_mask)
    indices = bitset.indices
    mask = bitset.mask
    @inbounds for i in 1:(bitset.last_ptr)
        idx = indices[i]
        mask[idx] &= intersect_mask[idx]
    end
end

function intersect_with_mask_feasible(bitset::RSparseBitSet)
    words = bitset.words
    mask = bitset.mask
    indices = bitset.indices
    @inbounds for i in (bitset.last_ptr):-1:1
        idx = indices[i]
        w = words[idx] & mask[idx]
        if w != UInt64(0)
            return true # is feasible
        end
    end
    return false
end

function intersect_with_mask(bitset::RSparseBitSet)
    words = bitset.words
    mask = bitset.mask
    indices = bitset.indices
    @inbounds for i in (bitset.last_ptr):-1:1
        idx = indices[i]
        w = words[idx] & mask[idx]
        if w != words[idx]
            words[idx] = w
            if w == zero(UInt64)
                indices[i] = indices[bitset.last_ptr]
                indices[bitset.last_ptr] = idx
                bitset.last_ptr -= 1
            end
        end
    end
end

function intersect_index(bitset::RSparseBitSet, mask)
    words = bitset.words
    indices = bitset.indices
    @inbounds for i in 1:(bitset.last_ptr)
        idx = indices[i]
        w = words[idx] & mask[idx]
        if w != zero(UInt64)
            return idx
        end
    end
    return 0
end

ext: jl


text:
a=[.4 .6; .3 .7]
b=[.3 .1 .6; .5 .2 .3]
data=[1,2,3]
model=hmm((a,b))

function naiveprob(data)
	filter=transpose(model.mu[:,data[1]])
	for t=2:length(data)
		filter=sum(filter*slice(model.m,:,data[t-1],:,data[t]),1)
	end
	sum(filter)
end

@test vec(sum(model.m,(3,4)))≈fill(1,6)
@test reshape(sum(model.m[:,1,:,:],4),2,2)≈a
@test loglikelihood(model,data)≈log(naiveprob(data))/length(data)
data=Array(Array,2)
data[1]=[1,2,3]
data[2]=[3,2,1]
@test loglikelihood(model,data)≈(log(naiveprob(data[1]))/length(data[1])+log(naiveprob(data[2]))/length(data[2]))/2

viterbi(model,data)

data=rand(model,10000)
abhat=em(model,data)

@test vecnorm(abhat[1]-a)<5e-2
@test vecnorm(abhat[2]-b)<5e-2


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct TensorProto <: ProtoType
    dtype::Int32
    tensor_shape::TensorShapeProto
    version_number::Int32
    tensor_content::Array{UInt8,1}
    half_val::Base.Vector{Int32}
    float_val::Base.Vector{Float32}
    double_val::Base.Vector{Float64}
    int_val::Base.Vector{Int32}
    string_val::Base.Vector{Array{UInt8,1}}
    scomplex_val::Base.Vector{Float32}
    int64_val::Base.Vector{Int64}
    bool_val::Base.Vector{Bool}
    dcomplex_val::Base.Vector{Float64}
    resource_handle_val::Base.Vector{ResourceHandleProto}
    variant_val::Base.Any
    uint32_val::Base.Vector{UInt32}
    uint64_val::Base.Vector{UInt64}
    TensorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorProto (has cyclic type dependency)
const __fnum_TensorProto = Int[1,2,3,4,13,5,6,7,8,9,10,11,12,14,15,16,17]
const __pack_TensorProto = Symbol[:half_val,:float_val,:double_val,:int_val,:scomplex_val,:int64_val,:bool_val,:dcomplex_val,:uint32_val,:uint64_val]
const __ftype_TensorProto = Dict(:variant_val => "Base.Vector{VariantTensorDataProto}")
meta(t::Type{TensorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_TensorProto, ProtoBuf.DEF_VAL, true, __pack_TensorProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, __ftype_TensorProto)

mutable struct VariantTensorDataProto <: ProtoType
    type_name::AbstractString
    metadata::Array{UInt8,1}
    tensors::Base.Vector{TensorProto}
    VariantTensorDataProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct VariantTensorDataProto (has cyclic type dependency)

export TensorProto, VariantTensorDataProto, TensorProto, VariantTensorDataProto


text:
module VarintTranslator # LoRDeckCodes

using ..LoRDeckCodes: ArgumentException

const AllButMSB = 0x7f
const JustMSB   = 0x80

function pop_varint(stream::IOBuffer)::Int
    result::UInt = 0
    currentShift = 0
    bytesPopped = 0
    while !eof(stream)
        (byte,) = read(stream, 1)
        current = byte & AllButMSB
        result |= current << currentShift
        (byte & JustMSB) != JustMSB && return Int(result)
        currentShift += 7
    end
    throw(ArgumentException("Byte array did not contain valid varints."))
end

function get_varint(value::Int)::Vector{UInt8}
    buff = zeros(UInt8, 10)
    currentIndex = 1
    iszero(value) && return UInt8[0]
    while !iszero(value)
        byteVal = value & AllButMSB
        value >>= 7
        if !iszero(value)
            byteVal |= JustMSB
        end
        buff[currentIndex] = byteVal
        currentIndex += 1
    end
    buff[1:currentIndex-1]
end

end # module LoRDeckCodes.VarintTranslator


text:
## ExoplanetsSysSim.jl
## (c) 2015 Eric B. Ford

module SimulationParameters

import Compat: @compat, readstring
import Pkg, LibGit2
using ExoplanetsSysSim

export SimParam, add_param_fixed, add_param_active, update_param, set_active, set_inactive, is_active
export get_any, get_real, get_int, get_bool, get_function

export make_vector_of_active_param_keys, make_vector_of_sim_param, get_range_for_sim_param, update_sim_param_from_vector!
export setup_sim_param_demo, test_sim_param_constructors
#export preallocate_memory!

const global version_id_str = "0.3.1"
const global version_id_pair = ("version",version_id_str)
const global julia_version_pair = ("version_julia",string(VERSION))

mutable struct SimParam
  param::Dict{String,Any}
  active::Dict{String,Bool}
end
copy(p::SimParam) = SimParam(copy(p.param),copy(p.active))

"""
    SimParam(p::Dict{String,Any})
Creates a SimParam object from the dictionary p, with all parameter defaulting to inactive.
"""
function SimParam(p::Dict{String,Any})   # By default all parameters are set as inactive (i.e., not allowed to be optimized)
  a = Dict{String,Bool}()
  for k in keys(p)
    a[k] = false
  end
  return SimParam(p,a)
end

"""
    SimParam()
Creates a nearly empty SimParam object, with just the version id and potentially other information about the code, system, runtime, etc.
"""
function SimParam()
  d = Dict{String,Any}([version_id_pair,julia_version_pair,("hostname",gethostname()), ("time",time())])
  try
    d["ExoplanetsSysSim version"] = Pkg.installed()["ExoplanetsSysSim"]
    #d["ExoplanetsSysSim directory"] = dirname(pathof(ExoplanetsSysSim))
    #d["ExoplanetsSysSim branch"] = LibGit2.headname(LibGit2.GitRepo(dirname(pathof(ExoplanetsSysSim))))
    #d["ExoplanetsSysSim head_oid"] = LibGit2.head_oid(LibGit2.GitRepo(dirname(pathof(ExoplanetsSysSim))))
  catch
    warn("# Couldn't add full information about version of SysSim to SimParam.")
  end
  SimParam(d)
end

"""
    add_param_fixed(sim::SimParam, key::String,val::Any)
Adds (or overwrites) key with value val to the SimParam object, sim, and sets the parameter set to inactive.
"""
function add_param_fixed(sim::SimParam, key::String,val::Any)
  sim.param[key] = val
  sim.active[key] = false
end

"""
### add_param_active(sim::SimParam, key::String,val::Any)
Adds (or overwrites) key with value val to the SimParam object, sim, and sets the parameter set to active.
"""
function add_param_active(sim::SimParam, key::String,val::Any)
  sim.param[key] = val
  sim.active[key] = true
end

"""
### update_param(sim::SimParam, key::String,val::Any)
Overwrites key with value val to the SimParam object, sim
"""
function update_param(sim::SimParam, key::String,val::Any)
  @assert haskey(sim.param,key)
  sim.param[key] = val
end

"""
### set_active(sim::SimParam, key::String)
Sets the key parameter to be active in sim.
"""
function set_active(sim::SimParam,key::String)
  @assert haskey(sim.param,key)
  sim.active[key] = true
end

"""
### set_active(sim::SimParam, keys::Vector{String})
Sets each of the key parameters to be active in sim.
"""
function set_active(sim::SimParam,keys::Vector{String})
  for k in keys
    set_active(sim,k)
  end
end

"""
### set_inactive(sim::SimParam, key::String)
Sets the key parameter to be inactive in sim.
"""
function set_inactive(sim::SimParam,key::String)
  @assert haskey(sim.param,key)
  sim.active[key] = false
end

"""
### set_inactive(sim::SimParam, keys::Vector{String})
Sets each of the key parameters to be inactive in sim.
"""
function set_inactive(sim::SimParam,keys::Vector{String})
  for k in keys
    set_inactive(sim,k)
  end
end

function is_active(sim::SimParam,key::String)
  @assert haskey(sim.active,key)
  sim.active[key]
end

import Base.get
function get(sim::SimParam, key::String, default_val::T) where T
  val::T = get(sim.param,key,default_val)::T
  return val
end

function get_any(sim::SimParam, key::String, default_val::Any)
  val = get(sim.param,key,default_val)
  return val
end

function get_real(sim::SimParam, key::String)
  val::Float64 = get(sim.param,key,convert(Float64,NaN) )::Float64
  @assert(val!=convert(Float64,NaN))
  return val
end

function get_int(sim::SimParam, key::String)
  val::Int64 = get(sim.param,key,zero(Int64))
  #@assert(val!=nan(zero(Int64)))
  #@assert(val!=oftype(x,NaN))
  return val
end

function get_bool(sim::SimParam, key::String)
  val::bool = get(sim.param,key,false)
  return val
end

function noop()
end

function get_function(sim::SimParam, key::String)
  val::Function = Base.get(sim.param,key,noop)::Function
  #val = Base.get(sim.param,key,null)
  @assert((val!=nothing) && (val!=noop))
  return val::Function
end

import Base.haskey
haskey(sim::SimParam, key::String) = return haskey(sim.param,key)

function make_vector_of_active_param_keys(sim::SimParam)
  sortedkeys = sort(collect(keys(sim.param)))
  sortedkeys_active = sortedkeys[map(k->get(sim.active,k,false),sortedkeys)]
  return sortedkeys_active
end

function make_vector_of_sim_param(sim::SimParam)
  param_vector = Float64[]      # QUERY: Currently we make separate vectors of Ints and Floats.  Does this make sense?
  sss = sort(collect(keys(sim.param)))
  for k in 1:length(sss)
    if(sim.active[sss[k]]==false)
      continue
    end
    if(length(sim.param[sss[k]])==1)
      if isa( sim.param[sss[k]], Real )
        push!(param_vector,sim.param[sss[k]])
      elseif(eltype( sim.param[sss[k]]) <: Real)
        append!(param_vector,vec(sim.param[sss[k]]))
      end
    elseif(eltype( sim.param[sss[k]]) <: Real)
      append!(param_vector, vec(sim.param[sss[k]]))
    else
      if eltype( sss[k]) <: Real
        append!(param_vector,sim.param[sss[k]])
      end
    end
  end
  return param_vector
end

function get_range_for_sim_param(key::String, sim::SimParam)
  sorted_keys = sort(collect(keys(sim.param)))
  i = 1
  for k in 1:length(sorted_keys)
    if(sim.active[sorted_keys[k]]==false)
      continue
    end
    param_len = length(sim.param[sorted_keys[k]])
    if sorted_keys[k]==key
       return i:(i+param_len)
    else
       #println("Didn't match >",sorted_keys[k],"< and >",key,"<.")
       i += param_len
    end
  end
  println("# ERROR: Never found range for param: ",key)
  return 0:0
end

function update_sim_param_from_vector!(param::Vector{Float64}, sim::SimParam)
  #println("# Input vector: ",param)
  sorted_keys = sort(collect(keys(sim.param)))
  i = 1
  for k in 1:length(sorted_keys)
    if(sim.active[sorted_keys[k]]==false)
      continue
    end
    param_len = length(sim.param[sorted_keys[k]])
    if param_len==1
      if isa( sim.param[sorted_keys[k]], Real )
        # println("# Replacing >",sorted_keys[k],"< with >",param[i],"<")
        sim.param[sorted_keys[k]] = param[i]
        i = i+1
      elseif eltype( sim.param[sorted_keys[k]]) <: Real
        # println("# Replacing >",sim.param[sorted_keys[k]],"< with >",reshape(param[i:i+param_len-1], size(sim.param[sorted_keys[k]]),"<")
        sim.param[sorted_keys[k]] = reshape(param[i:i+param_len-1], size(sim.param[sorted_keys[k]]))
        i = i+1
      end
    elseif param_len>1
      if eltype( sim.param[sorted_keys[k]]) <: Real
        # println("# Replacing >",sim.param[sorted_keys[k]],"< with >",reshape(param[i:i+param_len-1], size(sim.param[sorted_keys[k]])),"<")
        sim.param[sorted_keys[k]] = reshape(param[i:i+param_len-1], size(sim.param[sorted_keys[k]]))
        i = i+param_len
      end
    else
      println("# Don't know what to do with empty simulation parameter: ",sorted_keys[k])
    end
  end
  return sim
end

function preallocate_memory!(sim_param::SimParam)
  num_kepler_targets = get_int(sim_param,"num_kepler_targets")
  add_param_fixed(sim_param,"mem_kepler_target_obs", Array{KeplerTargetObs}(num_kepler_targets) )
end

function setup_sim_param_demo(args::Vector{String} = Array{String}(undef,0) )   # allow this to take a list of parameter (e.g., from command line)
  sim_param = SimParam()
  add_param_fixed(sim_param,"max_tranets_in_sys",7)
  add_param_fixed(sim_param,"num_targets_sim_pass_one",190000)                      # Note this is used for the number of stars in the simulations, not necessarily related to number of Kepler targets
  add_param_fixed(sim_param,"num_kepler_targets",190000)                           # Note this is used for the number of Kepler targets for the observational catalog
  add_param_fixed(sim_param,"generate_star",ExoplanetsSysSim.generate_star_dumb)
  #add_param_fixed(sim_param,"generate_planetary_system", ExoplanetsSysSim.generate_planetary_system_simple)
  add_param_fixed(sim_param,"generate_planetary_system", ExoplanetsSysSim.generate_planetary_system_uncorrelated_incl)

  # add_param_fixed(sim_param,"generate_kepler_target",ExoplanetsSysSim.generate_kepler_target_simple)
  add_param_fixed(sim_param,"generate_kepler_target",ExoplanetsSysSim.generate_kepler_target_from_table)
  add_param_fixed(sim_param,"star_table_setup",StellarTable.setup_star_table)
  add_param_fixed(sim_param,"stellar_catalog","q1q17_dr25_gaia_fgk.jld2")
  add_param_fixed(sim_param,"generate_num_planets",ExoplanetsSysSim.generate_num_planets_poisson)
  add_param_active(sim_param,"log_eta_pl",log(2.0))
  add_param_fixed(sim_param,"generate_planet_mass_from_radius",ExoplanetsSysSim.generate_planet_mass_from_radius_powerlaw)
  add_param_fixed(sim_param,"mr_power_index",2.0)
  add_param_fixed(sim_param,"mr_const",1.0)
  #add_param_fixed(sim_param,"generate_period_and_sizes",ExoplanetsSysSim.generate_period_and_sizes_log_normal)
  #add_param_active(sim_param,"mean_log_planet_radius",log(2.0*earth_radius))
  #add_param_active(sim_param,"sigma_log_planet_radius",log(2.0))
  #add_param_active(sim_param,"mean_log_planet_period",log(5.0))
  #add_param_active(sim_param,"sigma_log_planet_period",log(2.0))
  add_param_fixed(sim_param,"generate_period_and_sizes", ExoplanetsSysSim.generate_period_and_sizes_power_law)
  add_param_active(sim_param,"power_law_P",0.3)
  add_param_active(sim_param,"power_law_r",-2.44)
  add_param_fixed(sim_param,"min_period",1.0)
  add_param_fixed(sim_param,"max_period",100.0)
  add_param_fixed(sim_param,"min_radius",0.5*ExoplanetsSysSim.earth_radius)
  add_param_fixed(sim_param,"max_radius",10.0*ExoplanetsSysSim.earth_radius)
  add_param_fixed(sim_param,"generate_e_omega",ExoplanetsSysSim.generate_e_omega_rayleigh)
  add_param_fixed(sim_param,"sigma_hk",0.03)
  add_param_fixed(sim_param,"sigma_incl",2.0)   # degrees
  add_param_fixed(sim_param,"calc_target_obs_sky_ave",ExoplanetsSysSim.calc_target_obs_sky_ave)
  add_param_fixed(sim_param,"calc_target_obs_single_obs",ExoplanetsSysSim.calc_target_obs_single_obs)
  add_param_fixed(sim_param,"read_target_obs",ExoplanetsSysSim.simulated_read_kepler_observations)
  add_param_fixed(sim_param,"transit_noise_model",ExoplanetsSysSim.transit_noise_model_fixed_noise)
  # add_param_fixed(sim_param,"transit_noise_model",transit_noise_model_diagonal)
  # add_param_fixed(sim_param,"rng_seed",1234)   # If you want to be able to reproduce simulations

  # Do other initialization tasks belong here or elsewhere?
  # TODO OPT:  Try to preallocate memory for each target to see if this makes a performance difference
  # preallocate_memory!(sim_param)

  return sim_param
end

function test_sim_param_constructors()
  oldval = log(2.0)
  sim_param = SimParam( Dict([ ("version",version_id_str), ("num_kepler_targets",190000), ("log_eta_pl",oldval), ("max_tranets_in_sys",7)] ) )
  get(sim_param,"version","")
  set_active(sim_param,"log_eta_pl")
  sp_vec = make_vector_of_sim_param(sim_param)
  sp_vec .+= 0.1
  update_sim_param_from_vector!(sp_vec,sim_param)
  newval = get_real(sim_param,"log_eta_pl")
  isapprox(oldval+0.1,newval,atol=0.001)
end

end

#test_sim_param_constructors()


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

name = "Gnome_themes_extra"
version = v"3.28"

# Collection of sources required to build gnome-themes-extra
sources = [
    ArchiveSource("http://ftp.gnome.org/pub/gnome/sources/gnome-themes-extra/$(version.major).$(version.minor)/gnome-themes-extra-$(version.major).$(version.minor).tar.xz",
                  "7c4ba0bff001f06d8983cfc105adaac42df1d1267a2591798a780bac557a5819"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gnome-themes-extra-*/
apk add intltool

# Clear out `.la` files since they're often wrong and screw us up
rm -f ${prefix}/lib/*.la

# Remove a library from the host filesystem that is accidentally picked up by
# the build when doing a native compilation
rm /usr/lib/libexpat.so.1

./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-gtk2-engine
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
#platforms = supported_platforms()

# Limit to the same platforms as Gtk for now
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("share/themes/Adwaita/index.theme", :adwaita_index),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("GTK3_jll"),
    Dependency("Librsvg_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# this file will have the "master" function that's actually exposed to the user
include("pmg_graphs.jl")
include("smiles_rdkit_features.jl")

function featurize_inputs(input_data_folder, SMILES_atom_features, SMILES_bond_features, other_atom_features, graphs_folder, features_folder)
    # determine type of data in input folder (must be all SMILES or all CIFs)
    # check whether user has asked for features that make sense (may require some ancillary fcns)
    # loop over input files and farm out actual parsing to specialized functions for CIF vs. SMILES
    # eventually, save out adjacency matrices to graphs_folder and feature vectors/metadata to features_folder
end

# downstream methods for concatenating atom features generated by using SMILES
concat_smiles_atom_feat(atom_feat::OrderedDict) = vcat(values(atom_feat)...)
concat_smiles_bond_feat(bond_feat::OrderedDict) = vcat(values(bond_feat)...)
## Weave Featurization

# data structure for storing features returned from weave_featurizer
mutable struct WeaveMol
    atom_features
    pair_features
    num_atoms
    atom_feat_num
end

"""
weave featurizer

input:
    inputs(Array{String}): List of SIMLES strings
    atom_feature_list(Array{String}): list of atom features of interest
    bond_feature_list(Array{String}): list of bond features of interest
"""
function weave_featurize(inputs; atom_feature_list::Array{String}=atom_feat_list, bond_feature_list::Array{String}=bond_feat_list)
    outputs = Array{WeaveMol}(undef,length(inputs))
    for (i,smiles) in enumerate(inputs)
        mol = chem.MolFromSmiles(smiles)
        new_order = chem.rdmolfiles.CanonicalRankAtoms(mol)  # reorder the atoms so that they're always in the same canonical order.
        mol = chem.rdmolops.RenumberAtoms(mol, new_order)

        atoms = [atom for atom in mol.GetAtoms()]  # convert PyObject to Julia array
        atoms = sort!(atoms, by = x -> x.GetIdx()+1)  # sort by index to ensure same order as rdkit
        n_atom_feat = sum([length(atom_feat_bins[feat]) for feat in atom_feature_list])
        n_atoms = mol.GetNumAtoms()
        atom_feats = zeros(Float32,n_atoms,n_atom_feat)
        for (a_i,atom) in enumerate(atoms)
            atom_feats[a_i,:] = concat_smiles_atom_feat(smiles_atom_features(atom,feature_list=atom_feature_list))
        end

        bonds = [bond for bond in mol.GetBonds()]  # convert PyObject to Julia array
        bonds = sort!(bonds,by = x->(x.GetBeginAtomIdx()+1,x.GetEndAtomIdx()+1))
        bond_feats = OrderedDict()
        adj_list = [[] for n=1:n_atoms]
        for bond in bonds
            # sort index
            begin_atom_idx = bond.GetBeginAtomIdx()+1
            end_atom_idx = bond.GetEndAtomIdx()+1
            index = begin_atom_idx<end_atom_idx ? (begin_atom_idx,end_atom_idx) : (end_atom_idx, begin_atom_idx)
            bond_feats[index] = concat_smiles_bond_feat(smiles_bond_features(bond,feature_list=bond_feature_list))
            push!(adj_list[begin_atom_idx],end_atom_idx)
            push!(adj_list[end_atom_idx],begin_atom_idx)
        end
        # extract pair features
        bf_len = sum([length(bond_feat_bins[feat]) for feat in bond_feature_list])
        pair_feats = pair_features(mol,bond_feats,adj_list,bf_len = bf_len)
        outputs[i] = WeaveMol(atom_feats, pair_feats, mol.GetNumAtoms(), size(atom_feats,2))
    end
    return outputs
end


text:
@testset "pod of dolphins" begin
    # http://web.mit.edu/puzzle/www/2015/puzzle/pod_of_dolphins_meta/
    f = best_feature(["citygates", "impulsive", "clickspam", "baptistry", "leviathan", "policecar", "coupdetat", "sforzando", "cartwheel"])
    @test f.description == "has 8 unique letters"
end


text:
using RvSpectraKitLearn
using MultivariateStats

# Likely need to update this to point to directory containing the spectra you wish to analyze
include(joinpath(Pkg.dir("RvSpectraKitLearn"),"path_to_spectra.jl"))
cd(path_to_spectra)

(datafiles, phases) = get_filenames(); # num_max_files=4);
(lambda, obs_noisefree) = read_filelist(datafiles);

sampling = 3.0
planet_amplitude = 10.0
rvs_true = -planet_amplitude*cos.(2pi*phases)

snr_list = [ 300.0, 200.0, 150.0, 100.0, 80.0, 60.0, 40.0 ]
rms_simple_list = Array(Float64,0)
rms_gp_list = Array(Float64,0)
for snr in snr_list
 obs = make_noisy_spectra(obs_noisefree,snr,sampling=sampling)
 # First, use simple finite differences on averaged spectra to compute derivatives and resulting RVs
 doppler_comp_simple = calc_doppler_component_simple(lambda,obs);
 genpca_simple_out = fit_gen_pca_rv(obs,doppler_comp_simple)
 rvs_out_simple = est_rvs_from_pc(obs,genpca_simple_out[1],genpca_simple_out[2][:,1])
 # Now, use GP to compute derivatives and resulting RVs
 doppler_comp_gp = calc_doppler_component_gp(lambda,obs);
 genpca_gp_out = fit_gen_pca_rv(obs,doppler_comp_gp)
 # Evaluate quality of measured RVs
 rvs_out_gp = est_rvs_from_pc(obs,genpca_gp_out[1],genpca_gp_out[2][:,1])
 rms_simple = sqrt(sumabs2(rvs_out_simple.-rvs_true)/length(rvs_true))
 rms_gp = sqrt(sumabs2(rvs_out_gp.-rvs_true)/length(rvs_true))
 push!(rms_simple_list,rms_simple)
 push!(rms_gp_list,rms_gp)
 println("snr = ", snr, " rms = ", rms_simple, "  ", rms_gp)
end

using PyPlot
plt[:loglog](snr_list,rms_simple_list,"r.")
 plt[:loglog](snr_list,rms_gp_list,"g.")
 xlabel("SNR per Resolution Element")
 ylabel("RMS Velocity (m/s)")

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
using Libtask
using Test

# Test atomic values.
function f_cta()
  t = TRef(1);
  t[] = 0;
  while true
    produce(t[])
    t[]
    t[] += 1
  end
end

t = CTask(f_cta)

consume(t); consume(t)
a = copy(t);
consume(a); consume(a)

Base.@assert consume(t) == 2
Base.@assert consume(a) == 4

# Test dictionary functionality.
function dict_test()
  t = TRef(Dict("A" => 1, 5 => "B"));
  t["A"] = 0;
  while true
    produce(t["A"])
    t["A"]
    t["A"] += 1
  end
end

t = CTask(dict_test)

consume(t); consume(t)
a = copy(t);
consume(a); consume(a)

Base.@assert consume(t) == 2
Base.@assert consume(a) == 4

# Create a TRef storing a matrix.
x = TRef([1 2 3; 4 5 6])
x[1, 3] = 900
Base.@assert x[1,3] == 900

# TRef holding an array.
y = TRef([1,2,3])
y[2] = 19
Base.@assert y[2] == 19


text:
doc"""
    BetaBinomial(n,α,β)

A *Beta-binomial distribution* is the compound distribution of the [`Binomial`](:func:`Binomial`) distribution where the probability of success `p` is distributed according to the [`Beta`](:func:`Beta`). It has three parameters: `n`, the number of trials and two shape parameters `α`, `β`

$P(X = k) = {n \choose k} B(k + \alpha, n - k + \beta) / B(\alpha, \beta),  \quad \text{ for } k = 0,1,2, \ldots, n.$

```julia
BetaBinomial(n, a, b)      # BetaBinomial distribution with n trials and shape parameters a, b

params(d)       # Get the parameters, i.e. (n, a, b)
ntrials(d)      # Get the number of trials, i.e. n
```

External links:

* [Beta-binomial distribution on Wikipedia](https://en.wikipedia.org/wiki/Beta-binomial_distribution)
"""

immutable BetaBinomial{T<:Real} <: DiscreteUnivariateDistribution
    n::Int
    α::T
    β::T

    function (::Type{BetaBinomial{T}}){T}(n::Int, α::T, β::T)
        @check_args(BetaBinomial, n >= zero(n) && α >= zero(α) && β >= zero(β))
        new{T}(n, α, β)
    end
end

BetaBinomial{T<:Real}(n::Integer, α::T, β::T) = BetaBinomial{T}(n, α, β)
BetaBinomial(n::Integer, α::Real, β::Real) = BetaBinomial(n, promote(α, β)...)
BetaBinomial(n::Integer, α::Integer, β::Integer) = BetaBinomial(n, Float64(α), Float64(β))

@distr_support BetaBinomial 0 d.n
insupport(d::BetaBinomial, x::Real) = 0 <= x <= d.n

#### Conversions
function convert{T <: Real, S <: Real}(::Type{BetaBinomial{T}}, n::Int, α::S, β::S)
    BetaBinomial(n, T(α), T(β))
end
function convert{T <: Real, S <: Real}(::Type{BetaBinomial{T}}, d::BetaBinomial{S})
    BetaBinomial(d.n, T(d.α), T(d.β))
end

#### Parameters

ntrials(d::BetaBinomial) = d.n

params(d::BetaBinomial) = (d.n, d.α, d.β)
@inline partype{T<:Real}(d::BetaBinomial{T}) = T

#### Properties

mean(d::BetaBinomial) = (d.n * d.α) / (d.α + d.β)
function var(d::BetaBinomial)
    n, α, β = d.n, d.α, d.β
    numerator = n * α * β * (α + β + n)
    denominator = (α + β)^2 * (α + β + 1)
    return numerator / denominator
end

function skewness(d::BetaBinomial)
    n, α, β = d.n, d.α, d.β
    t1 = (α + β + 2n) * (β - α) / (α + β + 2)
    t2 = sqrt((1 + α +β) / (n * α * β * (n + α + β)))
    return t1 * t2
end

function kurtosis(d::BetaBinomial)
    n, α, β = d.n, d.α, d.β
    alpha_beta_sum = α + β
    alpha_beta_product = α * β
    numerator = ((alpha_beta_sum)^2) * (1 + alpha_beta_sum)
    denominator = (n * alpha_beta_product) * (alpha_beta_sum + 2) * (alpha_beta_sum + 3) * (alpha_beta_sum + n)
    left = numerator / denominator
    right = (alpha_beta_sum) * (alpha_beta_sum - 1 + 6n) + 3*alpha_beta_product * (n - 2) + 6n^2
    right -= (3*alpha_beta_product * n * (6 - n)) / alpha_beta_sum
    right -= (18*alpha_beta_product * n^2) / (alpha_beta_sum)^2
    return (left * right) - 3
end

function pdf(d::BetaBinomial, k::Int)
    n, α, β = d.n, d.α, d.β
    chooseinv = (n + 1) * beta(k + 1, n - k + 1)
    numerator = beta(k + α, n - k + β)
    denominator = beta(α, β)
    return numerator / (denominator * chooseinv)
end

function pdf(d::BetaBinomial)
    n, α, β = d.n, d.α, d.β
    denominator = beta(α, β)
    values = Array{eltype(denominator)}(n+1)
    for (i,k) in enumerate(0:n)
        chooseinv = (n + 1) * beta(k + 1, n - k + 1)
        numerator = beta(k + α, n - k + β)
        @inbounds values[i] = numerator / (denominator * chooseinv)
    end
    return values
end

entropy(d::BetaBinomial) = entropy(Categorical(pdf(d)))
median(d::BetaBinomial) = median(Categorical(pdf(d))) - 1
mode{T<:Real}(d::BetaBinomial{T}) = indmax(pdf(d)) - one(T)
modes(d::BetaBinomial) = [x - 1 for x in modes(Categorical(pdf(d)))]

quantile(d::BetaBinomial, p::Float64) = quantile(Categorical(pdf(d)), p) - 1


text:
module FuseApi

using CStructures: CStruct, CStructGuarded, CVector, CStructAccess, Cserialize
using CStructures: Layout, LFixedVector, LVarVector, LForwardReference

include("types.jl")
include("fuseapi.jl")
include("fusebridge.jl")
include("examplefs.jl")

end


text:
include("../functions/analysisPULSE.jl")

####################################################################
step_max_perturb_MTiXsens=step_max_perturb

Energysignal_MTiXsens,
rms_MTiXsens,
vdv_MTiXsens,
crest_factor_MTiXsens,
Energysignal_MTiXsens_polar,
rms_MTiXsens_polar,
vdv_MTiXsens_polar=analysisPULSE_MTiXsens(AccFree_MTiXsens_sync_perturb,
                                                    filenameMTiXsens[1:109],
                                                    sensorXsensNames,
                                                    Ts_MTiXsens,
                                                    step_max_perturb_MTiXsens)

Energysignal_MTiXsens_2,
rms_MTiXsens_2,
vdv_MTiXsens_2,
crest_factor_MTiXsens_2,
Energysignal_MTiXsens_polar_2,
rms_MTiXsens_polar_2,
vdv_MTiXsens_polar_2=analysisPULSE_MTiXsens(AccFree_MTiXsens_sync_perturb,
                                                filenameMTiXsens[90:end],
                                                sensorXsensNames,
                                                Ts_MTiXsens_40,
                                                step_max_perturb_MTiXsens)

merge!(Energysignal_MTiXsens,Energysignal_MTiXsens_2)
merge!(rms_MTiXsens,rms_MTiXsens_2)
merge!(vdv_MTiXsens,vdv_MTiXsens_2)
merge!(crest_factor_MTiXsens,crest_factor_MTiXsens_2)
merge!(Energysignal_MTiXsens_polar,Energysignal_MTiXsens_polar_2)
merge!(rms_MTiXsens_polar,rms_MTiXsens_polar_2)
merge!(vdv_MTiXsens_polar,vdv_MTiXsens_polar_2)

splitNameMTiXsens=[]
plot_order=[]
perturbLabel=[]
wheelType=[]
for k in filenameMTiXsens
        splitName=split(k,'_')
        append!(splitNameMTiXsens,[splitName])
        append!(plot_order,parse(Int,splitName[5][2]))
        append!(perturbLabel,[splitName[4]])
        append!(wheelType,[splitName[3]])
end
perturbLabel_=unique(perturbLabel)
perturbIndices=[]
for k in perturbLabel_
        append!(perturbIndices,findfirst(perturbLabel.==k))
end
for k in 1:size(perturbIndices,1)-1
        i1=perturbIndices[k]
        i2=perturbIndices[k+1]-1
        plot_order[i1:i2]=plot_order[i1:i2].+(0.05*(k-1))
end
plot_order[perturbIndices[end]:end]=plot_order[perturbIndices[end]:end].+(0.05*(size(perturbIndices,1)-1))

# include("../functions/analysisPULSE.jl")
plotAnalysisPULSE_MTiXsens(
        Energysignal_MTiXsens,
        rms_MTiXsens,
        vdv_MTiXsens,
        crest_factor_MTiXsens,
        filenameMTiXsens,
        sensorXsensNames,
        ["x-y";"z"],
        perturbLabel_,
        plot_order,
        perturbIndices,
        [:blue :red :green :orange :yellow],
        [:rect :rect :rect :rect :circle])

plotAnalysisPULSE_MTiXsens(
        Energysignal_MTiXsens_polar,
        rms_MTiXsens_polar,
        vdv_MTiXsens_polar,
        crest_factor_MTiXsens,
        filenameMTiXsens,
        sensorXsensNames,
        ["r";"\u03B8"],
        perturbLabel_,
        plot_order,
        perturbIndices,
        [:blue :red :green :orange :yellow],
        [:rect :rect :rect :rect :circle])


text:
origin(prim::HyperRectangle) = prim.origin
maximum(prim::HyperRectangle) = origin(prim) + widths(prim)
minimum(prim::HyperRectangle) = origin(prim)
length(prim::HyperRectangle{N, T}) where {T, N} = N
widths(prim::HyperRectangle) = prim.widths

"""
Splits an HyperRectangle into two along an axis at a given location.
"""
split(b::HyperRectangle, axis, value::Integer) = _split(b, axis, value)
split(b::HyperRectangle, axis, value::Number) = _split(b, axis, value)
function _split(b::H, axis, value) where H<:HyperRectangle
    bmin = minimum(b)
    bmax = maximum(b)
    b1max = setindex(bmax, value, axis)
    b2min = setindex(bmin, value, axis)

    return H(bmin, b1max-bmin),
           H(b2min, bmax-b2min)
end

# empty constructor such that update will always include the first point
function (HR::Type{HyperRectangle{N,T}})() where {T,N}
    HR(Vec{N,T}(typemax(T)), Vec{N,T}(typemin(T)))
end

# conversion from other HyperRectangles
function (HR::Type{HyperRectangle{N,T1}})(a::HyperRectangle{N,T2}) where {N,T1,T2}
    HR(Vec{N, T1}(minimum(a)), Vec{N, T1}(widths(a)))
end

function HyperRectangle(v1::Vec{N,T1}, v2::Vec{N,T2}) where {N,T1,T2}
    T = promote_type(T1,T2)
    HyperRectangle{N,T}(Vec{N,T}(v1), Vec{N,T}(v2))
end

function (HR::Type{HyperRectangle{N,T}})(a::GeometryPrimitive) where {N,T}
    HR(Vec{N, T}(minimum(a)), Vec{N, T}(widths(a)))
end
"""
```
HyperRectangle(vals::Number...)
```
HyperRectangle constructor for indidually specified intervals.
e.g. HyperRectangle(0,0,1,2) has origin == Vec(0,0) and
width == Vec(1,2)
"""
@generated function HyperRectangle(vals::Number...)
    # Generated so we get goodish codegen on each signature
    n = length(vals)
    @assert iseven(n)
    mid = div(n,2)
    v1 = Expr(:call, :Vec)
    v2 = Expr(:call, :Vec)
    # TODO this can be inbounds
    append!(v1.args, [:(vals[$i]) for i = 1:mid])
    append!(v2.args, [:(vals[$i]) for i = mid+1:length(vals)])
    Expr(:call, :HyperRectangle, v1, v2)
end

function HyperRectangle(r::SimpleRectangle{T}) where T
    HyperRectangle{2,T}(r)
end

function HyperRectangle{N,T}(r::SimpleRectangle) where {N,T}
    if N > 2
        return HyperRectangle(Vec{N, T}(T(r.x), T(r.y), Vec{N-2,T}(zero(T))...),
                              Vec{N, T}(T(r.w), T(r.h), Vec{N-2,T}(zero(T))...))
    else
        return HyperRectangle(Vec{N, T}(T(r.x), T(r.y)),
                              Vec{N, T}(T(r.w), T(r.h)))
    end
end

"""
Transform a `HyperRectangle` using a matrix. Maintains axis-align properties
so a significantly larger HyperRectangle may be generated.
"""
function *(m::Mat{N1,N1,T1}, h::HyperRectangle{N2, T2}) where {N1,N2,T1,T2}

    # TypeVar constants
    T = promote_type(T1, T2)
    D = N1 - N2

    # get all points on the HyperRectangle
    d = decompose(Point, h)
    # make sure our points are sized for the tranform
    pts = (Vec{N1, T}[vcat(pt, ones(Vec{D, T})) for pt in d]...,)::NTuple{2^N2,Vec{N1,T}}

    vmin = Vec{N1, T}(typemax(T))
    vmax = Vec{N1, T}(typemin(T))
    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * pt
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    HyperRectangle{N2, T}(vmin, vmax - vmin)
end

function -(h::HyperRectangle{N, T}, move::Vec{N}) where {N,T}
    HyperRectangle{N, T}(minimum(h) .- move, widths(h))
end
function +(h::HyperRectangle{N, T}, move::Vec{N}) where {N,T}
    HyperRectangle{N, T}(minimum(h) .+ move, widths(h))
end

function -(h::SimpleRectangle{T}, move::Vec{2}) where T
    SimpleRectangle{T}((minimum(h) .- move)..., widths(h)...)
end
function +(h::SimpleRectangle{T}, move::Vec{2}) where T
    SimpleRectangle{T}((minimum(h) .+ move)..., widths(h)...)
end

function *(m::Mat{N,N,T1}, h::HyperRectangle{N,T2}) where {N,T1,T2}
    # equal dimension case

    # TypeVar constants
    T = promote_type(T1,T2)

    # get all points on the HyperRectangle
    pts = decompose(Point, h)

    # make sure our points are sized for the tranform
    vmin = Vec{N, T}(typemax(T))
    vmax = Vec{N, T}(typemin(T))

    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * Vec(pt)
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    HyperRectangle{N,T}(vmin, vmax-vmin)
end

# fast path. TODO make other versions fast without code duplications like now
function *(m::Mat{4,4,T}, h::HyperRectangle{3,T}) where T
    # equal dimension case

    # get all points on the HyperRectangle
    pts = (
        Vec{4,T}(0.0,0.0,0.0, 1.0),
        Vec{4,T}(1.0,0.0,0.0, 1.0),
        Vec{4,T}(0.0,1.0,0.0, 1.0),
        Vec{4,T}(1.0,1.0,0.0, 1.0),
        Vec{4,T}(0.0,0.0,1.0, 1.0),
        Vec{4,T}(1.0,0.0,1.0, 1.0),
        Vec{4,T}(0.0,1.0,1.0, 1.0),
        Vec{4,T}(1.0,1.0,1.0, 1.0)
    )

    # make sure our points are sized for the tranform
    vmin = Vec{4, T}(typemax(T))
    vmax = Vec{4, T}(typemin(T))
    o, w = origin(h), widths(h)
    _o = Vec{4, T}(o[1], o[2], o[3], T(0))
    _w = Vec{4, T}(w[1], w[2], w[3], T(1))
    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * (_o + (pt .* _w))
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    _vmin = Vec{3, T}(vmin[1], vmin[2], vmin[3])
    _vmax = Vec{3, T}(vmax[1], vmax[2], vmax[3])
    HyperRectangle{3,T}(_vmin, _vmax - _vmin)
end

function HyperRectangle(geometry::AbstractArray{<: Point{N, T}}) where {N,T}
    HyperRectangle{N,T}(geometry)
end

@inline function minmax(p::StaticVector, vmin, vmax)
    isnan(p) && return (vmin, vmax)
    min.(p, vmin), max.(p, vmax)
end

# Annoying special case for view(Vector{Point}, Vector{Face})
@inline function minmax(tup::Tuple, vmin, vmax)
    for p in tup
        isnan(p) && continue
        vmin = min.(p, vmin)
        vmax = max.(p, vmax)
    end
    vmin, vmax
end

"""
Construct a HyperRectangle enclosing all points.
"""
function (t::Type{HyperRectangle{N1, T1}})(
        geometry::AbstractArray{PT}
    ) where {N1, T1, PT <: Point}
    N2, T2 = length(PT), eltype(PT)
    @assert N1 >= N2
    vmin = Point{N2, T2}(typemax(T2))
    vmax = Point{N2, T2}(typemin(T2))
    for p in geometry
        vmin, vmax = minmax(p, vmin, vmax)
    end
    o = vmin
    w = vmax - vmin
    if N1 > N2
        z = zero(Vec{N1-N2, T1})
        return HyperRectangle{N1, T1}(vcat(o, z),
                                     vcat(w, z))
    else
        return HyperRectangle{N1, T1}(o, w)
   end
end

xwidth(a::SimpleRectangle)  = a.w + a.x
width(a::SimpleRectangle)  = a.w
yheight(a::SimpleRectangle) = a.h + a.y
height(a::SimpleRectangle) = a.h
area(a::SimpleRectangle) = a.w*a.h
widths(a::SimpleRectangle{T}) where {T} = Point{2,T}(a.w, a.h)
maximum(a::SimpleRectangle{T}) where {T} = Point{2, T}(a.x + widths(a)[1], a.y +widths(a)[2])
minimum(a::SimpleRectangle{T}) where {T} = Point{2, T}(a.x, a.y)
origin(a::SimpleRectangle{T}) where {T} = Point{2, T}(a.x, a.y)

SimpleRectangle(val::Vec{2, T}) where {T} = SimpleRectangle{T}(0, 0, val...)
function SimpleRectangle(position::Vec{2,T}, width::Vec{2,T}) where T
    SimpleRectangle{T}(position..., width...)
end

function Base.to_indices(A::AbstractArray{T, 2}, I::Tuple{<: SimpleRectangle}) where T
    i = I[1]
    (i.x + 1 : (i.x + i.w), i.y + 1 : (i.y + i.h))
end

AbsoluteRectangle(mini::Vec{N,T}, maxi::Vec{N,T}) where {N,T} = HyperRectangle{N,T}(mini, maxi-mini)

AABB(a) = AABB{Float32}(a)

function (B::Type{AABB{T}})(a::Pyramid) where T
    w,h = a.width/T(2), a.length
    m = Vec{3,T}(a.middle)
    B(m-Vec{3,T}(w,w,0), m+Vec{3,T}(w, w, h))
end

(B::Type{AABB{T}})(a::Cube) where {T} = B(origin(a), widths(a))

(B::Type{AABB{T}})(a::AbstractMesh) where {T} = B(vertices(a))


text:
#=
run polynomial utilities tests
=#

using Test
using Printf
include(joinpath(@__DIR__, "polyutils.jl"))

real_types = [
    Float64,
    Float32,
    BigFloat,
    ]

@testset "polynomial utilities tests" begin
@testset "$T" for T in real_types
    println("$T ...")
    test_time = @elapsed begin
        test_interp_domain(T)
        test_complex_interp(T)
        test_quadrature(T)
    end
    @printf("%8.2e seconds\n", test_time)
end
end
;


text:
module MeltingParam

# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam, PhaseDiagram_LookupTable, AbstractMaterialParamsStruct
import Base.show, GeoParams.param_info
using ..MaterialParameters: MaterialParamsInfo

abstract type AbstractMeltingParam{T} <: AbstractMaterialParam end

export  compute_meltfraction,
        compute_meltfraction!,    # calculation routines
        param_info,
        MeltingParam_Caricchi     # constant

# Constant  -------------------------------------------------------
"""
    MeltingParam_Caricchi()

Implements the T-dependent melting parameterisation used by Caricchi et al
```math
    \\theta = (800.0 .- (T + 273.15))./23.0
```
```math
    \\phi_{solid} = 1.0 - {1.0 \\over (1.0 + e^\\theta)};
```

Note that T is in Kelvin.

"""
@with_kw_noshow struct MeltingParam_Caricchi{T,U} <: AbstractMeltingParam{T}
    a::GeoUnit{T,U}              =   800.0K
    b::GeoUnit{T,U}              =   23.0K
    c::GeoUnit{T,U}              =   273.15K # shift from C to K
end
MeltingParam_Caricchi(args...) = MeltingParam_Caricchi(convert.(GeoUnit,args)...)

function param_info(s::MeltingParam_Caricchi) # info about the struct
    return MaterialParamsInfo(Equation =  L"\phi = {1 \over {1 + \exp( {800-T[^oC] \over 23})}}")
end

# Calculation routine
function compute_meltfraction(p::MeltingParam_Caricchi{_T}, P::Quantity, T::Quantity) where _T
    a,b,c   = p.a, p.b, p.c

    θ       =   (a - (T - c))/b
    ϕ       =   1.0./(1.0 .+ exp.(θ))

    return ϕ
end

function compute_meltfraction!(ϕ, p::MeltingParam_Caricchi{_T}, P::_T, T::_T ) where _T
    @unpack_val a,b,c   = p

    θ       =   (a - (T - c))/b
    ϕ       =   1.0/(1.0 + exp(θ))

    return nothing
end

function compute_meltfraction!(ϕ::AbstractArray{_T}, p::MeltingParam_Caricchi{_T}, P::AbstractArray{_T},T::AbstractArray{_T}) where _T
    @unpack_val a,b,c   = p

    θ       =   (a .- (T .- c))./b
    ϕ      .=   1.0./(1.0 .+ exp.(θ))

    return nothing
end

# Print info
function show(io::IO, g::MeltingParam_Caricchi)
    print(io, "Caricchi et al. melting parameterization")
end
#-------------------------------------------------------------------------

"""
    ComputeMeltingParam(P,T, p::AbstractPhaseDiagramsStruct)

Computes melt fraction in case we use a phase diagram lookup table. The table should have the collum `:meltFrac` specified.
"""
function compute_meltfraction(p::PhaseDiagram_LookupTable, P::_T,T::_T) where _T
   return p.meltFrac.(T,P)
end

"""
    ComputeMeltingParam!(ϕ::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T:AbstractArray{<:AbstractFloat}, p::PhaseDiagram_LookupTable)

In-place computation of melt fraction in case we use a phase diagram lookup table. The table should have the collum `:meltFrac` specified.
"""
function compute_meltfraction!(ϕ::AbstractArray{_T}, p::PhaseDiagram_LookupTable, P::AbstractArray{_T}, T::AbstractArray{_T}) where _T
    ϕ[:]    =   p.meltFrac.(T,P)

    return nothing
end

"""
    ComputeMeltingParam!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
"""
function compute_meltfraction!(ϕ::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}, Phases::AbstractArray{<:Integer, N}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}) where N

    for i = 1:length(MatParam)

        if length(MatParam[i].Melting)>0

            # Create views into arrays (so we don't have to allocate)
            ind = Phases .== MatParam[i].Phase;
            ϕ_local   =   view(ϕ, ind )
            P_local   =   view(P  , ind )
            T_local   =   view(T  , ind )

            compute_meltfraction!(ϕ_local, MatParam[i].Melting[1], P_local, T_local)
        end

    end

    return nothing
end

"""
    ComputeMeltingParam!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
"""
function compute_meltfraction!(ϕ::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}) where {N,M}

    ϕ .= 0.0
    for i = 1:length(MatParam)

        ϕ_local  = zeros(size(ϕ))
        Fraction    = selectdim(PhaseRatios,M,i);
        if (maximum(Fraction)>0.0) & (length(MatParam[i].Melting)>0)

            compute_meltfraction!(ϕ_local, MatParam[i].Melting[1] , P, T)

            ϕ .= ϕ .+ ϕ_local.*Fraction
        end

    end

    return nothing
end

end