text:
#@static Sys.iswindows() ? (const thelib = "c:/v/build5/src/gmt_w64") : (const thelib = "libgmt")  # Name of GMT shared lib.
#@static Sys.iswindows() ? (const thelib = "C:/progs_cygw/GMTdev/gmt5/5.4/WIN64/bin/gmt_w64") : (const thelib = "libgmt")  # Name of GMT shared lib.
@static Sys.iswindows() ? (Sys.WORD_SIZE == 64 ? (const thelib = "gmt_w64") : (const thelib = "gmt_w32")) : (const thelib = "libgmt")  # Name of GMT shared lib.

function GMT_Create_Session(tag::String="GMT", pad=2, mode=0, print_func::Ptr{Cvoid}=C_NULL)
	ccall( (:GMT_Create_Session, thelib), Ptr{Cvoid}, (Ptr{UInt8}, UInt32, UInt32, Ptr{Cvoid}), tag, pad, mode, print_func)
end

function GMT_Create_Data(API::Ptr{Cvoid}, family, geometry, mode, dim=C_NULL, wesn=C_NULL,
                         inc=C_NULL, registration=0, pad=2, data::Ptr{Cvoid}=C_NULL)

	if (family == GMT_IS_DATASET)        ret_type = Ptr{GMT_DATASET}
	elseif (family == GMT_IS_TEXTSET)    ret_type = Ptr{GMT_TEXTSET}
	elseif (family == GMT_IS_GRID)       ret_type = Ptr{GMT_GRID}
	elseif (family == GMT_IS_CPT)        ret_type = Ptr{GMT_PALETTE}
	elseif (family == GMT_IS_IMAGE)      ret_type = Ptr{GMT_IMAGE}
	elseif (family == GMT_IS_MATRIX)     ret_type = Ptr{GMT_MATRIX}
	elseif (family == GMT_IS_MATRIX|GMT_VIA_MATRIX) ret_type = Ptr{GMT_MATRIX}
	elseif (family == GMT_IS_VECTOR)     ret_type = Ptr{GMT_VECTOR}
	elseif (family == GMT_IS_POSTSCRIPT) ret_type = Ptr{GMT_POSTSCRIPT}
	else                                 ret_type = Ptr{Cvoid}		# Should be error instead
	end

	ptr = ccall((:GMT_Create_Data, thelib), Ptr{Cvoid}, (Cstring, UInt32, UInt32, UInt32, Ptr{UInt64},
		Ptr{Cdouble}, Ptr{Cdouble}, UInt32, Cint, Ptr{Cvoid}), API, family, geometry, mode, dim, wesn, inc,
		registration, pad, data)

	convert(ret_type, ptr)
end

function GMT_Destroy_Session(API::Ptr{Cvoid})
	ccall( (:GMT_Destroy_Session, thelib), Cint, (Cstring,), API)
end

#=		Not used yet, so comment
function GMT_Read_Data(API::Ptr{Cvoid}, family, method, geometry, mode, wesn, input=C_NULL, data=C_NULL)

	if (family == GMT_IS_DATASET)
		ret_type = Ptr{GMT_DATASET}
	elseif (family == GMT_IS_TEXTSET)
		ret_type = Ptr{GMT_TEXTSET}
	elseif (family == GMT_IS_GRID)
		ret_type = Ptr{GMT_GRID}
	elseif (family == GMT_IS_CPT)
		ret_type = Ptr{GMT_PALETTE}
	elseif (family == GMT_IS_IMAGE)
		ret_type = Ptr{GMT_IMAGE}
	elseif (family == GMT_IS_MATRIX)
		ret_type = Ptr{GMT_MATRIX}
	elseif (family == GMT_IS_VECTOR)
		ret_type = Ptr{GMT_VECTOR}
	else
		ret_type = Ptr{Cvoid}			# Should be error instead
	end

	ptr = ccall( (:GMT_Read_Data, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
		Ptr{UInt8}, Ptr{Cvoid}), API, family, method, geometry, mode, wesn, input, data)

	convert(ret_type, ptr)
end

function GMT_Insert_Data(API::Ptr{Cvoid}, object_ID::Integer, data)
	ccall((:GMT_Insert_Data, thelib), Cint, (Ptr{Cvoid}, Cint, Ptr{Cvoid}), API, object_ID, data)
end

function GMT_Duplicate_Data(API::Ptr{Cvoid}, family::Integer, mode::Integer, data::Ptr{Cvoid})
	ccall( (:GMT_Duplicate_Data, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, family, mode, data)
end

function GMT_Get_Record(API::Ptr{Cvoid}, mode::Integer, retval::Ptr{Cint})
	ccall( (:GMT_Get_Record, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, Ptr{Cint}), API, mode, retval)
end

function GMT_Manage_Session(API::Ptr{Cvoid}, mode::Integer, args::Ptr{Cvoid})
	ccall( (:GMT_Manage_Session, thelib), Cint, (Ptr{Cvoid}, UInt32, Ptr{Cvoid}), API, mode, args)
end

function GMT_Register_IO(API, family::Integer, method::Integer, geometry::Integer, direction::Integer=0,
		wesn=C_NULL, resource=C_NULL)
	err = ccall((:GMT_Register_IO, thelib), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
		Ptr{Cvoid}), API, family, method, geometry, direction, wesn, resource)
end

function GMT_Init_IO(API::Ptr{Cvoid}, family::UInt32, geometry::UInt32, direction::UInt32, mode::UInt32, n_args::UInt32, args::Ptr{Cvoid})
	ccall( (:GMT_Init_IO, thelib), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, UInt32, Ptr{Cvoid}), API, family, geometry, direction, mode, n_args, args)
end
function GMT_Begin_IO(API::Ptr{Cvoid}, family::UInt32, direction::UInt32, header::UInt32)
	ccall( (:GMT_Begin_IO, thelib), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32), API, family, direction, header)
end
function GMT_Status_IO(API::Ptr{Cvoid}, mode::UInt32)
	ccall( (:GMT_Status_IO, thelib), Cint, (Ptr{Cvoid}, UInt32), API, mode)
end
function GMT_End_IO(API::Ptr{Cvoid}, direction::UInt32, mode::UInt32)
	ccall( (:GMT_End_IO, thelib), Cint, (Ptr{Cvoid}, UInt32, UInt32), API, direction, mode)
end

function GMT_Write_Data(API::Ptr{Cvoid}, family::Integer, method::Integer, geometry::Integer, mode::Integer,
	wesn, output::String, data)
	err = ccall((:GMT_Write_Data, thelib), Cint, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
		Ptr{UInt8}, Ptr{Cvoid}), API, family, method, geometry, mode, wesn, output, data)
end
=#

function GMT_Destroy_Data(API::Ptr{Cvoid}, object)
	ccall( (:GMT_Destroy_Data, thelib), Cint, (Cstring, Ptr{Cvoid}), API, object)
end

function GMT_Set_Comment(API::Ptr{Cvoid}, family::Integer, mode::Integer, arg::Ptr{Cvoid}, data::Ptr{Cvoid})
	ccall( (:GMT_Set_Comment, thelib), Cint, (Cstring, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}), API, family, mode, arg, data)
end

#=
function GMT_Put_Record(API::Ptr{Cvoid}, mode::UInt32, record::Ptr{Cvoid})
	ccall( (:GMT_Put_Record, thelib), Cint, (Ptr{Cvoid}, UInt32, Ptr{Cvoid}), API, mode, record)
end

function GMT_Encode_ID(API::Ptr{Cvoid}, fname::String, object_ID::Integer)
	err = ccall((:GMT_Encode_ID, thelib), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint), API, fname, object_ID)
end

function GMT_Get_Row(API::Ptr{Cvoid}, rec_no::Cint, G::Ptr{GMT_GRID}, row::Ptr{Cfloat})
	ccall( (:GMT_Get_Row, thelib), Cint, (Ptr{Cvoid}, Cint, Ptr{GMT_GRID}, Ptr{Cfloat}), API, rec_no, G, row)
end
function GMT_Put_Row(API::Ptr{Cvoid}, rec_no::Cint, G::Ptr{GMT_GRID}, row::Ptr{Cfloat})
	ccall( (:GMT_Put_Row, thelib), Cint, (Ptr{Cvoid}, Cint, Ptr{GMT_GRID}, Ptr{Cfloat}), API, rec_no, G, row)
end

function GMT_Get_ID(API::Ptr{Cvoid}, family::Integer, dir::Integer, resource=C_NULL)
	ccall((:GMT_Get_ID, thelib), Cint, (Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, family, dir, resource)
end

function GMT_Get_Family(API::Ptr{Cvoid}, dir::Integer, head::Ptr{GMT_OPTION})
	ccall((:GMT_Get_Family, thelib), Cint, (Ptr{Cvoid}, UInt32, Ptr{GMT_OPTION}), API, dir, head)
end

function GMT_Get_Index(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, row::Cint, col::Cint)
	ccall((:GMT_Get_Index, thelib), Clonglong, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Cint, Cint), API, header, row, col)
end
function GMT_Get_Coord(API::Ptr{Cvoid}, family::Integer, dim::Integer, container::Ptr{Cvoid})
	ccall((:GMT_Get_Coord, thelib), Ptr{Cdouble}, (Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, family, dim, container)
end

function GMT_Option(API::Ptr{Cvoid}, options)
	ccall((:GMT_Option, thelib), Cint, (Ptr{Cvoid}, Ptr{UInt8}), API, options)
end

function GMT_Get_Common(API::Ptr{Cvoid}, option::UInt32, par::Ptr{Cdouble})
	ccall((:GMT_Get_Common, thelib), Cint, (Ptr{Cvoid}, UInt32, Ptr{Cdouble}), API, option, par)
end

function GMT_Get_Value(API::Ptr{Cvoid}, arg::String, par::Ptr{Cdouble})
	ccall((:GMT_Get_Value, thelib), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Cdouble}), API, arg, par)
end
function GMT_Get_Values(API::Ptr{Cvoid}, arg::String, par::Ptr{Cdouble}, maxpar::Integer)
	ccall((:GMT_Get_Values, thelib), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Cdouble}, Cint), API, arg, par, maxpar)
end
=#

function GMT_Get_Default(API::Ptr{Cvoid}, keyword::String, value)
    ccall((:GMT_Get_Default, thelib), Cint, (Cstring, Ptr{UInt8}, Ptr{UInt8}), API, keyword, value)
end

function GMT_Call_Module(API::Ptr{Cvoid}, _module=C_NULL, mode::Integer=0, args=C_NULL)
	if (isa(args,String))	args = pointer(args)	end
	ccall((:GMT_Call_Module, thelib), Cint, (Cstring, Ptr{UInt8}, Cint, Ptr{Cvoid}), API, _module, mode, args)
end

function GMT_Create_Options(API::Ptr{Cvoid}, argc::Integer, args)
	# VERSATILIZAR PARA O CASO DE ARGS SER STRING OU ARRAY DE STRINGS
	#ccall((:GMT_Create_Options, thelib), Ptr{GMT_OPTION}, (Ptr{Cvoid}, Cint, Ptr{Cvoid}), API, argc, args)
	ccall((:GMT_Create_Options, thelib), Ptr{GMT_OPTION}, (Cstring, Cint, Cstring), API, argc, args)
end
#GMT_Create_Options(API::Ptr{Cvoid}, argc::Integer, args::String) =
#                   GMT_Create_Options(API, argc, convert(Ptr{Cvoid},pointer(args)))

function GMT_Destroy_Options(API::Ptr{Cvoid}, head::Ref{Ptr{GMT_OPTION}})
	ccall( (:GMT_Destroy_Options, thelib), Cint, (Cstring, Ref{Ptr{GMT_OPTION}}), API, head)
end

#=
function GMT_Make_Option(API::Ptr{Cvoid}, option::UInt8, arg::Ptr{UInt8})
	ccall((:GMT_Make_Option, thelib), Ptr{GMT_OPTION}, (Ptr{Cvoid}, UInt8, Ptr{UInt8}), API, option, arg)
end

function GMT_Find_Option(API::Ptr{Cvoid}, option::UInt8, head::Ptr{GMT_OPTION})
	ccall((:GMT_Find_Option, thelib), Ptr{GMT_OPTION}, (Ptr{Cvoid}, UInt8, Ptr{GMT_OPTION}), API, option, head)
end

function GMT_Append_Option(API::Ptr{Cvoid}, current::Ptr{GMT_OPTION}, head::Ptr{GMT_OPTION})
	ccall((:GMT_Append_Option, thelib), Ptr{GMT_OPTION}, (Ptr{Cvoid}, Ptr{GMT_OPTION}, Ptr{GMT_OPTION}), API, current, head)
end
function GMT_Create_Args(API::Ptr{Cvoid}, argc::Ptr{Int}, head::Ptr{GMT_OPTION})
	ccall( (:GMT_Create_Args, thelib), Ptr{Ptr{UInt8}}, (Ptr{Cvoid}, Ptr{Cint}, Ptr{GMT_OPTION}), API, argc, head)
end
function GMT_Create_Cmd(API::Ptr{Cvoid}, head::Ptr{GMT_OPTION})
	ccall( (:GMT_Create_Cmd, thelib), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{GMT_OPTION}), API, head)
end

function GMT_Destroy_Args(API::Ptr{Cvoid}, argc::Cint, argv::Ptr{Ptr{Ptr{UInt8}}})
	ccall( (:GMT_Destroy_Args, thelib), Cint, (Ptr{Cvoid}, Cint, Ptr{Ptr{Ptr{UInt8}}}), API, argc, argv)
end
function GMT_Destroy_Cmd(API::Ptr{Cvoid}, cmd::Ptr{Ptr{UInt8}})
	ccall( (:GMT_Destroy_Cmd, thelib), Cint, (Ptr{Cvoid}, Ptr{Ptr{UInt8}}), API, cmd)
end
function GMT_Update_Option(API::Ptr{Cvoid}, current::Ptr{GMT_OPTION}, arg::Ptr{UInt8})
	ccall( (:GMT_Update_Option, thelib), Cint, (Ptr{Cvoid}, Ptr{GMT_OPTION}, Ptr{UInt8}), API, current, arg)
end
function GMT_Delete_Option(API::Ptr{Cvoid}, current::Ptr{GMT_OPTION}, head::Ptr{Ptr{GMT_OPTION}})
	ccall( (:GMT_Delete_Option, thelib), Cint, (Ptr{Cvoid}, Ptr{GMT_OPTION}, Ptr{Ptr{GMT_OPTION}}), API, current, head)
end
function GMT_Parse_Common(API::Ptr{Cvoid}, given_options::Ptr{UInt8}, options::Ptr{GMT_OPTION})
	ccall( (:GMT_Parse_Common, thelib), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{GMT_OPTION}), API, given_options, options)
end

function GMT_FFT_Option(API::Ptr{Cvoid}, option::UInt8, dim::UInt32, string::Ptr{UInt8})
	ccall( (:GMT_FFT_Option, thelib), UInt32, (Ptr{Cvoid}, UInt8, UInt32, Ptr{UInt8}), API, option, dim, string)
end
function GMT_FFT_Parse(API::Ptr{Cvoid}, option::UInt8, dim::UInt32, args::Ptr{UInt8})
	ccall( (:GMT_FFT_Parse, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, UInt8, UInt32, Ptr{UInt8}), API, option, dim, args)
end
function GMT_FFT_Create(API::Ptr{Cvoid}, X::Ptr{Cvoid}, dim::UInt32, mode::UInt32, F::Ptr{Cvoid})
	ccall( (:GMT_FFT_Create, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt32, UInt32, Ptr{Cvoid}), API, X, dim, mode, F)
end
function GMT_FFT_Wavenumber(API::Ptr{Cvoid}, k::UInt16, mode::UInt32, K::Ptr{Cvoid})
	ccall( (:GMT_FFT_Wavenumber, thelib), Cdouble, (Ptr{Cvoid}, UInt16, UInt32, Ptr{Cvoid}), API, k, mode, K)
end
function GMT_FFT(API::Ptr{Cvoid}, X::Ptr{Cvoid}, direction::Cint, mode::UInt32, K::Ptr{Cvoid})
	ccall( (:GMT_FFT, thelib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, UInt32, Ptr{Cvoid}), API, X, direction, mode, K)
end
function GMT_FFT_Destroy(API::Ptr{Cvoid}, K::Ptr{Cvoid})
	ccall( (:GMT_FFT_Destroy, thelib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), API, K)
end
function GMT_FFT_1D(API::Ptr{Cvoid}, data::Ptr{Cfloat}, n::UInt16, direction::Cint, mode::UInt32)
	ccall( (:GMT_FFT_1D, thelib), Cint, (Ptr{Cvoid}, Ptr{Cfloat}, UInt16, Cint, UInt32), API, data, n, direction, mode)
end
function GMT_FFT_2D(API::Ptr{Cvoid}, data::Ptr{Cfloat}, nx::UInt32, ny::UInt32, direction::Cint, mode::UInt32)
	ccall( (:GMT_FFT_2D, thelib), Cint, (Ptr{Cvoid}, Ptr{Cfloat}, UInt32, UInt32, Cint, UInt32), API, data, nx, ny, direction, mode)
end
=#

function GMT_Report(API, vlevel::Int, txt)
	ccall((:GMT_Report, thelib), Cvoid, (Cstring, Cint, Ptr{UInt8}), API, vlevel, txt)
end

function GMT_Encode_Options(V_API::Ptr{Cvoid}, _module, n_argin::Integer, head::Ref{Ptr{GMT_OPTION}}, n::Ptr{Int})
	ccall((:GMT_Encode_Options, thelib), Ptr{GMT_RESOURCE}, (Cstring, Ptr{UInt8}, Int32, Ref{Ptr{GMT_OPTION}},
					Ptr{UInt32}), V_API, _module, n_argin, head, n)
end
function GMT_Encode_Options(V_API::Ptr{Cvoid}, _module, n_argin::Integer, head, n::Ptr{Int})
	ccall((:GMT_Encode_Options, thelib), Ptr{GMT_RESOURCE}, (Cstring, Ptr{UInt8}, Int32, Ptr{Ptr{Cvoid}}, Ptr{UInt32}),
					V_API, _module, n_argin, head, n)
end

function GMT_Expand_Option(V_API::Ptr{Cvoid}, opt::Ptr{GMT_OPTION}, arg)
	ccall((:GMT_Expand_Option, thelib), Cint, (Cstring, Ptr{GMT_OPTION}, Ptr{UInt8}), V_API, opt, arg)
end

#=
function gmt_core_module_info(API, candidate)
	ccall((:gmt_core_module_info, thelib), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{UInt8}), API, candidate)
end

function gmtlib_grd_flip_vertical(gridp, n_cols::Integer, n_rows::Integer, n_stride::Integer=0, cell_size::Integer=1)
	ccall((:gmtlib_grd_flip_vertical, thelib), Cvoid, (Ptr{Cvoid}, UInt32, UInt32, UInt32, Csize_t),
				 gridp, n_cols, n_rows, n_stride, cell_size)
end

function GMT_set_mem_layout(API, mem_layout)
	ccall((:GMT_set_mem_layout, thelib), Cvoid, (Ptr{Cvoid}, Ptr{UInt8}), API, mem_layout)
end
=#

function GMT_Set_Default(API::Ptr{Cvoid}, keyword, value)
	ccall((:GMT_Set_Default, thelib), Cvoid, (Cstring, Ptr{UInt8}, Ptr{UInt8}), API, keyword, value)
end

function GMT_blind_change_struct(API::Ptr{Cvoid}, X, what, keyword::String, off::Integer)
	ccall((:GMT_blind_change_struct, thelib), Cint, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}, Csize_t),
				 API, X, what, keyword, off)
end

function GMT_Convert_Data(API::Ptr{Cvoid}, In::Ptr{Cvoid}, family_in::Integer, out::Ptr{Cvoid}, family_out::Integer, flag)
	ccall((:GMT_Convert_Data, thelib), Ptr{Cvoid}, (Cstring, Ptr{Cvoid}, UInt32, Ptr{Cvoid}, UInt32, Ptr{UInt32}), API, In,
				 family_in, out, family_out, flag)
end

#=
function GMT_blind_change_struct_(API::Ptr{Cvoid}, X, what, keyword::String)
	ccall((:GMT_blind_change_struct_, thelib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{UInt8}),
				 API, X, what, keyword)
end

function GMT_Read_Group(API::Ptr{Cvoid}, family::Integer, method::Integer, geometry::Integer, mode::Integer,
												wesn, sources::Ptr{Cvoid}, n_items::Ptr{UInt32}, data::Ptr{Cvoid})
	ccall((:GMT_Read_Group, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, UInt32, UInt32, UInt32, Ptr{Cdouble},
				 Ptr{Cvoid}, Ptr{UInt32}, Ptr{Cvoid}), API, family, method, geometry, mode, wesn, sources, n_items, data)
end

function GMT_Get_Index(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, row::Integer, col::Integer)
	ccall((:GMT_Get_Index, thelib), UInt64, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Cint, Cint), API, header, row, col)
end

function GMT_Get_Pixel(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, row::Integer, col::Integer, layer::Integer)
	ccall((:GMT_Get_Pixel, thelib), UInt64, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Cint, Cint, Cint), API, header, row, col, layer)
end

function GMT_Set_Index(API::Ptr{Cvoid}, header::Ptr{GMT_GRID_HEADER}, code)
	ccall((:GMT_Set_Index, thelib), Cint, (Ptr{Cvoid}, Ptr{GMT_GRID_HEADER}, Ptr{UInt8}), API, header, code)
end
=#

function GMT_Alloc_Segment(API::Ptr{Cvoid}, family::Integer, n_rows::Integer, n_columns::Integer, header, S::Ptr{Cvoid})
	if (family == GMT_IS_DATASET || family == GMT_WITH_STRINGS)
		ret_type = Ptr{GMT_DATASEGMENT}
	elseif (family == GMT_IS_TEXTSET)
		ret_type = Ptr{GMT_TEXTSEGMENT}
	else
		error("Bad family type")
	end

	ptr = ccall((:GMT_Alloc_Segment, thelib), Ptr{Cvoid}, (Cstring, UInt32, UInt64, UInt64, Ptr{UInt8}, Ptr{Cvoid}),
							 API, family, n_rows, n_columns, header, S)

	convert(ret_type, ptr)
end

#=
function GMT_Set_Columns(API::Ptr{Cvoid}, n_columns::Integer, mode)
	ccall((:GMT_Set_Columns, thelib), Cint, (Ptr{Cvoid}, UInt32, UInt32), API, n_columns, mode)
end

function GMT_Destroy_Group(API::Ptr{Cvoid}, obj::Ptr{Cvoid}, n_items::Integer)
	ccall((:GMT_Destroy_Group, thelib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cuint), API, obj, n_items)
end

function GMT_Change_Layout(API::Ptr{Cvoid}, family::Integer, code, mode::Integer, obj::Ptr{Cvoid}, data::Ptr{Cvoid}=C_NULL, alpha::Ptr{Cvoid}=C_NULL)
	ccall((:GMT_Change_Layout, thelib), Cint, (Ptr{Cvoid}, Cuint, Ptr{UInt8}, Cuint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
	      API, family, code, mode, obj, data, alpha)
end
=#

function GMT_Duplicate_String(API::Ptr{Cvoid}, str)
	ccall((:GMT_Duplicate_String, thelib), Ptr{UInt8}, (Cstring, Ptr{UInt8}), API, str)
end

function GMT_Open_VirtualFile(API::Ptr{Cvoid}, family::Integer, geometry::Integer, dir::Integer, data, name)
	ccall((:GMT_Open_VirtualFile, thelib), Cint, (Cstring, UInt32, UInt32, UInt32, Ptr{Cvoid}, Ptr{UInt8}), API, family, geometry, dir, data, name)
end
function GMT_Close_VirtualFile(API::Ptr{Cvoid}, str)
	ccall((:GMT_Close_VirtualFile, thelib), Cint, (Cstring, Ptr{UInt8}), API, str)
end

function GMT_Read_VirtualFile(API::Ptr{Cvoid}, str)
	ccall((:GMT_Read_VirtualFile, thelib), Ptr{Cvoid}, (Cstring, Ptr{UInt8}), API, str)
end

#=
function GMT_Init_VirtualFile(API::Ptr{Cvoid}, mode::Integer, name)
	ccall((:GMT_Init_VirtualFile, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, UInt32, Ptr{UInt8}), API, mode, name)
end

function GMT_Put_Vector(API::Ptr{Cvoid}, V::GMT_VECTOR, col::Integer, tipo::Integer, vector)
	ccall((:GMT_Put_Vector, thelib), Cint, (Ptr{Cvoid}, GMT_VECTOR, UInt32, UInt32, Ptr{Cvoid}), API, V, col, tipo, vector)
end
function GMT_Get_Vector(API::Ptr{Cvoid}, V::Ptr{GMT_VECTOR}, col::Integer)
	ccall((:GMT_Get_Vector, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{GMT_VECTOR}, UInt32), API, V, col)
end

function GMT_Put_Matrix(API::Ptr{Cvoid}, M::GMT_MATRIX, tipo::Integer, matrix)
	ccall((:GMT_Put_Matrix, thelib), Cint, (Ptr{Cvoid}, GMT_MATRIX, UInt32, Ptr{Cvoid}), API, M, tipo, matrix)
end
function GMT_Get_Matrix(API::Ptr{Cvoid}, M::Ptr{GMT_MATRIX})
	ccall((:GMT_Get_Matrix, thelib), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{GMT_MATRIX}), API, M)
end

function GMT_Error_Message(API::Ptr{Cvoid})
	ccall((:GMT_Error_Message, thelib), Ptr{UInt8}, (Ptr{Cvoid},), API)
end
=#

# ------------------ Development function in 5.4.0 ------------------------------------------------------------

# -------------------------------------------------------------------------------------------------------------
function GMT_Set_AllocMode(API::Ptr{Cvoid}, family::Integer, object)
	ccall((:GMT_Set_AllocMode, thelib), Cint, (Cstring, UInt32, Ptr{Cvoid}), API, family, object)
end

function gmt_manage_workflow(API::Ptr{Cvoid}, mode::Integer, texto)
	ccall((:gmt_manage_workflow, thelib), Cint, (Cstring, UInt32, Cstring), API, mode, texto)
end

function GMT_Get_Version()
	ver = ccall((:GMT_Get_Version, thelib), Cfloat, (Ptr{Cvoid}, Ptr{Cuint}, Ptr{Cuint}, Ptr{Cuint}), C_NULL, C_NULL, C_NULL, C_NULL)
end
function GMT_Get_Version(major, minor, patch)
	ver = ccall((:GMT_Get_Version, thelib), Cfloat, (Ptr{Cvoid}, Ptr{Cuint}, Ptr{Cuint}, Ptr{Cuint}), C_NULL, major, minor, patch)
end

function GMT_Get_Ctrl(API::Ptr{Cvoid})
	ccall((:GMT_Get_Ctrl, thelib), Ptr{Cvoid}, (Cstring,), API)
end

function gmt_getpen(API::Ptr{Cvoid}, buffer, P)
	GMT_ = GMT_Get_Ctrl(API)
	ccall((:gmt_getpen, thelib), Cint, (Cstring, Ptr{Cuint}, Ref{GMT_PEN}), GMT_, buffer, P)
end

function gmt_ogrread(API::Ptr{Cvoid}, fname::String)
	GMT_ = GMT_Get_Ctrl(API)
	ccall((:gmt_ogrread, thelib), Ptr{OGR_FEATURES}, (Cstring, Ptr{UInt8}), GMT_, fname)
end

#=
function gmt_put_history(API::Ptr{Cvoid})
	ccall((:gmt_put_history, thelib), Cint, (Cstring,), GMT_Get_Ctrl(API))
end
=#


text:
module Trello
using OrderedCollections
using HTTP
using JSON2

export TrelloCred
export get_boards, get_lists, get_labels, get_cards
export create_board, create_list, create_card
export delete_board

include("auth.jl")
include("api_components.jl")
include("api.jl")

end #module

ext: jl


text:
using Toys
using Documenter

makedocs(;
    modules=[Toys],
    authors="QuantEx team",
    repo="https://git.ichec.ie/quantex/Toys.jl/blob/{commit}{path}#L{line}",
    sitename="Toys.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://quantex.gitlab.io/Toys.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
module BaslerCamera

export Camera, get_data, connect, start_acquisition, stop_acquisition

export start_ffmpeg, end_ffmpeg, change_ffmpeg_folder

using Libdl, FFMPEG

const mydl = string(dirname(Base.source_path()),"/../deps/BaslerCamera.so")

function __init__()
    try
        init_pylon()
    catch
        println("Could not successfully load Pylon libraries. Are they installed?")
    end
end

mutable struct Camera
    cam::Ptr{Nothing}
    connected::Bool
    acquiring::Bool
    recording::Bool
    w::Int64
    h::Int64
    bytes::Int64
    ffmpeg_path::String
    ffmpeg_input_opts::String
    ffmpeg_output_opts::String
    output_folder::String
    num_cam::Int64
end

function Camera(h,w,_bytes,num_cam=1)
    cam = init_camera(num_cam)
    Camera(cam,false,false,false,w,h,_bytes,FFMPEG.FFMPEG_jll.ffmpeg_path,"","",".",num_cam)
end

function get_data(cam::Camera)
    get_camera_data(cam.cam,cam.w,cam.h,cam.num_cam)
end

function connect(cam::Camera)
    connect_camera(cam.cam)
    sleep(1.0)
    change_resolution(cam.cam,cam.w,cam.h)
end

function start_acquisition(cam::Camera)
    start_acquisition(cam.cam)
end

function stop_acquisition(cam::Camera)
    stop_acquisition(cam.cam)
end

function init_pylon()
    ccall((:initPylon,mydl),Nothing,())
    nothing
end

function init_camera(num_cam=1)
    cam = ccall((:newMyCamera,mydl),Ptr{Nothing},(Int32,),num_cam)
end

function connect_camera(cam::Ptr{Nothing})
    ccall((:MyCamera_Connect,mydl),Nothing,(Ptr{Nothing},),cam)
    nothing
end

function connect_camera(cam::Ptr{Nothing},serial::String)
    ccall((:MyCamera_ConnectBySerial,mydl),Nothing,(Ptr{Nothing},Ptr{UInt8}),cam,serial)
    nothing
end

function start_acquisition(cam::Ptr{Nothing})
    ccall((:MyCamera_StartAcquisition,mydl),Nothing,(Ptr{Nothing},),cam)
    nothing
end

function stop_acquisition(cam::Ptr{Nothing})
    ccall((:MyCamera_StopAcquisition,mydl),Nothing,(Ptr{Nothing},),cam)
    nothing
end

function get_camera_data(cam::Ptr{Nothing},w,h,n_cam=1)
    ccall((:MyCamera_GrabFrames,mydl),Nothing,(Ptr{Nothing},),cam)

    #If there is no TTL trigger, the data matrix
    #Associated with the camera may not be initialized
    #and we will just return a blank matrix
    grabbed = ccall((:MyCamera_GetFramesGrabbed,mydl),Bool,(Ptr{Nothing},),cam)
    if (grabbed)
        a = ccall((:MyCamera_GetData,mydl),Ptr{UInt8},(Ptr{Nothing},),cam)
        b=unsafe_wrap(Array,a,w*h*n_cam)
        c=reshape(b,(convert(Int64,h),convert(Int64,w*n_cam)))
    else
        c = zeros(UInt8,w,h*n_cam)
    end

    (c, grabbed)
end

function start_ffmpeg(cam::Ptr{Nothing})
    ccall((:MyCamera_StartFFMPEG,mydl),Nothing,(Ptr{Nothing},),cam)
    nothing
end

function end_ffmpeg(cam::Ptr{Nothing})
    ccall((:MyCamera_EndFFMPEG,mydl),Nothing,(Ptr{Nothing},),cam)
    nothing
end

function change_resolution(cam,w,h)
    ccall((:MyCamera_changeResolution,mydl),Nothing,(Ptr{Nothing},Int32,Int32),cam,w,h)
    nothing
end

function change_ffmpeg_folder(cam,folder)

    ccall((:MyCamera_ChangeFolder,mydl),Nothing,(Ptr{Nothing},Ptr{UInt8}),cam,pointer(string(folder,"/")))

    nothing
end

function change_camera_config(cam,path)
    ccall((:MyCamera_ChangeCameraConfig,mydl),Nothing,(Ptr{Nothing},Ptr{UInt8}),cam,pointer(string(path)))

    nothing
end

function change_ffmpeg_path(cam,path)
    ccall((:MyCamera_ChangeFFMPEG,mydl),Nothing,(Ptr{Nothing},Ptr{UInt8}),cam,pointer(string(path)))

    nothing
end

function change_ffmpeg_input_opts(cam,opts)
    ccall((:MyCamera_ChangeFFMPEGInputOptions,mydl),Nothing,(Ptr{Nothing},Ptr{UInt8}),cam,pointer(string(opts)))

    nothing
end

function change_ffmpeg_output_opts(cam,opts)
    ccall((:MyCamera_ChangeFFMPEGOutputOptions,mydl),Nothing,(Ptr{Nothing},Ptr{UInt8}),cam,pointer(string(opts)))

    nothing
end

function change_bytes(cam,_bytes)
    ccall((:MyCamera_ChangeBytes,mydl),Nothing,(Ptr{Nothing},Int32),cam,_bytes)

    nothing
end

end


text:
using Test

@static if VERSION ≥ v"1.6"
    @testset "prefs" begin
        include("prefs.jl")
    end
end

include("basic.jl")
include("configs.jl")
include("dataframe.jl")


text:
struct ConvergenceException <: Exception
end


text:
using Test
using YaoArrayRegister, Random, LinearAlgebra, SparseArrays, BitBasis

@testset "broadcast register" begin
    reg = rand_state(5; nbatch = 3)
    @test typeof.(reg)[1] <: ArrayReg{2,ComplexF64,<:SubArray}
    @test [reg...] |> length == 3
    @test [rand_state(3; nbatch=2)...] |> length == 2
end

@testset "arithmetics" begin
    nbit = 5
    reg1 = zero_state(5)
    reg2 = arrayreg(bit"00100")
    @test reg1 != reg2
    @test statevec(reg2) == onehot(ComplexF64, nbit, 4)
    reg3 = reg1 + reg2
    reg4 = (reg1 + reg2)'

    @test statevec(reg3) == onehot(ComplexF64, nbit, 4) + onehot(ComplexF64, nbit, 0)
    @test statevec(reg3 |> normalize!) ==
          (onehot(ComplexF64, nbit, 4) + onehot(ComplexF64, nbit, 0)) / sqrt(2)
    @test statevec(reg4 |> normalize!) ==
          (onehot(ComplexF64, nbit, 4) + onehot(ComplexF64, nbit, 0))' / sqrt(2)
    @test (reg1 + reg2 - reg1) == reg2
    @test reg1' + reg2' - reg1' == reg2'
    @test isnormalized(reg4)
    @test isnormalized(reg3)

    @test statevec(-reg4) == -statevec(reg4)
    @test statevec(-reg3) == -statevec(reg3)
    reg = rand_state(4)
    @test all(state(reg + (-reg)) .== 0)
    @test all(state(reg * 2 - reg / 0.5) .== 0)

    reg = rand_state(3)
    @test norm(reg) ≈ 1
    @test reg' * reg ≈ 1

    @test state(reg1 * 2) == state(reg1) * 2
    @test state(reg1' * 2) == state(reg1') * 2
    @test reg1 * 2 == 2 * reg1
    @test reg1' * 2 == 2 * reg1'
    reg = rand_state(3; nbatch = 2) * 2
    @test norm(reg) ≈ [2.0, 2.0]
end

@testset "partial ⟨bra|ket⟩" begin
    bra = arrayreg(bit"10")
    ket = arrayreg(bit"100") + 2 * arrayreg(bit"110") + 3 * arrayreg(bit"111")

    focus!(ket, 2:3)
    t = bra' * ket
    relax!(t, 1)
    @test state(t) ≈ [1, 0]

    relax!(ket, 2:3)
    focus!(ket, 1)
    focus!(bra, 2)
    @test_throws ErrorException bra' * ket

    reg1 = rand_state(5; nbatch = 10)
    reg2 = rand_state(5; nbatch = 10)
    @test reg1' * reg2 ≈ reg1' .* reg2
    reg1 = rand_state(2; nbatch = 10)
    reg2 = rand_state(5; nbatch = 10)
    focus!(reg2, 2:3)
    @test all(reg1' * reg2 .≈ reg1' .* reg2)
end

@testset "inplace funcs" begin
    for reg in [
            rand_state(5; nbatch = NoBatch()),
            rand_state(5; nbatch = 3),
            transpose_storage(rand_state(5; nbatch = 3)),
        ]
        reg0 = copy(reg)
        @test regscale!(reg, 0.3) ≈ 0.3 * reg0
        reg1 = rand_state(5; nbatch = nbatch(reg))
        reg2 = rand_state(5; nbatch = nbatch(reg))
        reg10 = copy(reg1)
        regsub!(reg1, reg2)
        @test reg1 ≈ reg10 - reg2

        reg1 = rand_state(5; nbatch = nbatch(reg))
        reg2 = rand_state(5; nbatch = nbatch(reg))
        reg10 = copy(reg1)
        regadd!(reg1, reg2)
        @test reg1 ≈ reg10 + reg2
    end
end

@testset "more (push test coverage)" begin
    reg1 = focus!(rand_state(5; nbatch=5), (2, 3))
    reg2 = focus!(rand_state(5), (2, 3))
    @test fidelity(reg1, reg2) ≈ fidelity(reg2, reg1)
    @test join(reg1) == reg1
    us = uniform_state(3; nlevel=3, nbatch=2)
    @test nbatch(us) == 2
    @test nlevel(us) == 3
    println(reg1)
    println(reg2)
    @test von_neumann_entropy(repeat(reg2, 3), [2,1]) ≈ fill(von_neumann_entropy(reg2, [1,2]), 3)
end


text:
using PyPlot, Polynomials, CoolProp, Roots
using JFVM
include("../functions/rel_perms_real.jl")

# some data
MW_water = 0.018 # [kg/mol]
MW_oil   = 0.142 # [kg/mol]
MW_DME   = 0.04607 # [kg/mol]

# Fig 4, SPE-177919, T = 70 degC, p = 2000 psi
x_DME_oil = [0.0, 0.04, 0.10, 0.39, 0.62]         # [mol frac]
w_DME_oil = MW_DME*x_DME_oil./(MW_DME*x_DME_oil+MW_oil*(1-x_DME_oil))
ρ_oil     = [0.748, 0.746, 0.744, 0.722, 0.697]*1000   # [g/cm^3]*1000 = [kg/m^3]

rho_fit = polyfit(w_DME_oil, ρ_oil, 2)

# plot(w_DME_oil, ρ_oil, "o", w_DME_oil, polyval(rho_fit, w_DME_oil))
# xlabel("DME mass fraction in oil [-]")
# ylabel("Oil density [g/cm^3]")

# experimental data from SPE-179771
# data for 1.4% brine and 71 degC and 2000 psia
x_DME_oil = [0.15, 0.30, 0.5, 0.6]   # [mol frac] [0.5, 13] is a madeup point
w_DME_oil = MW_DME*x_DME_oil./(MW_DME*x_DME_oil+MW_oil*(1-x_DME_oil)) # [mass frac]
K_value   = [17.0, 15.0, 13.0, 12.5]  # [-]
x_DME_water = x_DME_oil./K_value # [mole frac]
w_DME_water = MW_DME*x_DME_water./(MW_DME*x_DME_water+MW_water*(1-x_DME_oil))
K_value_mass = w_DME_oil./w_DME_water
K_fit = polyfit(w_DME_water, K_value_mass, 2)
w_DME_plot = collect(linspace(minimum(w_DME_water), maximum(w_DME_water), 50))
# plot(w_DME_water, K_value_mass, "o", w_DME_plot, polyval(K_fit,w_DME_plot))
# xlabel("DME mass fraction in water [-]")
# ylabel("K-value [-]")

# Fig 4, SPE-177919, T = 70 degC, p = 2000 psi
x_DME_oil = [0.0, 0.04, 0.10, 0.39, 0.62] # [mol frac]
w_DME_oil = MW_DME*x_DME_oil./(MW_DME*x_DME_oil+MW_oil*(1-x_DME_oil))
μ_oil = [1.04, 0.96, 0.84, 0.52, 0.34]/1000.0    # [cP]/1000 = [Pa.s]
mu_oil_fit = polyfit(w_DME_oil, μ_oil, 2)
w_DME_plot = collect(linspace(minimum(w_DME_oil), maximum(w_DME_oil), 50))
# plot(w_DME_oil, μ_oil, "o", w_DME_plot, polyval(mu_oil_fit, w_DME_plot))
# xlabel("DME mass fraction in oil [-]")
# ylabel("Oil viscosity [Pa.s]")

T0 = 70 + 273.15 # [K]
p0 = 2000/14.7*1e5   # [Pa]
mu_water = PropsSI("viscosity", "T", T0, "P", p0, "water") # [Pa.s]
mu_water_DME = [mu_water, 2*mu_water]
w_DME_water = [0.0, maximum(w_DME_oil)/maximum(K_value_mass)]
mu_water_fit = polyfit(w_DME_water, mu_water_DME, 1)
# plot(w_DME_water, mu_water_DME, "o", w_DME_water, mu_water_fit(w_DME_water))
# xlabel("DME mass fraction in water [-]")
# ylabel("Water-DME viscosity [Pa.s]")

# T = 70 + 273.15 # [K]
# p = 2000/14.7*1e5   # [Pa]
rho_water = PropsSI("D", "T", T0, "P", p0, "water")
# Corey rel-perm parameters
krw0 = 0.2
kro0 = 0.8
n_o  = 2.0
n_w  = 2.0
swc  = 0.08
sor  = 0.3
sorting_factor = 2.4
pce = 100 # [Pa]
pc0 = 1e5 # [Pa]
contact_angle = deg2rad(20) # [radian]

perm_val  = 0.01e-12 # [m^2] permeability
poros_val = 0.40     # [-] porosity

KRW  = sw -> krw.(sw, krw0, sor, swc, n_w)
dKRW = sw -> dkrwdsw.(sw, krw0, sor, swc, n_w)
KRO  = sw -> kro.(sw, kro0, sor, swc, n_o)
dKRO = sw -> dkrodsw.(sw, kro0, sor, swc, n_o)
# PC   = sw -> pc_imb(sw, pce, swc, sor, teta=contact_angle,
#                 labda=sorting_factor, b = 0.0, pc01=pc0, pc02=pc0)
# dPC  = sw -> dpc_imb(sw, pce, swc, sor, teta=contact_angle,
#                 labda1=sorting_factor, labda2=sorting_factor,
#                 b = 1.0, pc01=pc0, pc02=pc0)
# PC2  = sw -> pc_imb2(sw, pce, swc, sor, teta=contact_angle,
#                 labda=sorting_factor, b = 0.0, pc_star1=pc0,
#                 pc_star2=pc0)
# dPC2 = sw -> dpc_imb2(sw, pce, swc, sor, teta=contact_angle,
#                 labda1=sorting_factor, labda2=sorting_factor,
#                 b = 1.0, pc_star1=pc0, pc_star2=pc0)
PC  = sw -> pc_imb3.(sw, pce, swc, sor, teta=contact_angle,
                labda=sorting_factor, b = 0.0, pc_star1=pc0,
                pc_star2=pc0)
dPC = sw -> dpc_imb3.(sw, pce, swc, sor, teta=contact_angle,
                labda1=sorting_factor, labda2=sorting_factor,
                b = 1.0, pc_star1=pc0, pc_star2=pc0)
d2PC = sw -> d2pc_imb3.(sw, pce, swc, sor, teta=contact_angle,
                labda1=sorting_factor, labda2=sorting_factor,
                b = 1.0, pc_star1=pc0, pc_star2=pc0)
# PCdrain2 = sw -> pc_drain2(sw, pce, swc, labda=sorting_factor, pc_star=pc0)
# PCdrain3 = sw -> pc_drain3(sw, pce, swc, labda=sorting_factor, pc_star=pc0)
ρ_oil = w_DME_oil -> polyval(rho_fit, w_DME_oil)
dρ_oil = w_DME_oil -> polyval(polyder(rho_fit), w_DME_oil)
μ_water = w_DME_water -> polyval(mu_water_fit, w_DME_water)
dμ_water = w_DME_water -> polyval(polyder(mu_water_fit), w_DME_water)
μ_oil = w_DME_oil -> polyval(mu_oil_fit, w_DME_oil)
dμ_oil = w_DME_oil -> polyval(polyder(mu_oil_fit), w_DME_oil)
Kval = w_DME_water -> polyval(K_fit, w_DME_water)
dKval = w_DME_water -> polyval(polyder(K_fit), w_DME_water)
w_oil = w_DME_water -> polyval(Poly([0.0,1.0])*K_fit, w_DME_water)
dw_oil = w_DME_water -> polyval(polyder(Poly([0.0,1.0])*K_fit), w_DME_water)

sw_plot = collect(linspace(0,1,1000))
# plot(sw_plot, KRW.(sw_plot), sw_plot, KRO.(sw_plot))
# axis([0,1,0,1])

# plot(sw_plot, PC.(sw_plot))
# plot([swc, swc], [-pc0, pc0], [1-sor, 1-sor], [-pc0, pc0],
#     [0,1], [0,0], "--")
# xlabel("Water saturation (Sw) [-]")
# ylabel("Capillary pressure (Pc) [Pa]")
sw_imb_end = fzero(PC, [swc, 1-sor])
# plot(sw_imb_end, 0, "o")
# axis([0, 1, -pc0, pc0])
# yscale("symlog")

Lx   = 1.0 # [m]
Nx  = 200  # number of grids in the x direction
m   = createMesh1D(Nx, Lx)

u_inj = 1.0/(3600*24) # 1 m/day to m/s injection velocity
c_inj = 0.2 # DME mass fraction in the injected water

BCp = createBC(m) # pressure boundary condition
BCs = createBC(m) # saturation boundary condition
BCc = createBC(m) # concentration (DME in water) boundary condition

BCp.right.a[:] = 0.0
BCp.right.b[:] = 1.0
BCp.right.c[:] = p0

BCp.left.a[:]  = perm_val/μ_water(c_inj)
BCp.left.b[:]  = 0.0
BCp.left.c[:]  = -u_inj

BCs.left.a[:]  = 0.0
BCs.left.b[:]  = 1.0
BCs.left.c[:]  = 1.0

# BCs.right.a[:]  = 0.0
# BCs.right.b[:]  = 1.0
# BCs.right.c[:]  = sw_imb_end

BCc.left.a[:]  = 0.0
BCc.left.b[:]  = 1.0
BCc.left.c[:]  = c_inj

# discretize
M_bc_p, RHS_bc_p = boundaryConditionTerm(BCp)
M_bc_s, RHS_bc_s = boundaryConditionTerm(BCs)
M_bc_c, RHS_bc_c = boundaryConditionTerm(BCc)

c0  = 0.0 # [mass frac] initial DME concentration in the water phase
sw0 = swc+0.05 # [vol frac] initial water saturation

p_init  = createCellVariable(m, p0, BCp)
c_init  = createCellVariable(m, c0, BCc)
sw_init = createCellVariable(m, sw0, BCs)

# new values of each variable
p_val  = createCellVariable(m, p0, BCp)
c_val  = createCellVariable(m, c0, BCc)
sw_val = createCellVariable(m, sw0, BCs)

Δsw_init = sw_init - sw_val # we solve for this variable
Δc_init  = c_init - c_val   # we solve for this variable

# Other cell variables
k = createCellVariable(m, perm_val)
ϕ = createCellVariable(m, poros_val)

n_pv    = 0.4 # number of injected pore volumes
t_final = n_pv*Lx/(u_inj/poros_val) # [s] final time
dt0     = t_final/n_pv/Nx/10.0 # [s] time step
dt      = dt0

# outside the loop: initialization
uw       = gradientTerm(p_val) # only for initialization of the water velocity vector
k_face   = harmonicMean(k)     # permeability on the cell faces

# this whole thing goes inside two loops (Newton and time)
tol_s = 1e-7
tol_c = 1e-7
max_change_c = 0.1 # 10 % relative change
max_change_s = 0.1 # 10 % relative change
max_int_loop = 5
t = dt
while t<t_final
  error_s = 2*tol_s
  error_c = 2*tol_c
  loop_countor = 0
  while error_s>tol_s || error_c>tol_c
    loop_countor += 1
    if loop_countor > max_int_loop
      sw_val = copyCell(sw_init)
      c_val  = copyCell(c_init)
      p_val  = copyCell(p_init)
      dt = dt/3.0
      break
    end
    ∇p0      = gradientTerm(p_val)
    ∇s0      = gradientTerm(sw_val)

    c_face   = upwindMean(c_val, uw)
    c_oil         = cellEval(w_oil, c_val) # [mass frac] DME concentration in oil
    c_oil_face    = faceEval(w_oil, c_face)
    sw_face       = upwindMean(sw_val, uw)

    mu_water_face = faceEval(μ_water, c_face)
    dμ_water_face = faceEval(dμ_water, c_face)
    mu_oil_face   = faceEval(μ_oil, c_oil_face)
    dμ_oil_face   = faceEval(dw_oil, c_face).*faceEval(dμ_oil,c_oil_face)
    rho_oil_cell  = cellEval(ρ_oil, c_oil)
    rho_oil_face  = faceEval(ρ_oil, c_oil_face)
    dρ_oil_cell   = cellEval(dw_oil, c_val).*cellEval(dρ_oil, c_oil)
    dρ_oil_face   = faceEval(dw_oil, c_face).*faceEval(dρ_oil, c_oil_face)
    dρc_oil_cell = dρ_oil_cell.*c_oil+cellEval(dw_oil, c_val).*rho_oil_cell
    dρc_oil_face = dρ_oil_face.*c_oil_face+faceEval(dw_oil, c_face).*rho_oil_face
    dpc_face      = faceEval(dPC, sw_face)
    d2pc_face      = faceEval(d2PC, sw_face)

    dρ_dμ_oil = (dρ_oil_face.*mu_oil_face-dμ_oil_face.*rho_oil_face)./(mu_oil_face.*mu_oil_face)
    dρc_dμ_oil = (dρc_oil_face.*mu_oil_face-dμ_oil_face.*rho_oil_face.*c_oil_face)./(mu_oil_face.*mu_oil_face)

    krw_face      = faceEval(KRW,sw_face)
    kro_face      = faceEval(KRO,sw_face)

    λ_w_face      = k_face.*krw_face./mu_water_face
    λ_o_face      = k_face.*kro_face./mu_oil_face
    uw            = -λ_w_face.*∇p0
    uo            = -λ_o_face.*∇p0
    ut            = uw + uo

    dλ_w_face     = k_face.*faceEval(dKRW,sw_face)./mu_water_face
    dλ_o_face     = k_face.*faceEval(dKRO,sw_face)./mu_oil_face

    # λ_w_face     = harmonicMean(λ_w) # Wrong!
    # λ_o_face = harmonicMean(λ_o)     # Wrong!

    # Δsw_init = sw_init - sw_val # we solve for this variable
    # Δc_init  = c_init - c_val   # we solve for this variable

    # water mass balance (wmb)
    M_t_s_wmb, RHS_t_s_wmb = transientTerm(sw_init, dt, rho_water*ϕ)
    M_d_p_wmb = diffusionTerm(-rho_water*λ_w_face)
    M_a_s_wmb = convectionUpwindTerm(-rho_water*dλ_w_face.*∇p0, ut)
    M_a_c_wmb = convectionUpwindTerm(rho_water*dμ_water_face./mu_water_face.*λ_w_face.*∇p0, ut)

    # oil mass balance (omb)
    M_t_s_omb, RHS_t_s_omb = transientTerm(sw_init, dt, -rho_oil_cell.*ϕ)
    M_t_c_omb, RHS_t_c_omb = transientTerm(c_init, dt, (1.0-sw_val).*dρ_oil_cell.*ϕ)
    M_d_p_omb = diffusionTerm(-rho_oil_face.*λ_o_face)
    # M_a_c_omb = convectionUpwindTerm(-dρ_dμ_oil.*k_face.*kro_face.*(∇p0+dpc_face.*∇s0))
    M_a_c_omb = convectionUpwindTerm(-dρ_dμ_oil.*k_face.*kro_face.*∇p0, ut)
    M_a_s_omb = convectionUpwindTerm(-rho_oil_face.*(dλ_o_face.*∇p0), ut)
                                    # +(dλ_o_face.*dpc_face+λ_o_face.*d2pc_face).*∇s0))
    M_d_s_omb = 0.0 #diffusionTerm(-rho_oil_face.*λ_o_face.*dpc_face)

    # DME mass balance (dme)
    M_t_s_dme, RHS_t_s_dme = transientTerm(sw_init, dt, (rho_water*c_val-rho_oil_cell.*c_oil).*ϕ)
    M_t_c_dme, RHS_t_c_dme = transientTerm(c_init, dt, (rho_water*sw_val+
                                           (1.0-sw_val).*dρc_oil_cell).*ϕ)
    M_d_p_dme  = diffusionTerm(-rho_water*c_face.*λ_w_face-rho_oil_face.*c_oil_face.*λ_o_face)
    M_a_s_dme  = convectionUpwindTerm(-(rho_water*c_face.*dλ_w_face+
                                      rho_oil_face.*c_oil_face.*dλ_o_face).*∇p0, ut)
                                    #   -rho_oil_face.*c_oil_face.*
                                    #   (dλ_o_face.*dpc_face+λ_o_face.*d2pc_face).*∇s0)
    M_a_c_dme  = convectionUpwindTerm(-(rho_water*(mu_water_face-dμ_water_face.*c_face)
                                      ./mu_water_face.*λ_w_face+
                                        dρc_dμ_oil.*k_face.*kro_face).*∇p0, ut)
                                    #   -dρc_dμ_oil.*k_face.*kro_face.*dpc_face.*∇s0)
    M_d_s_dme  = 0.0 # diffusionTerm(-c_oil_face.*rho_oil_face.*λ_o_face.*dpc_face)

    # create the PDE system M [p;s;c]=RHS
    # x_val = [p_val.value[:]; sw_val.value[:]; c_val.value[:]]
    M = [M_bc_p+M_d_p_wmb   M_t_s_wmb+M_a_s_wmb    M_a_c_wmb;
         M_d_p_omb   M_bc_s+M_t_s_omb+M_a_s_omb+M_d_s_omb    M_t_c_omb+M_a_c_omb;
         M_d_p_dme    M_t_s_dme+M_a_s_dme+M_d_s_dme    M_bc_c+M_t_c_dme+M_a_c_dme;]

    RHS = [RHS_bc_p+RHS_t_s_wmb+M_a_s_wmb*sw_val.value[:]+M_a_c_wmb*c_val.value[:];
           RHS_bc_s+RHS_t_s_omb+RHS_t_c_omb+(M_a_s_omb+M_d_s_omb)*sw_val.value[:]+M_a_c_omb*c_val.value[:];
           RHS_bc_c+RHS_t_s_dme+RHS_t_c_dme+(M_a_s_dme+M_d_s_dme)*sw_val.value[:]+M_a_c_dme*c_val.value[:];]

    # x_sol = solveLinearPDE(m, M, RHS)
    x_sol = M\RHS

    p_new = x_sol[1:Nx+2]
    s_new = x_sol[Nx+3:2*Nx+4]
    c_new = x_sol[2*Nx+5:end]
    error_s = sumabs(s_new[2:end-1]-sw_val.value[2:end-1])
    error_c = sumabs(c_new[2:end-1]-c_val.value[2:end-1])
    # println(error_s)
    # println(error_c)
    p_val.value[:] = p_new[:]
    sw_val.value[:] = s_new[:]
    c_val.value[:] = c_new[:]
  end
  if loop_countor<max_int_loop
    p_init = copyCell(p_val)
    sw_init = copyCell(sw_val)
    c_init = copyCell(c_val)
    t +=dt
    dt = dt0
    # println("time is $t [s]")
    println("progress is $(t/t_final*100) [%]")
  end
end
visualizeCells(sw_init)
visualizeCells(c_init)
savefig("profile.png")

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
function open_vault!(account, output; rep)
    account.vault_pct -= 100.0
    account.wc_count[2] += 3
    account.wc_count[3] += 2
    account.wc_count[4] += 1
    output.wcs_from_vault[rep,2] += 3
    output.wcs_from_vault[rep,3] += 2
    output.wcs_from_vault[rep,4] += 1
    account
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
# get the expression of the inverse function
function dual_func(m::Module, fname, args, ts, body)
    :(function $(:(~$fname))($(args...)) where {$(ts...)};
            $(dual_body(m, body)...);
        end)
end

# get the function name of the inverse function
function dual_fname(op)
    @smatch op begin
        :($x::$tp) => :($x::$invtype($tp))
        :(~$x) => x
        _ => :($(gensym("~$op"))::$_typeof(~$op))
    end
end
_typeof(x) = typeof(x)
_typeof(x::Type{T}) where T = Type{T}

"""
    dual_ex(m::Module, ex)

Get the dual expression of `ex`.
"""
function dual_ex(m::Module, ex)
    @smatch ex begin
        :(($t1=>$t2)($x)) => :(($t2=>$t1)($x))
        :(($t1=>$t2).($x)) => :(($t2=>$t1).($x))
        :($x ↔ $y) => dual_swap(x, y)
        :($s[end+1] ← $x) => :($s[end] → $x)
        :($s[end] → $x) => :($s[end+1] ← $x)
        :($x → $val) => :($x ← $val)
        :($x ← $val) => :($x → $val)
        :($f($(args...))) => startwithdot(f) ? :($(getdual(removedot(sym))).($(args...))) : :($(getdual(f))($(args...)))
        :($f.($(args...))) => :($(getdual(f)).($(args...)))
        :($a += $b) => :($a -= $b)
        :($a .+= $b) => :($a .-= $b)
        :($a -= $b) => :($a += $b)
        :($a .-= $b) => :($a .+= $b)
        :($a *= $b) => :($a /= $b)
        :($a .*= $b) => :($a ./= $b)
        :($a /= $b) => :($a *= $b)
        :($a ./= $b) => :($a .*= $b)
        :($a ⊻= $b) => :($a ⊻= $b)
        :($a .⊻= $b) => :($a .⊻= $b)
        Expr(:if, _...) => dual_if(m, copy(ex))
        :(while ($pre, $post); $(body...); end) => begin
            Expr(:while, :(($post, $pre)), Expr(:block, dual_body(m, body)...))
        end
        :(for $i=$start:$step:$stop; $(body...); end) => begin
            Expr(:for, :($i=$stop:(-$step):$start), Expr(:block, dual_body(m, body)...))
        end
        :(for $i=$start:$stop; $(body...); end) => begin
            j = gensym("j")
            Expr(:for, :($j=$start:$stop), Expr(:block, :($i ← $stop-$j+$start), dual_body(m, body)..., :($i → $stop-$j+$start)))
        end
        :(for $i=$itr; $(body...); end) => begin
            Expr(:for, :($i=Base.Iterators.reverse($itr)), Expr(:block, dual_body(m, body)...))
        end
        :(@safe $line $subex) => Expr(:macrocall, Symbol("@safe"), line, subex)
        :(@cuda $line $(args...)) => Expr(:macrocall, Symbol("@cuda"), line, args[1:end-1]..., dual_ex(m, args[end]))
        :(@launchkernel $line $(args...)) => Expr(:macrocall, Symbol("@launchkernel"), line, args[1:end-1]..., dual_ex(m, args[end]))
        :(@inbounds $line $subex) => Expr(:macrocall, Symbol("@inbounds"), line, dual_ex(m, subex))
        :(@simd $line $subex) => Expr(:macrocall, Symbol("@simd"), line, dual_ex(m, subex))
        :(@threads $line $subex) => Expr(:macrocall, Symbol("@threads"), line, dual_ex(m, subex))
        :(@avx $line $subex) => Expr(:macrocall, Symbol("@avx"), line, dual_ex(m, subex))
        :(@invcheckoff $line $subex) => Expr(:macrocall, Symbol("@invcheckoff"), line, dual_ex(m, subex))
        :(begin $(body...) end) => Expr(:block, dual_body(m, body)...)
        :(nothing) => ex
        ::LineNumberNode => ex
        ::Nothing => ex
        :() => ex
        _ => error("can not invert target expression $ex")
    end
end

function dual_if(m::Module, ex)
    _dual_cond(cond) = @smatch cond begin
        :(($pre, $post)) => :(($post, $pre))
    end
    if ex.head == :if
        ex.args[1] = _dual_cond(ex.args[1])
    elseif ex.head == :elseif
        ex.args[1].args[2] = _dual_cond(ex.args[1].args[2])
    end
    ex.args[2] = Expr(:block, dual_body(m, ex.args[2].args)...)
    if length(ex.args) == 3
        if ex.args[3].head == :elseif
            ex.args[3] = dual_if(m, ex.args[3])
        elseif ex.args[3].head == :block
            ex.args[3] = Expr(:block, dual_body(m, ex.args[3].args)...)
        end
    end
    ex
end

function dual_swap(x, y)
    e1 = isemptyvar(x)
    e2 = isemptyvar(y)
    if e1 && !e2 || !e1 && e2
        :($(_dual_swap_var(x)) ↔ $(_dual_swap_var(y)))
    else
        :($y ↔ $x)
    end
end

_dual_swap_var(x) = @smatch x begin
    :($s[end+1]) => :($s[end])
    :($x::∅) => :($x)
    :($s[end]) => :($s[end+1])
    _ => :($x::∅)
end

export @code_reverse

"""
    @code_reverse ex

Get the reversed expression of `ex`.

```jldoctest; setup=:(using NiLangCore)
julia> @code_reverse x += exp(3.0)
:(x -= exp(3.0))
```
"""
macro code_reverse(ex)
    QuoteNode(dual_ex(__module__, ex))
end

getdual(f) = @smatch f begin
    :(~$f) => f
    _ => :(~$f)
end

function dual_body(m::Module, body)
    out = []
    # fix function LineNumberNode
    if length(body) > 1 && body[1] isa LineNumberNode && body[2] isa LineNumberNode
        push!(out, body[1])
        start = 2
    else
        start = 1
    end
    ptr = length(body)
    # reverse the statements
    len = 0
    while ptr >= start
        if  ptr-len==0 || body[ptr-len] isa LineNumberNode
            ptr-len != 0 && push!(out, body[ptr-len])
            for j=ptr:-1:ptr-len+1
                push!(out, dual_ex(m, body[j]))
            end
            ptr -= len+1
            len = 0
        else
            len += 1
        end
    end
    return out
end


text:
module Sessions

import SHA, HTTP, Dates, Logging
import Genie

const HTTP = HTTP

"""
    mutable struct Session

Represents a session object
"""
mutable struct Session
  id::String
  data::Dict{Symbol,Any}
end

Session(id::String) = Session(id, Dict{Symbol,Any}())

export Session

struct InvalidSessionIdException <: Exception
  msg::String
end
InvalidSessionIdException() =
  InvalidSessionIdException("Can't compute session id - make sure that secret_token!(token) is called in config/secrets.jl")

"""
    id() :: String

Generates a new session id.
"""
function id() :: String
  if isempty(Genie.secret_token())
    if !Genie.Configuration.isprod()
      @error "Empty Genie.secret_token(); using a temporary token"
      Genie.secret_token!()
    else
      throw(InvalidSessionIdException())
    end
  end

  try
    return join([
      Genie.secret_token(),
      bytes2hex(SHA.sha1(string(Dates.now()))),
      string(rand()),
      string(hash(Genie))
    ], ":") |> SHA.sha256 |> bytes2hex
  catch ex
    @error ex
    throw(InvalidSessionIdException())
  end
end

"""
    id(payload::Union{HTTP.Request,HTTP.Response}) :: String

Attempts to retrieve the session id from the provided `payload` object.
If that is not available, a new session id is created.
"""
function id(payload::Union{HTTP.Request,HTTP.Response}) :: String
  (Genie.Cookies.get(payload, Genie.config.session_key_name) !== nothing) &&
    ! isempty(Genie.Cookies.get(payload, Genie.config.session_key_name)) &&
      return Genie.Cookies.get(payload, Genie.config.session_key_name)

  id()
end

"""
    id(req::HTTP.Request, res::HTTP.Response) :: String

Attempts to retrieve the session id from the provided request and response objects.
If that is not available, a new session id is created.
"""
function id(req::HTTP.Request, res::HTTP.Response) :: String
  for r in [req, res]
    val = Genie.Cookies.get(r, Genie.config.session_key_name)
    (val !== nothing) && ! isempty(val) &&
      return val
  end

  id()
end

"""
    init() :: Nothing

Sets up the session functionality, if configured.
"""
function init() :: Nothing
  @eval Genie.config.session_storage === nothing && (Genie.config.session_storage = :File)
  @eval Genie.config.session_storage == :File && include(joinpath(@__DIR__, "session_adapters", "FileSession.jl"))

  push!(Genie.Router.pre_match_hooks, Genie.Sessions.start)
  push!(Genie.Router.pre_response_hooks, Genie.Sessions.persist)

  nothing
end

"""
    start(session_id::String, req::HTTP.Request, res::HTTP.Response; options = Dict{String,String}()) :: Tuple{Session,HTTP.Response}

Initiates a new HTTP session with the provided `session_id`.

# Arguments
- `session_id::String`: the id of the session object
- `req::HTTP.Request`: the request object
- `res::HTTP.Response`: the response object
- `options::Dict{String,String}`: extra options for setting the session cookie, such as `Path` and `HttpOnly`
"""
function start(session_id::String, req::HTTP.Request, res::HTTP.Response; options::Dict{String,Any} = Genie.config.session_options) :: Tuple{Session,HTTP.Response}
  Genie.Cookies.set!(res, Genie.config.session_key_name, session_id, options)

  load(session_id), res
end

"""
    start(req::HTTP.Request, res::HTTP.Response; options::Dict{String,String} = Dict{String,String}()) :: Session

Initiates a new default session object, generating a new session id.

# Arguments
- `req::HTTP.Request`: the request object
- `res::HTTP.Response`: the response object
- `options::Dict{String,String}`: extra options for setting the session cookie, such as `Path` and `HttpOnly`
"""
function start(req::HTTP.Request, res::HTTP.Response, params::Dict{Symbol,Any} = Dict{Symbol,Any}(); options::Dict{String,Any} = Genie.config.session_options) :: Tuple{HTTP.Request,HTTP.Response,Dict{Symbol,Any}}
  session, res = start(id(req, res), req, res; options = options)

  params[Genie.PARAMS_SESSION_KEY]   = session
  params[Genie.PARAMS_FLASH_KEY]     = begin
                                          if session !== nothing
                                            s = get(session, Genie.PARAMS_FLASH_KEY)
                                            if s === nothing
                                              ""
                                            else
                                              unset!(session, Genie.PARAMS_FLASH_KEY)
                                              s
                                            end
                                          else
                                            ""
                                          end
                                        end

  req, res, params
end

"""
    set!(s::Session, key::Symbol, value::Any) :: Session

Stores `value` as `key` on the `Session` object `s`.
"""
function set!(s::Session, key::Symbol, value::Any) :: Session
  s.data[key] = value

  s
end

"""
    get(s::Session, key::Symbol) :: Union{Nothing,Any}

Returns the value stored on the `Session` object `s` as `key`, wrapped in a `Union{Nothing,Any}`.
"""
function get(s::Session, key::Symbol) :: Union{Nothing,Any}
  haskey(s.data, key) ? (s.data[key]) : nothing
end

"""
    get(s::Session, key::Symbol, default::T) :: T where T

Attempts to retrive the value stored on the `Session` object `s` as `key`.
If the value is not set, it returns the `default`.
"""
function get(s::Session, key::Symbol, default::T) :: T where T
  val = get(s, key)

  val === nothing ? default : val
end

"""
    unset!(s::Session, key::Symbol) :: Session

Removes the value stored on the `Session` `s` as `key`.
"""
function unset!(s::Session, key::Symbol) :: Session
  delete!(s.data, key)

  s
end

"""
    isset(s::Session, key::Symbol) :: Bool

Checks wheter or not `key` exists on the `Session` `s`.
"""
function isset(s::Union{Session,Nothing}, key::Symbol) :: Bool
  s !== nothing && haskey(s.data, key)
end

"""
    persist(s::Session) :: Session

Generic method for persisting session data - delegates to the underlying `SessionAdapter`.
"""
function persist end

"""
    load(session_id::String) :: Session

Loads session data from persistent storage - delegates to the underlying `SessionAdapter`.
"""
function load end

"""
    session(params::Dict{Symbol,Any}) :: Sessions.Session

Returns the `Session` object associated with the current HTTP request.
"""
function session(params::Dict{Symbol,Any}) :: Sessions.Session
  ( (! haskey(params, Genie.PARAMS_SESSION_KEY) || params[Genie.PARAMS_SESSION_KEY] === nothing) ) &&
      (params[Genie.PARAMS_SESSION_KEY] = Sessions.start(params[Genie.PARAMS_REQUEST_KEY], params[Genie.PARAMS_RESPONSE_KEY])[1])

  params[Genie.PARAMS_SESSION_KEY]
end

end


text:
module JuliaCipher

using HDF5, JLD

""" saveChart(chart::Dict, chartLoc::String)

Handles saving the chart to file, given the location. """

function saveChart(chart::Dict, chartLoc::String)
	#saves a chart to file in a nice easy format. -NEEDS PERMISSIONS
	#handle bad permissions? where should it go?
	save(chartLoc, "chart", chart)

end

""" loadChart(chartLoc::String)

Handles loading a chart from given location, returns as a Dict. """

function loadChart(chartLoc::String)
	#grabs a chart from a file, and loads it properly
	#do a check if you even can first, then pop back a new chart? flag to allow new charting? just outright reject?
	return load(chartLoc)["chart"]
end

""" createChart(chartLoc::String; codeSpace::String="ADFGLTVXYZ", alphabetPool::String="abcde...")

Handles creation of chart based on a given codespace and alphabet. Assumes standard if none given. """

function createChart(chartLoc::String; codeSpace::String="ADFGLTVXYZ", alphabetPool::String="abcdefghiklmnopqrstuvwxyzj0123456789{}()[],.:;+-=*/@%_~!?#^&'<>|ABCDEFGHIJKLMNOPQRSTUVWXYZ `\u0000\u000D\u0020\u005C\u0024\u0022\u0009\u000A")
	chart=Dict()
	charArray=split(alphabetPool, "")
	for i in 1:length(codeSpace)
		for j in 0:length(codeSpace)-1
			randSpot=rand(Int8)%length(charArray)
			randSpot=(randSpot==0?1:(randSpot<0?randSpot*-1:randSpot))
			chart[charArray[randSpot][1]] = string(codeSpace[i])*string(codeSpace[j+1])
			deleteat!(charArray, randSpot)
		end
	end
	saveChart(chart, chartLoc)
	return chart
end

""" encryptFile(fileName::String, passPhrase::String, chart::String)

Handles encryption of file-loads a file from given location, uses given phrase, loads a chart from given location. """

function encryptFile(fileName::String, passPhrase::String, chart::String)
	chart = loadChart(chart) #catch an error if no chart(?)
	#encrypt phase
	preString=""
	inFile = open(fileName,"r")
	fname=split(fileName, '.')\
	fnameOut=fname[1]*".adfgvx"*"."*fname[end]
	outFile = open(fnameOut, "w")
	inString = readstring(inFile)
	for char in inString
		preString*=get(chart, char, 'x')
	end
	fillAmt=length(preString)%length(passPhrase)
	if fillAmt > 0
		for i in 1:(length(passPhrase)-fillAmt)
			preString*=chart['#']
		end
	end

	#phraseShifting phase- take phrase, arrange letters below, sort by phrase alphabetically
	#-turn code letters into array with len(phrase) cols and however many rows
	#-sort array, carry the letters array along with it.
	#=
	codeArray=fill(' ', Int(ceil(length(preString)/length(passPhrase))), length(passPhrase))
	#codeArray=[ceil(length(outString)/length(prePhrase))][length(prePhrase)]
	k=1
	for i in 1:Int(ceil(length(preString)/length(passPhrase)))
		for j in 1:length(passPhrase)
			codeArray[i, j]=preString[k]
			k+=1
		end
	end
	#sortedString=phraseSort(codeArray, passPhrase, true)
	=#
	#temp fix/bypass
	sortedString=preString

	#output phase
	outString=""
	spacing=0; lines=0
	for char in sortedString
		if spacing == 8
			outString*=" "
			spacing=0
			lines+=1
		end
		if lines == 8
			outString*="\n"
			lines=0
		end
		outString*=char
		spacing+=1
	end
	print(outFile, outString)
	close(outFile); close(inFile)
end

""" decryptFile(fileName::String, passPhrase::String, chartLoc::String)

Handles decryption of file- given location of file (.adfgvx.ext), given phrase, and given the chart location. """

function decryptFile(fileName::String, passPhrase::String, chartLoc::String)
	chart=loadChart(chartLoc)
	#reverse the Dict
	decodeChart=Dict()
	for(key, value) in chart
		decodeChart[value]=key
	end
	#reading from file
	fname=split(fileName, '.')
	fnameOut = fname[1]*".test"*"."*fname[end]
	outFile = open(fnameOut, "w")
	spacing=0; lines=0
	inFile = open(fileName,"r")
	fileString = readstring(inFile)

	#phraseshifting-
	#take filestring, unshift it according to phrase, carry on decoding
	fileArray = split(fileString)
	codeArray=fill(' ', length(passPhrase), Int(ceil(length(fileString)/length(passPhrase))))
	#codeArray=[ceil(length(outString)/length(prePhrase))][length(prePhrase)]
	#=
	k=1
	for i in 1:Int(ceil(length(fileString)/length(passPhrase)))
		for j in 1:length(passPhrase)
			codeArray[j, i]=fileArray[k]
			k+=1
		end
	end
	sortedArray=phraseSort(codeArray, passPhrase, false)
	=#

	#output phase
	#can do this better- allow for custom block sizes, and just grab 2 chars at a time, throw them at the chart, keep on going
	for block in fileArray
		for char in 1:2:length(block)
			code = string(block[char])*string(block[char+1])
			print(outFile, get(decodeChart, code, 'x'))
		end
	end
	close(outFile); close(inFile)
end

""" phraseSort(message::Any, passPhrase::String, encoding::Bool)

Handles sorting the entire 'message' by a given phrase. also takes a flag to decide if shifting for encoding or decoding. """

function phraseSort(message::Any, passPhrase::String, encoding::Bool)
	postPhrase=[]
	prePhrase=split(passPhrase,"")
	stringOut=message #just because i dont want to do this YET
	if encoding == true
		#encoding
		#sort the phrase into alphabetical order, and drag the array 'rows' with you
		#for loops? per char in that char of phrase? seems slow-
		#get offset library - phrase => aehprs == [3,1,2,-3,1,-4], and just apply the offsets when re-assigning chars to new/duped array
		#grab chars in phrase, and math out their ascii offsets? then simply put all chars in newArray[j+offset[j], i]

		#return as finished string
	else
		#decoding
		#like above, but do opposite- subtract offsets when putting in out array- keep single by putting in newArray[i+(j-1)*len(phrase)]

		#return as a single dimensional code char array
	end

	return stringOut #prolly not stringout
end

""" testModule()

Runs a very basic test of the programs features- creates chart, encrypts, decrypts. """

function testModule()
	fname = ""
	dname = ""
	#chart = createChart("chart.jld")
	#print(chart)
	encryptFile(fname, "phrase", "chart.jld")
	decryptFile(dname, "phrase", "chart.jld")
end

end


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools, Random
Random.seed!(1)

u0 = zeros(2)
v0 = ones(2)
γ = 1

f1_harmonic(v,u,p,t) = -u
f2_harmonic(v,u,p,t) = v
g(u,p,t) = 1

f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t)
f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t)
g_iip(du,u,p,t) = du .= g(u,p,t)

ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0))
sol1 = solve(prob1,BAOAB(gamma=γ);dt=1/10,save_noise=true)

prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
sol2 = solve(prob2,BAOAB(gamma=γ);dt=1/10)

@test sol1[:] ≈ sol2[:]

dts = (1/2) .^ (8:-1:4)

# Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v.
sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
@test abs(sim1.𝒪est[:weak_final]-1) < 0.3


text:
################################################################################
#                           TEST: Compose Rotations
################################################################################

for i = 1:samples
    # DCMs
    # ==========================================================================

    # Sample 4 DCMs.
    for sym in [:D1, :D2, :D3, :D4]
        @eval $sym = angle_to_dcm(EulerAngles(-pi + 2*pi*rand(),
                                              -pi + 2*pi*rand(),
                                              -pi + 2*pi*rand(),
                                              rand($rot_seq_array)))
    end

    # Test the function `compose_rotation`.
    @test          D1 === compose_rotation(D1)
    @test       D2*D1  ≈  compose_rotation(D1,D2)       atol=1e-14
    @test    D3*D2*D1  ≈  compose_rotation(D1,D2,D3)    atol=1e-14
    @test D4*D3*D2*D1  ≈  compose_rotation(D1,D2,D3,D4) atol=1e-14

    # Euler angle and axis
    # ==========================================================================

    # Sample 4 Euler angle and axis.
    for sym in [:ea1, :ea2, :ea3, :ea4]
        α   = -π + 2π*rand()
        v   = randn(3)
        v   = v/norm(v)

        @eval $sym = EulerAngleAxis($α,$v)
    end

    ear1 = ea2*ea1
    eac1 = compose_rotation(ea1,ea2)

    ear2 = ea3*ea2*ea1
    eac2 = compose_rotation(ea1,ea2,ea3)

    ear3 = ea4*ea3*ea2*ea1
    eac3 = compose_rotation(ea1,ea2,ea3,ea4)

    # Test the function `compose_rotation`.
    @test ea1 === compose_rotation(ea1)
    @test ear1 ≈ eac1 atol = 1e-14
    @test ear2 ≈ eac2 atol = 1e-14
    @test ear3 ≈ eac3 atol = 1e-14

    # Euler angles
    # ==========================================================================

    # Sample 4 Euler angles.
    for sym in [:Θ1, :Θ2, :Θ3, :Θ4]
        θx      = -π + 2π*rand()
        θy      = -π + 2π*rand()
        θz      = -π + 2π*rand()
        rot_seq = Meta.quot(rand([:XYZ,:ZYX]))

        @eval $sym = EulerAngles($θx, $θy, $θz, $rot_seq)
    end

    Θr1 = Θ2*Θ1
    Θc1 = compose_rotation(Θ1,Θ2)

    Θr2 = Θ3*Θ2*Θ1
    Θc2 = compose_rotation(Θ1,Θ2,Θ3)

    Θr3 = Θ4*Θ3*Θ2*Θ1
    Θc3 = compose_rotation(Θ1,Θ2,Θ3,Θ4)

    # Test the function `compose_rotation`.
    @test Θ1 === compose_rotation(Θ1)
    @test Θr1 ≈ Θc1 atol = 1e-14
    @test Θr2 ≈ Θc2 atol = 1e-14
    @test Θr3 ≈ Θc3 atol = 1e-14

    # Quaternions
    # ==========================================================================

    # Sample 4 quaternions.
    for sym in [:q1, :q2, :q3, :q4]
        @eval $sym = angle_to_quat(EulerAngles(-pi + 2*pi*rand(),
                                               -pi + 2*pi*rand(),
                                               -pi + 2*pi*rand(),
                                               rand($rot_seq_array)))
    end

    # Test the function `compose_rotation`.
    @test           q1     === compose_rotation(q1)
    @test       (q1*q2)[:]  ≈  compose_rotation(q1,q2)[:]       atol=1e-14
    @test    (q1*q2*q3)[:]  ≈  compose_rotation(q1,q2,q3)[:]    atol=1e-14
    @test (q1*q2*q3*q4)[:]  ≈  compose_rotation(q1,q2,q3,q4)[:] atol=1e-14
end


max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# SPDX-License-Identifier: X11
# 2020-11-20
# Wandering (400pt)

function main()
  N = parse(Int, readline())
  A = parse.(Int, (split ∘ readline)())

  step = 0
  curr = 0
  mstep = 0
  m = 0
  for x ∈ A
    step += x
    mstep = max(mstep, step)
    m = max(m, curr + mstep)
    curr += step
  end

  println(m)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using HierarchicalUtils: encode, stringify

"""
    AbstractMillModel

Supertype for any model defined in `Mill.jl`.
"""
abstract type AbstractMillModel end

const MillFunction = Union{Dense, Chain, Function}

_make_array_model(m::Union{MillFunction, Function}) = ArrayModel(m)
_make_array_model(m) = m

include("arraymodel.jl")
include("bagmodel.jl")
include("productmodel.jl")
include("lazymodel.jl")

function Base.show(io::IO, @nospecialize(m::AbstractMillModel))
    print(io, nameof(typeof(m)))
    if !get(io, :compact, false)
        _show_submodels(io, m)
    end
end

_show_submodels(io, m::ArrayModel) = print(io, "(", m.m, ")")
_show_submodels(io, m::BagModel) = print(io, " … ↦ ", m.a, " ↦ ", m.bm)
_show_submodels(io, m::ProductModel) = print(io, " … ↦ ", m.m)
_show_submodels(io, m::LazyModel{Name}) where {Name} = print(io, "{", Name, "}")
_show_submodels(io, _) = print(io)

include("reflectinmodel.jl")


text:
using Flux, Flux.Data.MNIST, LinearAlgebra
using Flux: onehotbatch

imgs   = Flux.Data.MNIST.images()
labels = Flux.Data.MNIST.labels()
nTrain = length(imgs)
trainData = vcat([hcat(float.(imgs[i])...) for i in 1:nTrain]...)
trainLabels = labels[1:nTrain]
testImgs = Flux.Data.MNIST.images(:test)
testLabels = Flux.Data.MNIST.labels(:test)
nTest = length(testImgs)
testData = vcat([hcat(float.(testImgs[i])...) for i in 1:nTest]...)

A = [ones(nTrain) trainData]
Adag = pinv(A)
tfPM(x) = x ? +1 : -1
yDat(k) = tfPM.(onehotbatch(trainLabels,0:9)'[:,k+1])
bets = [Adag*yDat(k) for k in 0:9]

classify(input) = findmax([([1 ; input])'*bets[k] for k in 1:10])[2]-1

predictions = [classify(testData[k,:]) for k in 1:nTest]
confusionMatrix = [sum((predictions .== i) .& (testLabels .== j))
				for i in 0:9, j in 0:9]
accuracy = sum(diag(confusionMatrix))/nTest

println("Accuracy: ", accuracy, "\nConfusion Matrix:")
show(stdout, "text/plain", confusionMatrix)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, Test, DiffEqBase, DiffEqCallbacks

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

prob = prob_ode_linear

cb = AutoAbstol()
integrator = init(prob,BS3(),callback=cb)
at1 = integrator.opts.abstol
step!(integrator)
at2 = integrator.opts.abstol
@test at1 < at2
step!(integrator)
at3 = integrator.opts.abstol
@test at2 < at3

solve(prob,BS3(),callback=cb)

prob2D = prob_ode_2Dlinear
cb2D = AutoAbstol(init_curmax=zeros(4,2))
integrator2D_1 = init(prob2D,BS3(),callback=cb,abstol=1e-6)
integrator2D_2 = init(prob2D,BS3(),callback=cb,abstol=fill(1e-6,4,2))
@test all(integrator2D_2.opts.abstol .== integrator2D_1.opts.abstol)
step!(integrator2D_1)
step!(integrator2D_2)
@test all(integrator2D_2.opts.abstol .== integrator2D_1.opts.abstol)
step!(integrator2D_1)
step!(integrator2D_2)
@test all(integrator2D_2.opts.abstol .== integrator2D_1.opts.abstol)

sol1 = solve(prob2D,BS3(),callback=cb2D,abstol=fill(1e-6,4,2))
sol2 = solve(prob2D,BS3(),callback=cb2D,abstol=fill(1e-6,4,2),reltol=fill(1e-3,4,2))
@test sol1.t == sol2.t && sol1.u == sol2.u
@test_throws MethodError solve(prob,BS3(),callback=cb2D,abstol=1e-6)


text:
# This file includes modified source code from https://github.com/odow/SDDP.jl
# as at 5b3ba3910f6347765708dd6e7058e2fcf7d13ae5

#  Copyright 2017, Oscar Dowson and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################

uberror(x) =  error("You must provide an upper bound on $(x).
    The lower the bound, the smaller the statespace.")

function anonymousstate!(sp::JuMP.Model, nstates::Int, initial::Vector{<:Real},
                             nameout::AbstractString, namein::AbstractString)
    statein  = @variable(sp, [i=1:nstates], basename="_bin_$namein", start=initial[i])
    stateout = @variable(sp,   [1:nstates], basename="_bin_$nameout")
    SDDP.statevariable!(sp, statein, stateout)
    stateout, statein
end

function binarystate!(sp::JuMP.Model, xout::JuMP.Variable, xin::JuMP.Variable)

    ub = sp.colUpper[xout.col]
    isinf(ub) && uberror(xout)

    nstates = bitsrequired(Int(ub))
    init = binexpand(Int(getvalue(xin)), length=nstates)

    v, v0 = anonymousstate!(sp, nstates, init, JuMP.getname(xout), JuMP.getname(xin))
    @constraint(sp, xout == bincontract(v))
    @constraint(sp, xin  == bincontract(v0))
end

function binarystate!(sp::JuMP.Model, xout::JuMP.Variable, xin::JuMP.Variable, eps::Float64)

    ub = sp.colUpper[xout.col]
    isinf(ub) && uberror(xout)

    nstates = bitsrequired(ub, eps)
    init = binexpand(float(getvalue(xin)), eps, length=nstates)

    v, v0 = anonymousstate!(sp, nstates, init, JuMP.getname(xout), JuMP.getname(xin))
    @constraint(sp, xout == bincontract(Float64, v, eps))
    @constraint(sp, xin  == bincontract(Float64, v0, eps))
end

function binarystate!(sp::JuMP.Model, xout::Array{JuMP.Variable}, xin::Array{JuMP.Variable})
    @assert size(xout) == size(xin)
    for i = 1:length(xin)
        binarystate!(sp, xout[i], xin[i])
    end
end

function binarystate!(sp::JuMP.Model, xout::Array{JuMP.Variable}, xin::Array{JuMP.Variable}, eps::Float64)
    @assert size(xout) == size(xin)
    for i = 1:length(xin)
        binarystate!(sp, xout[i], xin[i], eps)
    end
end

function binarystate!(sp::JuMP.Model, xout::JuMP.JuMPArray, xin::JuMP.JuMPArray)
    @assert length(keys(xin)) == length(keys(xout))
    for key in keys(xin)
        binarystate!(sp, xout[key...], xin[key...])
    end
end

function binarystate!(sp::JuMP.Model, xout::JuMP.JuMPArray, xin::JuMP.JuMPArray, eps::Float64)
    @assert length(keys(xin)) == length(keys(xout))
    for key in keys(xin)
        binarystate!(sp, xout[key...], xin[key...], eps)
    end
end

"""
    @binarystate(sp, stateleaving, stateentering, vartype)

# Description

Define a new state variable in the SDDiP subproblem `sp`.

# Arguments
 * `sp`               the subproblem
 * `stateleaving`     any valid JuMP `@variable` syntax to define the value of the state variable at the end of the stage
 * `stateentering`    any valid JuMP `@variable` syntax to define the value of the state variable at the beginning of the stage
 *  `vartype`         type of variable: choices are Bin, Int, or Cont
 *  `epsilon`         the accuracy of binary expansion for continuous variables (default is always 0.1 if not specified)

# Examples
 * `@binarystate(sp, 0 <= x[i=1:3] <= 1,  x0==rand([0 1], 3)[i], Bin )`
 * `@binarystate(sp,      y        <= 1,  y0==0.5              , Cont, 0.01)`
 * `@binarystate(sp,      z        <= 10, z0==5                , Int )`
"""
macro binarystate(sp, x, x0, vtype_args...)

    # Model
    sp = esc(sp)

    # eps and type
    if length(vtype_args) == 0
        vtype = :(Bin)
    elseif length(vtype_args) == 1
        vtype = vtype_args[1]
        (vtype == :Cont) && (eps = 0.1)
    elseif length(vtype_args) == 2
        vtype = vtype_args[1]
        @assert vtype == :Cont
        eps = vtype_args[2]
    else
        error("Type of variable $(x) not properly specified in @binarystate.")
    end

    # If it's already a binary variable, just call @state
    if vtype == :Bin
        return quote
            # Ordinary @state call
            $(Expr(:macrocall, Symbol("@state"), sp, esc(x), esc(x0)))
        end

    # Otherwise prepare for reformulation
    elseif vtype in [:Int, :Cont]

        # Get variables exactly as in SDDP.jl
        @assert x0.head == :call && x0.args[1] == :(==) # must be ==
        _, symin, init = x0.args
        if SDDP.is_comparison(x)
            if length(x.args) == 5
                xin = SDDP._copy(x.args[3])
            elseif length(x.args) == 3
                xin = SDDP._copy(x.args[2])
            else
                error("Unknown format for $(x)")
            end
        else
            xin = SDDP._copy(x)
        end
        if isa(xin, Expr)                   # x has indices
            xin.args[1] = symin             # so just change the name
        else                                # its just a Symbol
            xin = symin                     # so change the Symbol
        end

        # The names the user chose will be dummy variables in the model
        varin, varout = gensym(), gensym()

        code = quote
            $varout = $(Expr(:macrocall, Symbol("@variable"), sp, esc(x)))
            $varin  = $(Expr(:macrocall, Symbol("@variable"), sp, esc(xin), Expr(SDDP.KW_SYM, SDDP.START, esc(init))))
        end

        if vtype == :Int
            # Integer case
            push!(code.args, quote
                binarystate!($sp, $varout, $varin)
            end)
        else
            # Continuous case
            push!(code.args, quote
                binarystate!($sp, $varout, $varin, $(esc(eps)))
            end)
        end

        push!(code.args, quote
            $varout, $varin
        end)
        return code

    else
        error("You must specify a variable type for $(x) because you are calling @binarystate.
                Choices are Bin, Int, or Cont.")
    end # vtype
end # function

"""
    @binarystates(sp, begin
        stateleaving1, stateentering1, vartype
        stateleaving2, stateentering2, vartype
    end)
# Description
Define a new state variables in the subproblem `sp`.
# Arguments
* `sp`               the subproblem
* `stateleaving`     any valid JuMP `@variable` syntax to define the value of the state variable at the end of the stage
* `stateentering`    any valid JuMP `@variable` syntax to define the value of the state variable at the beginning of the stage
*  `vartype`         type of variable: choices are Bin, Int, or Cont
*  `epsilon`         the accuracy of binary expansion for continuous variables (default is always 0.1 if not specified)

# Usage
    @binarystates(sp, begin
        0 <= x[i=1:3] <= 1,  x0==rand([0 1], 3)[i], Bin
             y        <= 1,  y0==0.5              , Cont, 0.01
             z        <= 10, z0==5                , Int
     end)
"""
macro binarystates(m, b)
    @assert b.head == :block || error("Invalid syntax for @binarystates.")
    code = quote end
    for line in b.args
        if !Base.Meta.isexpr(line, :line)
            if line.head == :tuple && 3 <= length(line.args) <= 4
                escargs = []
                for ex in line.args
                    push!(escargs, esc(ex))
                end
                push!(code.args, Expr(:macrocall, Symbol("@binarystate"), esc(m), escargs...))
            else
                error("Unknown arguments in @binarystates.")
            end
        end
    end
    push!(code.args, :(nothing))
    return code
end


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this        = parse(Int64, this)

#
# Load packages and functions
include("../src/init.jl")
include("../src/fct.jl")
include("../src/worker_fct.jl")

#
# Output paths and filter settings
plots_base_path = "../../waveforms/sim/plots/"
base_path_raw   = "../../waveforms/sim/raw_wf/"
base_path       = "../../waveforms/sim/wf/"
plots_bl_path   = "../../waveforms/baselines/plots/"
base_bl_path    = "../../waveforms/baselines/wf/"
# base_bl_path_filter = "../../waveforms/baselines/wf_filter/"
sampling_time   = 1u"ns";
parameters = JSON.parsefile("../dicts/electronics_parameters.json")

channels = []
for c in 0:1:36
    files = glob(base_path_raw * lpad(c, 2, "0") * "-" * channel_to_name[c] * "/w_filter/raw*/*.h5");
    if length(files) > 0 && channel_to_bege[c]
        push!(channels, c)
    end
end
channels = get_share_for_worker(channels, num_workers, this);

for ch in channels
#     IJulia.clear_output(true)
    Base.run(`clear`)
    str_ch = lpad(ch, 2, "0");

    @info("Start Ch$str_ch | " * channel_to_name[ch])
    @info(">---------------------------------------------------------------<")
    @info("Read in baselines")
    bl_files  = glob(joinpath(base_bl_path, "*/" * str_ch*"-"*channel_to_name[ch]*"/*.h5"));

    bl = Table(run = [], datasetID = [], waveform = [])
    @showprogress 1 "Collecting baselines for Ch$str_ch ..." for file in bl_files
        tmp_data = HDF5.h5open(file, "r") do h5f
            LegendHDF5IO.readdata(h5f, "data")
        end
        if typeof(tmp_data.run) == Array{Int64,1}
            append!(bl, Table(run = tmp_data.run, datasetID = tmp_data.datasetID, waveform = tmp_data.waveform))
        else
            println(typeof(tmp_data.run))
        end
    end
    bl_dict = JSON.parsefile("../dicts/bl_rms_limits.json");
    bl = bl |> @filter(bl_dict[channel_to_name[ch]]["low_rms_limit"] < rms(_.waveform.value) < bl_dict[channel_to_name[ch]]["high_rms_limit"]) |> Table

    @info("Number of baselines: " * string(size(bl,1)))
    @info(">---------------------------------------------------------------<")
    files = glob(joinpath(base_path_raw, str_ch * "-" * channel_to_name[ch] * "/w_filter/*/*.h5"));
    number_of_pulses = 0;
    key_str = []
    for file in files
        push!(key_str, split(basename(file), ".root_evts_")[1])
    end
    key_str = unique(key_str)
    @showprogress 1 "Counting simulated waveforms..." for k in key_str
        last_file = basename(glob(joinpath(base_path_raw, str_ch * "-" * channel_to_name[ch] * "/w_filter/*/" * String(k) * "*.h5"))[end])
        number_of_pulses += parse(Int64, split(split(last_file, ".h5")[1], "-")[end])
    end

    @info("Found number of simulation pulses: " * string(number_of_pulses))
    @info("Found number of noise baselines: " * string(size(bl,1)))

    # get first_range
    if number_of_pulses <= size(bl,1)
        first_range = [1]
    else
        mult_factor = ceil(number_of_pulses/size(bl,1))
        first = 1
        last = 201
        step = Int(round((last-first)/(mult_factor)))
        first_range = first:step:last
    end

    @info("This results in "*string( size(bl, 1) *length(first_range) )*" baselines for "*string(number_of_pulses)*" pulses")
    @info(">---------------------------------------------------------------<")
    @info("Indexing the baselines")
    indices = Table(bl=[], first=[])
    for i in 0:1:number_of_pulses
        bl_id = i%size(bl,1) + 1
        first = first_range[Int(ceil((i + 1) / size(bl,1)))]
        append!(indices, Table(bl=[bl_id], first=[first]))
    end
    pulse_id = 1
    @info("Add baselines to pulses")

    @showprogress 1 "Adding baselines..." for file in files
        create_file = false
        ch_str = lpad(ch, 2, "0");
        filename = joinpath(
            joinpath(
                joinpath(base_path,
                    str_ch * "-" * channel_to_name[ch]),
                split(basename(file), ".root")[1]),
            basename(file)
        )
        !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""

        if isfile(filename)
            if stat(filename).size <= stat(file).size/2
                create_file = true
            end
        else
            create_file = true
        end
        if create_file
            data_w_bl = Table(energy=[], multiplicity=[], timestamp=[], run=[], channel=[], waveform=[])

            data = HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "data")
            end
            for event in data
                time = event.waveform.time
                pulse = event.waveform.value ./ get_avg_maximum(event.waveform.value, 10)
                pulse .*= sum(event.energy) * 1000 # Energy is stored in MeV

                bl_id = indices[pulse_id].bl
                first = indices[pulse_id].first

                last = first + length(pulse) -1
                pulse .+= bl[bl_id].waveform.value[first:1:last]

                append!(data_w_bl, Table(energy=[event.energy],
                                        multiplicity=[event.multiplicity],
                                        timestamp=[event.timestamp],
                                        run=[event.run],
                                        channel=[event.channel],
                                        waveform=[RDWaveform(time, pulse)]));
                pulse_id += 1
            end
            HDF5.h5open(filename, "w") do h5f
                LegendHDF5IO.writedata(h5f, "data", Table(energy        = VectorOfArrays(Array{typeof(data_w_bl[1].energy),1}(data_w_bl.energy)),
                                                        multiplicity = Int32.(data_w_bl.multiplicity),
                                                        timestamp    = Int32.(data_w_bl.timestamp),
                                                        run          = Int32.(data_w_bl.run),
                                                        channel      = Int32.(data_w_bl.channel),
                                                        waveform     = ArrayOfRDWaveforms( Array{typeof(data_w_bl[1].waveform), 1}(data_w_bl.waveform) )))
            end
        else
            pulse_id += size(HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "data")
            end, 1)

        end
    end
end


text:
### Types

# Represents the number to fill into a cell, i.e. 1 to 9.
typealias Number Int

# A subset of {1, 2, ..., 9}, which represents the leftover number that can be use to fill the corresponding cells.
immutable LeftoverNumbers
    val_bits::Uint16
end

# Represents a cell.
# In natural order, (1, 1) is the top left, (9, 1) the bottom left, and (9, 9) the bottom right.
immutable Cell
    row::Int
    col::Int
end

# Represent the state of all the cells
immutable CellsState
    data::Matrix{Number}
end

# Represent the state of components like the blocks, rows, and columns.
# The function `findComponentIndex` is used to find the index in the component of the given cell.
immutable ComponentState
    data::Vector{LeftoverNumbers}
    findComponentIndex::Function
end

# Represents the solution sequence.
# By the end of the program, this sequence shows the order of filling in of number for the last found correct solution.
immutable SequenceState
    data::Vector{Cell}
end

# Represents the number of times we change the values with that amount of entries filled.
immutable LevelCount
    data::Vector{Int}
end

# The state of the game currently.
# This will evolve as the Sudoku puzzle is being solved.
type GameState
    cells::CellsState

    blocks::ComponentState
    rows::ComponentState
    cols::ComponentState

    sequence::SequenceState

    solutions::Vector{Matrix{Number}}
    levelCount::LevelCount
end

### Constructors

# Each cell is left blank initially.
CellsState() = CellsState([BLANK for i in 1:9, j in 1:9])

# Each leftover numbers in the component is initialize with the full set {1, 2,..., 9} inside.
ComponentState(findComponentIndex::Function) = ComponentState([LeftoverNumbers(ONES) for i in 1:9], findComponentIndex)

# We allocate cells to the sequence in the order (1, 1), ..., (1, 9), (2, 1), ..., ..., (9, 9).
function SequenceState()
    data = Array(Cell, 81)
    index = 1
    for i in 1:9, j in 1:9
        data[index] = Cell(i, j)
        index += 1
    end
    SequenceState(data)
end

LevelCount() = LevelCount([zero(Int) for i in 1:81])

# Constructing `GameState` is as simple as calling constructor for all its components.
function GameState()
    return GameState(CellsState(),

                     ComponentState(getBlock),
                     ComponentState(getRow),
                     ComponentState(getCol),

                     SequenceState(),

                     Array(Matrix{Number}, 0),
                     LevelCount())
end

max_stars_repo_name: Teo-ShaoWei/SudokuSolverSG

max_issues_repo_name: Teo-ShaoWei/SudokuSolverSG

max_forks_repo_name: Teo-ShaoWei/SudokuSolverSG


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Distributions, SpecialFunctions, QuadGK, Plots
pyplot(fmt = :svg)

α = range(0.1, 10; length=101)
# θ = 2

EX(α) = 2α
ElogX(α) = log(2) + digamma(α) # digamma(x) = ψ(x) = d(log Γ(x))/dx

F(α) = quadgk(x -> x*pdf(Gamma(α, 2.0), x), 0, Inf)[1]
G(α) = quadgk(x -> log(x)*pdf(Gamma(α, 2.0), x), 0, Inf)[1]

y = EX.(α)
z = ElogX.(α)
y_numint = F.(α)
z_numint = G.(α)

P = plot(α, y; label="E[X|α]", xlabel="α", legend=:topleft)
plot!(α, y_numint; label="num. int.", ls=:dash)

Q = plot(α, z; label="E[log X|α]", xlabel="α", legend=:bottomright)
plot!(α, z_numint; label="num. int.", ls=:dash)

plot(P, Q; size=(800, 300))

# %%

max_stars_repo_path: 0007/gamma dist.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0007/gamma dist.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0007/gamma dist.jl

max_forks_repo_name: genkuroki/public


text:
"""
    Node(id::Int, MTG<:AbstractNodeMTG, attributes)
    Node(name::String, id::Int, MTG<:AbstractNodeMTG, attributes)
    Node(id::Int, parent::Node, MTG<:AbstractNodeMTG, attributes)
    Node(name::String, id::Int, parent::Node, MTG<:AbstractNodeMTG, attributes)
    Node(
        name::String,
        id::Int,
        parent::Node,
        children::Union{Nothing,Dict{Int,Node}},
        siblings::Union{Nothing,Dict{Int,Node}},
        MTG<:AbstractNodeMTG,
        attributes
    )

Type that defines an MTG node (*i.e.* an element) with the name of the node, its unique id,
its parent (only if its not the root node), children, MTG encoding (see [`NodeMTG`](@ref) or
[`MutableNodeMTG`](@ref)) and attributes.

The node is an entry point to a Mutli-Scale Tree Graph, meaning we can move through the MTG from any
of its node. The root node is the node without parent. A leaf node is a node without any children.
Root and leaf nodes are used with their computer science meaning throughout the package, not in the
biological sense.
"""
mutable struct Node{T<:AbstractNodeMTG,A}
    "Name of the node. Should be unique in the MTG."
    name::String
    "Node unique ID"
    id::Int
    "Parent node."
    parent::Union{Nothing,Node}
    "Dictionary of children nodes, or Nothing if no children."
    children::Union{Nothing,Dict{Int,Node}}
    "Dictionary of sibling(s) nodes if any, or else Nothing. Can be Nothing if not computed too."
    siblings::Union{Nothing,Dict{Int,Node}}
    "MTG encoding (see [`NodeMTG`](@ref) or [`MutableNodeMTG`](@ref))."
    MTG::T
    "Node attributes. Can be anything really."
    attributes::A
end

# Shorter way of instantiating a Node:

# - for the root:
Node(name::String, id::Int, MTG::T, attributes) where {T<:AbstractNodeMTG} = Node(name, id, nothing, nothing, nothing, MTG, attributes)

# If the name is not given, we compute one from the id:
Node(id::Int, MTG::T, attributes) where {T<:AbstractNodeMTG} = Node(join(["node_", id]), id, MTG, attributes)
# If the id is not given, it is the root node, so we use 1
Node(MTG::T, attributes) where {T<:AbstractNodeMTG} = Node(1, MTG, attributes)

# Special case for the NamedTuple and MutableNamedTuple, else it overspecializes and we
# can't mutate attributes, i.e. we get somthing like
# Node{NodeMTG,MutableNamedTuple{(:a,), Tuple{Base.RefValue{Int64}}}} instead of just:
# Node{NodeMTG,MutableNamedTuple}
function Node(name::String, id::Int, MTG::M, attributes::T) where {M<:AbstractNodeMTG,T<:MutableNamedTuple}
    Node{typeof(MTG),MutableNamedTuple}(name, id, nothing, nothing, nothing, MTG, attributes)
end

function Node(name::String, id::Int, MTG::M, attributes::T) where {M<:AbstractNodeMTG,T<:NamedTuple}
    Node{typeof(MTG),NamedTuple}(name, id, nothing, nothing, nothing, MTG, attributes)
end

# - for all others:
function Node(name::String, id::Int, parent::Node, MTG::M, attributes) where {M<:AbstractNodeMTG}
    node = Node(name, id, parent, nothing, nothing, MTG, attributes)
    addchild!(parent, node)
    return node
end

# Idem for MutableNamedTuple here:
function Node(name::String, id::Int, parent::Node, MTG::M, attributes::T) where {M<:AbstractNodeMTG,T<:MutableNamedTuple}
    node = Node{typeof(MTG),MutableNamedTuple}(name, id, parent, nothing, nothing, MTG, attributes)
    addchild!(parent, node)
    return node
end

# Idem, if the name is not given, we compute one from the id:
Node(id::Int, parent::Node, MTG::T, attributes) where {T<:AbstractNodeMTG} = Node(join(["node_", id]), id, parent, MTG, attributes)
# If the id is not given, it is the root node, so we use 1
Node(parent::Node, MTG::T, attributes) where {T<:AbstractNodeMTG} = Node(new_id(get_root(parent)), parent, MTG, attributes)

#  Next lines are adapted from either:
# <https://github.com/JuliaCollections/AbstractTrees.jl>
# <https://github.com/dellison/ConstituencyTrees.jl/blob/master/src/trees.jl>
# <https://github.com/vh-d/DataTrees.jl/blob/master/src/indexing.jl>

Base.eltype(::Type{<:TreeIterator{Node{T,D}}}) where {T,D} = Node{T,D}
Base.IteratorEltype(::Type{<:TreeIterator{Node{T,D}}}) where {T,D} = Base.HasEltype()

# Help Julia infer what's inside a Node when doing iteration (another node)
Base.eltype(::Type{Node{T,D}}) where {T,D} = Node{T,D}

# Iteration over the immediate children:
function Base.iterate(node::T) where {T<:Node}
    isleaf(node) ? nothing : (node[1], 1)
end

function Base.iterate(node::T, state::Int) where {T<:Node}
    state += 1
    state > length(children(node)) ? nothing : (node[state], state)
end

Base.IteratorSize(::Type{Node{T,D}}) where {T,D} = Base.SizeUnknown()

## Things we need to define to leverage the native iterator from AbstractTrees over children

# Set the traits of this kind of tree
AbstractTrees.parentlinks(::Type{Node{T,D}}) where {T,D} = AbstractTrees.StoredParents()
AbstractTrees.siblinglinks(::Type{Node{T,D}}) where {T,D} = AbstractTrees.StoredSiblings()
AbstractTrees.children(node::Node) = MultiScaleTreeGraph.children(node)
AbstractTrees.nodetype(::Node) = Node

Base.parent(node::Node) = isdefined(node, :parent) ? node.parent : nothing
Base.parent(root::Node, node::Node) = isdefined(node, :parent) ? node.parent : nothing

function AbstractTrees.nextsibling(tree::Node, child::Node)
    MultiScaleTreeGraph.nextsibling(child)
end

# We also need `pairs` to return something sensible.
# If you don't like integer keys, you could do, e.g.,
#   Base.pairs(node::BinaryNode) = BinaryNodePairs(node)
# and have its iteration return, e.g., `:left=>node.left` and `:right=>node.right` when defined.
# But the following is easy:
Base.pairs(node::Node) = enumerate(node)


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Qualification Round Africa 2010
# Problem B: Reverse Words

# Reverses the order of words in a string.
revwords(str) = join(reverse(split(str)), " ")

function main()
	n = int(readline())
	for i = 1:n
		str = revwords(readline())
		println("Case #$i: ", str)
	end
end

main()


text:
"""
$(TYPEDEF)

A system of partial differential equations.

# Fields
$(FIELDS)

# Example

```julia
using ModelingToolkit

@parameters x
@variables t u(..)
Dxx = Differential(x)^2
Dtt = Differential(t)^2
Dt = Differential(t)

#2D PDE
C=1
eq  = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))

# Initial and boundary conditions
bcs = [u(t,0) ~ 0.,# for all t > 0
       u(t,1) ~ 0.,# for all t > 0
       u(0,x) ~ x*(1. - x), #for all 0 < x < 1
       Dt(u(0,x)) ~ 0. ] #for all  0 < x < 1]

# Space and time domains
domains = [t ∈ (0.0,1.0),
           x ∈ (0.0,1.0)]

@named pde_system = PDESystem(eq,bcs,domains,[t,x],[u])
```
"""
struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
    "The equations which define the PDE"
    eqs::Any
    "The boundary conditions"
    bcs::Any
    "The domain for the independent variables."
    domain::Any
    "The independent variables"
    ivs::Any
    "The dependent variables"
    dvs::Any
    "The parameters"
    ps::Any
    """
    defaults: The default values to use when initial conditions and/or
    parameters are not supplied in `ODEProblem`.
    """
    defaults::Dict
    """
    type: type of the system
    """
    connector_type::Any
    """
    systems: The internal systems. These are required to have unique names.
    """
    systems::Vector
    """
    name: the name of the system
    """
    name::Symbol
    @add_kwonly function PDESystem(eqs, bcs, domain, ivs, dvs,
                                   ps = SciMLBase.NullParameters();
                                   defaults = Dict(),
                                   systems = [],
                                   connector_type = nothing,
                                   checks::Bool = true,
                                   name)
        if checks
            all_dimensionless([dvs; ivs; ps]) || check_units(eqs)
        end
        eqs = eqs isa Vector ? eqs : [eqs]
        new(eqs, bcs, domain, ivs, dvs, ps, defaults, connector_type, systems, name)
    end
end

function Base.getproperty(x::PDESystem, sym::Symbol)
    if sym == :indvars
        return getfield(x, :ivs)
        Base.depwarn("`sys.indvars` is deprecated, please use `get_ivs(sys)`", :getproperty,
                     force = true)

    elseif sym == :depvars
        return getfield(x, :dvs)
        Base.depwarn("`sys.depvars` is deprecated, please use `get_dvs(sys)`", :getproperty,
                     force = true)

    else
        return getfield(x, sym)
    end
end

Base.summary(prob::PDESystem) = string(nameof(typeof(prob)))
function Base.show(io::IO, ::MIME"text/plain", sys::PDESystem)
    println(io, summary(sys))
    println(io, "Equations: ", get_eqs(sys))
    println(io, "Boundary Conditions: ", get_bcs(sys))
    println(io, "Domain: ", get_domain(sys))
    println(io, "Dependent Variables: ", get_dvs(sys))
    println(io, "Independent Variables: ", get_ivs(sys))
    println(io, "Parameters: ", get_ps(sys))
    print(io, "Default Parameter Values", get_defaults(sys))
    return nothing
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 HopLink2 <: SwaggerModel
    nextHopId::Any # spec type: Union{ Nothing, String } # spec name: nextHopId
    linkType::Any # spec type: Union{ Nothing, String } # spec name: linkType
    properties::Any # spec type: Union{ Nothing, HopLinkProperties2 } # spec name: properties
    issues::Any # spec type: Union{ Nothing, Vector{ConnectivityIssue2} } # spec name: issues
    context::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: context
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId

    function HopLink2(;nextHopId=nothing, linkType=nothing, properties=nothing, issues=nothing, context=nothing, resourceId=nothing)
        o = new()
        validate_property(HopLink2, Symbol("nextHopId"), nextHopId)
        setfield!(o, Symbol("nextHopId"), nextHopId)
        validate_property(HopLink2, Symbol("linkType"), linkType)
        setfield!(o, Symbol("linkType"), linkType)
        validate_property(HopLink2, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(HopLink2, Symbol("issues"), issues)
        setfield!(o, Symbol("issues"), issues)
        validate_property(HopLink2, Symbol("context"), context)
        setfield!(o, Symbol("context"), context)
        validate_property(HopLink2, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        o
    end
end # type HopLink2

const _property_map_HopLink2 = Dict{Symbol,Symbol}(Symbol("nextHopId")=>Symbol("nextHopId"), Symbol("linkType")=>Symbol("linkType"), Symbol("properties")=>Symbol("properties"), Symbol("issues")=>Symbol("issues"), Symbol("context")=>Symbol("context"), Symbol("resourceId")=>Symbol("resourceId"))
const _property_types_HopLink2 = Dict{Symbol,String}(Symbol("nextHopId")=>"String", Symbol("linkType")=>"String", Symbol("properties")=>"HopLinkProperties2", Symbol("issues")=>"Vector{ConnectivityIssue2}", Symbol("context")=>"Dict{String, String}", Symbol("resourceId")=>"String")
Base.propertynames(::Type{ HopLink2 }) = collect(keys(_property_map_HopLink2))
Swagger.property_type(::Type{ HopLink2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_HopLink2[name]))}
Swagger.field_name(::Type{ HopLink2 }, property_name::Symbol) =  _property_map_HopLink2[property_name]

function check_required(o::HopLink2)
    true
end

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


text:
@testset "NTT 1D" begin
    g = 2
    q = 31
    x = [1:5;]
    @test intt(ntt(x, g, q), g, q) == x

    g = 3
    q = 257
    x = [1:256;]
    @test intt(ntt(x, g, q), g, q) == x
end

@testset "NTT 1D convolution" begin
    using DSP

    g = 9
    q = 271
    # padding inputs with zeros to get non circular convolution
    x = [1:8; zeros(Int64, 7);]
    h = [1:8; zeros(Int64, 7);]

    X = ntt(x, g, q)
    H = ntt(h, g, q)
    Y = X .* H

    y = intt(Y, g, q)

    @test y == DSP.conv(1:8, 1:8)
end

@testset "NTT 2D" begin
    g = 16
    q = 257
    x = reshape([1:16;] , (4,4))
    y = ntt(x, g, q)
    @test intt(y, g, q) == x

    g = 4
    q = 17
    x = reshape([1:16;] , (4,4))
    y = ntt(x, g, q)
    @test intt(y, g, q) == x

    g = 7
    q = 19
    x = reshape([1:9;] , (3,3))
    y = ntt(x, g, q)
    @test intt(y, g, q) == x
end

@testset "NTT 2D convolution" begin
    using DSP

    (g, q, n) = (7, 4733, 7)

    #zero padding to get non circular convolution
    x = reshape([1:16;] , (4,4))
    x_padded = zeros(Int64, 7, 7)
    x_padded[1:4, 1:4] = x

    h = reshape([4 3 2 0; 4 1 0 0; 0 0 1 2; 0 0 2 3] , (4,4))
    h_padded = zeros(Int64, 7, 7)
    h_padded[1:4, 1:4] = h

    X = ntt(x_padded, g, q)
    H = ntt(h_padded, g, q)
    Y = X .* H

    y = intt(Y, g, q)

    @test y == DSP.conv(x, h)
end