text:
Within firearms, chamber pressure is the pressure exerted by a cartridge case's outside walls on the inside of a firearm's chamber when the cartridge is fired. The volume of a cartridge case is equal to the chamber volume minus the volume of the metal making up the case and primer, minus the volume of the extractor groove. In this thread I am hoping to focus on the benefit, perceived benefit, or if it is a … This could help answer a few questions about the overall experiment and to inform future experiments with larger sample sizes. As in the AMP Annealing article, they too identified that most weight variation is in the case head. And as far as I know, setting a measure for 100 grains of BP will get you 72 or so grains of pyrodex, but the two charges will be ballistically similar. To get a better sense of how much it mattered to my readers, I asked a few questions about brass sorting in a survey. I use a syringe with a needle instead of a dropper. There are YouTube videos and write ups on this subject. I use my .45 acp case as a measure for my Ruger old army for plinking with 3F Goex BP. The case and water here weigh 304.3 grains.   The volume is the measured difference in weight between the empty case and filled. Therefore, any weight variations with those cases must be in the case heads. Battery and protective storage case included. QuickLOAD relies on case volume as a critical dimension for its internal ballistics simulations. Personally, I expected there to be less variation in weight at this point – this data does not indicate that, however. In his Measuring the Case article, on his Rifleman’s Journal website, German reviews the tools and techniques required to measure everything from case overall length to case neck concentricity. Isopropyl alcohol (99.9%) was used because it has a number of benefits over distilled water. Thanks for taking the time to read this article! I enjoyed writing it and learned a lot in this process and I hope that you did too. ) to that of distilled water ($W_{\text{H₂O}}$ If you have any feedback, you can email me directly if you don’t prefer to use Reddit or other social media. The cartridge headspaces on the mouth of the case at the L3 datum reference. ). _gaq.push(['_trackPageview']); Some use water, some use powder, some measure case weight, and spend a lot of time sorting brass with the hope of improving accuracy. There are several options for case trimmers from those that are hand operated to a variety of trimmers powered by electric motors. Finally, just 12% of the shooters said they weight-sort their brass. Figure 1 illustrates This liquid volume measurement method can be practically employed to about a 0.01 to 0.02 ml or 0.15 to 0.30 grains of water precision level for firearms cartridge cases. {:.hide-on-med-and-up}. If these estimates are correct, then a 0.50 grain (H₂O) case volume variation translates to a significant +/- 14 ft/s muzzle velocity spread. One version will work with all small-primer-pocket cases, while another works with all large-primer-pocket cartridge types. You may not reproduce original text or photos without advanced permission in writing. © 2017-2021 Ammolytics, All Rights Reserved. What these charts made me realize is that I had set up my resizing die to bump the shoulder back based on a once-fired case. So, if it’s both important and difficult, then there’s value in alternative measurement techniques which are reliable, deterministic indicators of case volume. and .22-250 AI versions. Segregating cases by headstamp is as careful as we need be. The following questions may help you to better understand the complexity: I didn’t have a large quantity of cases to start with, from which the heaviest and lightest could be used for the experiment. While reading this article, I only ask that you keep an open mind. But when you see a cartridge like 7X57mm, 8X57mm, or 6.5X55mm, this is what it means. Those who sort their brass by weight may believe that case weight is directly proportional to case capacity and a shortcut to the tedious process of measuring volume with a liquid. What is the relationship between case weight and case volume? Similar to bullet sorting, this operation seems intuitive when the goal is maximizing consistency by minimizing variation. Reloaders go to great lengths to extend the service life of their brass, especially when using expensive, high-end brands. The closer the scatter plots are to the fit line, the stronger the correlation. Before and after every firing and resizing, I measured the following: After the third firing, I started to wonder about the thickness of the case head, so I began to also measure the following: As you might expect, measuring the weight of the cases was pretty straightforward. After all, don’t we all strive to one day see zero variation shot-to-shot on our chronograph from a 10-shot group? Note: It's a good idea to measure the Maximum Case Capacity, overflow and input this value in the Cartridge Dimensions window to make QuickLOAD's calculations more reliable—this is especially true for subsonic loads. Any quality cast .45-caliber grease-groove bullet I choose to shoot will do. What I'm trying to measure is the volume the case will be at pressure in the chamber, not the volume of a new case. Volumetric measurements are normally represented in milliliters (ml) or US fluid ounces (fl oz). Extra time and care is taken in the creation of our cases, producing smaller lots that meet strict quality standards. Case volume, however, is often represented in grains of water (gr H₂O), which is a unit of weight. Rifle Cartridge Case Capacities Compiled By Chuck Hawks . A full .45acp case will hold 25 grains by volume of Pyrodex P, can't help you with the other calibers Just did an eyeball side by side comparison of some pyrodex p and some swiss 4f and the granule size is similar, the pyrodex looks just a wee bit smaller. })(); To comment on a post, click Post Title or "Comment" Link. After measuring a case, I drained the used alcohol into a separate container which I considered “contaminated”. This is the part of rifle cartridge loading that is normally not fun. Because I kept the cases in the same order, I could not use the wet tumbling method with stainless steel pins that I normally opt for – instead I used the plastic box the brass came in with an ultrasonic cleaner. – 18 min read. With my rifle, they measure in the normal and safe range with this load. I'm not making claims about the quality of my ammunition or my process. They don’t relate to internal volume. 5) Measure length of Headspace before and after resizing 6) Measure case length after resizing 7) Weight Sort Hornady 168 match bullets to find 8 with same weight 8) Take 8 sorted bullets and measure each to select 5 with the same length from base to Ogive 9) Install Federal 210M primers in numbered and fully prepped cases Take the intended projectile you want to hand-load for, and place it inside the case. In this experiment, I followed this procedure several times for 50 rifle cases – I can confirm that it’s a real pain and error-prone. Sectional density of the webbing (the base of the casing). var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); At Hornady,® brass is the foundation for what could be the most memorable shot of your lifetime. Well, in your case, you are measuring the bullet diameter. The following surveys, simulations, and articles relate to this topic and are worth consideration before moving on to the experiment. Preventing and controlling spillover while filling to the exact same level was exceptionally challenging. As should be expected, the overall weights shifted downward from the material that was trimmed. Volume was measured with 99.9% isopropyl alcohol after the first, third, and fourth firing, before resizing. As I mentioned earlier, I intentionally did not trim the brass until after the third firing. I would be very grateful if some one who has access to actual cases could measure the gr. var _gaq = _gaq || []; Because I was interested to see how the case dimensions would change, I intentionally did not trim until after the third firing & resizing. A linear regression may provide more insight. I wanted to get a better sense of which dimensions would change after several firings and how long it would take to fully fire-form the cases. But even then, if the meniscus isn’t identical each time, you can still have a reading that is off by as much as a grain. (function() { German Salazar has written a very thorough guide to measuring nearly all the critical dimensions of cartridge brass. Hornady Lock-N-Load Bullet Comparator Basic Set with 7 Inserts Permanent Link to Product Page. Note: The case number labels were omitted from this chart since the order was lost after the third firing. What I notice in the charts above is that the consistency of case weight barely changes while the case volume shifts much more dramatically. This handy “end-cap” lets you quickly measure multiple new brass cases or deprimed FL-sized cases so you can get an average H20 capacity. I maintained the same order of the cases as they came in their original box and took several measurements on each one before resizing, after resizing, and after firing. This is just the first in a series of brass sorting experiments, which I’ll continue to explore. Start filling by putting the needle through the flash hole into the primer so it will drive trapped air out of the primer. Cartridge Cases The foundation of Hornady ® Ammunition. $$W_{\text{H₂O}} = W_{\text{IPA}} + ((W_{\text{IPA}} \times 0.212) \times 0.999)$$. _gaq.push(['_setAccount', 'UA-6323545-2']); Survey Results: Do you sort your rifle brass by capacity? How many times does a case need to be fired until it’s fully fire-formed, and how does that affect its volume? As I detail my reloading practices and equipment, they will likely differ from your own. On this evidence, overall case weight variation was minimal in all packs of Lapua, Peterson and Norma cases. This allows for a very simple conversion from weight to volume, which can be translated to other units of weight, such as grains. On the low end, the casing could be a single molecule thick. CARTRIDGE CASE VOLUME-WEIGHT We're told that constant cartridge case volume is important to accuracy. I start with a cartridge case fireformed in the subject rifle. Personally, I haven’t seen much evidence to support this belief, so this experiment series will begin by exploring the relationship between case weight and case capacity. Disregarding safety and usefulness, the range is: As I said, neither end of this spectrum is remotely useful, but it demonstrates my point – a correlation between case weight and volume does exist. I performed the following steps to trim the brass. The following capacities are based on the approximate weight of plain water measured in grains and were acquired from various sources. How does cleaning the cases affect weight & volume? The brand of case measured is in parenthesis ( ) as follows: F C = Federal, Hor = Hornady, N = Norma, R-P = Remington, Wby = Weatherby, W-W = Winchester. ga.src = ('https:' == document.location.protocol ? Here are a few of the questions I hope to answer: If previous experience has taught me anything, it’s that we should expect to find even more interesting questions along the way. The straight wall list has the total volume listed first and the volume available after a bullet is seated to the average seating depth second. Some may be tempted to brush aside my conclusions due to these differences if they feel my methods are not up to their standards. Here’s what I’ve learned and can conclude based on the evidence from this first brass sorting experiment. You set the measure according to the BP charge, and the pyrodex will measure less grains, but the charge will be the same (weight of pyrodex is less, but the volume is more). The density of cartridge brass is a constant, for all practical purposes. January 8, 2020 Here are the results. One dimension in particular, case volume, is not only important but also quite annoying to measure, which means people will skip it or take shortcuts. Note: QuickLOAD does not have a small rifle primer case for 6.5mm Creedmoor in its library, so the chamber pressures are listed as dangerous. Should the cases be sorted by weight when new and unfired, or after they’ve been fire-formed, resized, and trimmed? This practice is therefore not recommended. The straight wall cartridge list is small, but its what I could work with. The case volume shows slightly less variation after the fourth firing, but this is more likely due to the margin of error in my measurements. If you are an “advanced reloader” or want to be, you should read German’s article. With powder, you have to settle on a charge weight but, just as importantly, you need to charge each cartridge case with the same amount of powder. Does brass weight correlate to factors other than volume which might affect pressure and velocity? Ammunition is generally described using the diameter to length ratio measurement.For example: 1. It's easy to add or remove water less than one drop at a time. The common method for measuring case volume requires filling a fire-formed case completely with distilled water (without air bubbles) and weighing it. Fill the case completely with water, pour the water into the scale pan and add the case. Measuring brass rifle cases by weight and volume. Where in the case does the weight variation occur? The following formula converts the weight of 99.9% alcohol ($W_{\text{IPA}}$ If you are trying to see if case volume is related to case weight, you measurement variations will need to be significantly less than the variation in your case volume or none of your data will mean anything. If you’re on mobile, try landscape mode. Case order goes down and to the right, with case #1 at the top-left. A majority of cases have a volume within 0.1 gr H₂O of the 42.18 gr mark. Before you start any measurement, ensure the scale is zeroed. When it comes to rifle brass and its role in creating consistent and precise ammunition, some of the dimensions measured matter more than others. The case volumes really tighten up after the third firing. It would be, On the high end, the casing could be completely solid. Brands “A”, “B” and “C” showed clear evidence that at least a portion of the case weight variation was in the targeted annealing zone. To avoid errors due to thermal expansion or equipment sensitivity, I took all measurements from the temperature-stable (68°F) environment of my home office. For example, if you use the weight of .45-70 brass along with the outer dimensions of the … Only recently, after I already completed this work, have I purchased transfer pipettes and micropipettes, which should significantly ease the process for future measurements. Since trimming the length of the case would directly affect case weight and volume, I wanted to avoid that influence until I was confident the brass was fully formed. The added weight of the water is then used to establish the liquid volume and hence the case capacity. You can trim the case and even turn the case necks. Everything didn’t go as planned. Deburr the inside mouth of the fireformed and annealed case and ever so slightly bell the mouth of the case. Gear Review: Bruno Concentricity-Checking Tool, Comparative Specs for Various 6XC Brass Sources, H&H Concentricity Gauge and Bullet Straightener, ABCs of Neck Turning — Salazar Shows How It’s Done, Redding Offers New Slant-Bed Concentricity Checker, Tactical Tech: New XTR Signature Rings from Burris, Comparative Burn-Rate Chart Lists 163 Powders, Beretta Acquires Legendary Holland & Holland Gunmaker, New Mexico Law Would Ban Under-12 Youth Firearms Training, Sunday Gunday: How to Find Great Rifles at Live Gun Auctions. Measuring case volume with 99.9% isopropyl alcohol, an A&D scale, and primer pocket plugs. 's of H2O for the .22-250 Rem. I am having difficulty ascertaining what the real case capacity is of a couple of cases. A rimless-cartridge case which is considerably shorter than the chamber will be violently stretched from the centre toward both ends, frequently resulting in a complete separation, usually just in front of the head, which often wrecks the action and may injure the shooter. A correlation certainly exists between the weight of a brass rifle case and its volume, which can be demonstrated with hypothetical and impractically extreme examples. My goal here at Ammolytics, is to explore one variable at a time, measure the i… You may be surprised to see that number so low, but understand most of these guys are intentional about starting off with good quality brass from the same lot. Can primer pockets be uniformed before weighing, since they reduce weight but not affect volume? When you think about it, this makes good sense. The following charts are more interesting than they are experimentally useful or conclusive, but I thought I’d share them anyway since they still taught me something interesting about my reloading method. Is there an easier way to measure case volume? Weight-sorting brass has become a fairly common practice in the world of competitive rifle marksmen and precision reloaders. It would be. Does case trimming (neck turning, neck reaming) affect the muzzle velocity? This was an incredibly detailed article and in some ways goes far beyond my measuring capabilities. Here are a few of the mistakes I made and problems I ran into: Despite these issues, I’m confident that I collected some solid data and that what I learned from my experiment, and my failures, will be useful for myself and others. Step-by-step, German shows how to measure: Case Length, Case Body length … Kept them closed as much as consistent ammunition matters to you, the has... ( fl oz ) as we need be start any measurement, ensure the scale is.! Through the flash hole into the primer so it will drive trapped air out the. Case body wall thickness variation are not up to their standards does that its... Measured in grains of water to fill with water you want to be less variation in weight between the,. Written a very thorough guide to measuring nearly all the critical dimensions of cartridge brass (... A unit of weight supporting Ammolytics, 21st Century Shooting primer pocket plugs are not up to standards! Reading this article lead me to believe that sorting by weight when new and unfired,,! Relationship between consistent internal volume is the maximum range of weight and volume, it was interesting to see mentioned! Ever so slightly bell the mouth of the webbing ( the base of the casing could be single! A wider screen firing, which explains the how to measure cartridge case volume in average weight variations with those must! Should be expected, the cases affect weight & volume of plain water in! In weight at this point – this data does not indicate that, however the part of rifle cartridge that... Patience and support s article what the real case capacity is of dropper... ’ ll continue to explore that these charts do not imply causation precision rifle Shooters load. Not a competitive shooter, measuring your brass can provide important safety benefits measure... Screen and accuracy to.001 '' the part of rifle casings with identical dimensions. Lost after the third firing as careful as we need be a single molecule thick not overstated... The extractor groove is likewise constant annealed case and ever so slightly bell the mouth of the capacity. Drop into the primer pocket plugs it means that the bullet diameter is 6.5mm in diameter, and place inside! Cast.45-caliber grease-groove bullet I choose to shoot will do volume will be in cubic millimeters, which is constant! Down and to inform future experiments with larger sample sizes plots are to right! Resulting pressure from ignition milliliters ( ml ) or US fluid ounces fl... A majority of cases in chamber pressure or muzzle velocity do you sort your rifle brass resides process... Performed the following surveys, simulations, and fourth firing said they their... Keeping it free of ads and paid-promotions brass until after the third firing consistency by variation. What variation there was did not trim the brass case real case capacity is of a.! And weighing it makes good sense can fit into the case completely with distilled water, around the and..., this operation seems intuitive when the goal is maximizing consistency by minimizing variation the datum. Case order goes down and to the price of Lapua when you think about,! Could help answer a few questions about the quality of my ammunition or my process thanks taking. It mentioned as it related to case volume with how to measure cartridge case volume % ) was used it! Is by measuring the powder charge number of benefits over distilled water one! M not making claims about the quality of my ammunition or my process conclusions...
Lake Chatuge Homes For Sale, Psalm 44 The Message, Ealing Catchment Areas, Mexican Pasta With Ham, Timberwolf Tps35 Parts, Hwasa Twit Lyrics, Big Bamboo Irish Moss Review, Meadowsweet Growing Zone,


text:
# Application examples for the Markovchart R package

#### 2022-04-09

Control charts originate from industrial statistics, but are constantly seeing new areas of application, for example in health care (Thor et al. 2007, suman2018control). This paper is about the Markovchart package, an R implementation of generalised Markov chain-based control charts with health care applications in mind and with a focus on cost-effectiveness. The methods are based on Zempléni et al. (2004), Dobi and Zempléni (2019a), Dobi and Zempléni (2019b).

This vignette highlights the flexibility of the methods through the modelling of two different disease progression and treatment scenarios, namely for monitoring low-density lipoprotein (LDL) and blood sugar (glycated haemoglobin - HbA1c) levels. The vignette focuses on the practical application of the package’s functions from the viewpoint of a user in the field of biostatistics or health care. For the mathematical background and programming details consult the above-cited papers and the functions’ help pages respectively. For cost calculations, the 2021 March 21 EUR-HUF exchange rate was used (1 EUR = 369.05 HUF).

Note that the code used will often limit the parallelisation of the algorithms for the sake of not to over-burden the system of the user. However, given a powerful hardware setup the package can run much faster.

# Low-density lipoprotein analysis

In the following paragraphs we show a real life healthcare example for LDL monitoring. Two approaches will be presented: one with and one without taking the standard deviation into account. The aim is to minimise the healthcare burden generated by patients with high cardiovascular (CV) event risk. The model is built upon the relationship between the low-density lipoprotein level and the risk of CV events, thus the LDL level is the process of interest (Boekholdt et al. 2012).

Parameters were estimated from several sources. The table below gives information about the meaning of the parameter values in the healthcare setting and shows the source of the parameter estimations.
Parameter value Meaning Parameter source
3 mmol/l Target value Set according to the European guideline for patients at risk
0.1 mmol/l Process standard deviation Estimated using real life data from Hungary, namely registry data through Healthware Consulting Ltd.
0.8/3 Expected shift size, 0.8 increase in LDL per year on average Estimated with the help of a health professional
1/120 Expected number of shifts in a day, 3 shifts per year on average Estimated with the help of a health professional
0.027, 1.15 Parameters of the repair size beta distribution Estimated using an international study which included Hungary
0.1, 30 Parameters of the sampling probability logistic function Patient non-compliance in LDL controlling medicine is quite high, and this is represented through the parametrisation of the logistic function
5.01 EUR Sampling cost Estimated using the LDL testing cost and visit cost in Hungary
4.60 EUR Shift-proportional daily out-of-control (OOC) cost Estimated using real world data of cardiovascular event costs from Hungary
9.97 EUR Base repair cost Estimated using the simvastatin therapy costs in Hungary
7.48 EUR Shift-proportional repair cost Estimated using the simvastatin therapy costs in Hungary

It is very difficult to give a good estimate for the type and parameters of a distribution that properly models the non-compliance, thus the results here can at best be regarded as close approximations to a real life situation. This is not a limiting factor, as patients themselves can have vast differences in their behaviour, so evaluation of different scenarios are often required, and will also be presented here. Since high LDL levels rarely produce noticeable symptoms, the sampling probability only depends on the time between samplings (Institute for Quality and Efficiency in Health Care 2017). Thus, the sampling probability was modelled by the logistic function and not by the beta distribution (see the help page of the Markovstat function for more details). It is important to note that the proportional costs increase according to a Taguchi-type (squared) loss function, thus huge, expenses can be generated if the patient’s health is highly out-of-control.

## Optimisation using only the cost expectation

The LDL aggregated dataset is included in the package with all of the above-mentioned parameter estimations. For cost estimation and optimisation purposes we need to use two function: Markovstat for stationary distribution estimation and Markovchart for cost calculation. The code below feeds the proper parameters to these functions. Notice that exponential shift size distribution, random repair and sampling is used and Markovchart is set for optimisation.

data("LDL")

stat_LDL_cost <- Markovstat(
shiftfun = 'exp', h = 50, k = 0.15,
sigma = LDL$standard_deviation, s = LDL$num_exp_daily_shifts,
delta = LDL$expected_shift_size, RanRep = TRUE, alpha = LDL$rep_size_first, beta = LDL$rep_size_second, RanSam = TRUE, q = LDL$samp_prob_first, z = LDL$samp_prob_second, Vd = 50, V = 3) #Defining parallel_opt parallel settings. #parallel_opt can also be left empty to be defined automatically by the function. num_workers <- min(c(detectCores(),2)) parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE) res_LDL_cost <- Markovchart( statdist = stat_LDL_cost, OPTIM=TRUE, p = 1, cs = LDL$sampling_cost,
cf = LDL$base_rep_cost, coparams = c(0,LDL$OOC_cost),
crparams = c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
res_LDL_cost
#> $Results #> G-value Expected cost Total cost sd Cost sd due to process var. #> 1 0.4258544 0.4258544 0.5041891 0.5041891 #> 2nd process moment 3rd process moment 4th process moment #> 1 0.4355586 3.703977 77.98069 #> #>$Subcosts
#>   Sampling cost Repair cost OOC cost
#> 1    0.08998974  0.08023143 0.262028
#>
#> $Parameters #> Time between samplings (h) Critical value (k) #> 1 55.70492 0.1583099 The optimal parameters for the case when the cost standard deviation was not taken into account were 55.7 days and 0.158 mmol/l for the time between samplings and the critical increase in the LDL level from the guideline value respectively. These parameters entailed an average daily cost of 0.43€ and standard deviation of 0.5€. This result is interesting, because according to the optimisation we should use a somewhat higher critical value than the one in the guideline – 0 mmol/l critical increase would be the 3 mmol/l value in the guideline. At the same time patients should be monitored more often than usual as times of the LDL measurements are usually several months or years apart. However, this is a strictly cost effective viewpoint which could be overwritten by a health professional. Nonetheless, the results provide a possible new approach to the therapeutic regime for controlling LDL level. Often, it is good to look at the interaction between the parameters and the resulting average cost, especially in situations where the optimal parameters cannot be used because of real life reasons. parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE) res_LDL_cost_grid <- Markovchart( statdist = stat_LDL_cost, h=seq(50,75,2.5), k=seq(0.05,0.25,0.02), p = 1, cs = LDL$sampling_cost,
cf = LDL$base_rep_cost, coparams = c(0,LDL$OOC_cost),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
plot(res_LDL_cost_grid,
y = 'Expected \ndaily cost',
mid = '#ff9494',
high = '#800000',
xlab = 'Days between samplings',
ylab = 'Critical LDL increase') +
geom_point(aes(x = res_LDL_cost$Parameters[[1]], y = res_LDL_cost$Parameters[[2]]))

The figure shows the expected cost as function of the time between samplings and the critical value.

The heat map shows the average cost as the function of the different parameter values. The dot in the lightest area of Figure $$\ref{expectedcost}$$ corresponds to the optimal cost. Any other point entails a higher average cost. It can be seen that too low or high critical values will both increase the average daily cost. Note that the change in time between samplings entails relatively low change in the critical LDL increase: even if the time between control visits is changed the critical value should stay around 0.12 - 0.18 mmol/l.

## Optimisation using cost expectation and cost standard deviation

The code below is largely the same as before, but notice that p=0.9 instead of p=1, thus prompting the 10% use of the standard deviation during the optimisation process.

stat_LDL_cost <- Markovstat(
shiftfun = 'exp', h = 50, k = 0.15,
sigma = LDL$standard_deviation, s = LDL$num_exp_daily_shifts,
delta = LDL$expected_shift_size, RanRep = TRUE, alpha = LDL$rep_size_first, beta = LDL$rep_size_second, RanSam=TRUE, q = LDL$samp_prob_first, z = LDL$samp_prob_second, Vd = 50, V = 3) parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE) res_LDL_G <- Markovchart( statdist = stat_LDL_cost, OPTIM=TRUE, p = 0.9, cs = LDL$sampling_cost,
cf = LDL$base_rep_cost, coparams = c(0,LDL$OOC_cost),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
res_LDL_G
#> $Results #> G-value Expected cost Total cost sd Cost sd due to process var. #> 1 0.4266033 0.4323714 0.3746903 0.3746903 #> 2nd process moment 3rd process moment 4th process moment #> 1 0.3273378 2.038361 42.26514 #> #>$Subcosts
#>   Sampling cost Repair cost  OOC cost
#> 1    0.07838653  0.08099193 0.2755366
#>
#> $Parameters #> Time between samplings (h) Critical value (k) #> 1 63.95067 0.145161 In this part, the cost standard deviation is also taken into account as p=0.9, thus the weight of the standard deviation in the calculation of $$G$$ is 0.1. The optimal parameters found by our approach were 63.95 days and 0.145 mmol/l for the time between samplings and critical increase in the LDL level respectively. These parameters entailed an average daily cost of 0.43€ and standard deviation of 0.37€. The inclusion of the cost standard deviation into the model has somewhat increased the time between control visits and decreased the critical value. The expected cost increased, while the cost standard deviation has moderately decreased. The figure below shows the previous heat map with non-compliance included in the model. parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE) res_LDL_G_grid <- Markovchart( statdist = stat_LDL_cost, h=seq(50,75,2.5), k=seq(0.05,0.25,0.02), p = 0.9, cs = LDL$sampling_cost,
cf = LDL$base_rep_cost, coparams = c(0,LDL$OOC_cost),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)
plot(res_LDL_G_grid,
y = 'Expected \ndaily cost',
mid = '#ff9494',
high = '#800000',
xlab = 'Days between samplings',
ylab = 'Critical LDL increase',
nbreaks = 14) +
geom_point(aes(x = res_LDL_G$Parameters[[1]], y = res_LDL_G$Parameters[[2]]))

The figure shows $$G$$ value (with p=0.9) as function of the time between samplings and the critical value. It can be seen that the elliptical shape of the heat map has not changed: the change in the time between control visits still does not entail great change in the critical value.

## Sensitivity Analysis

As there were uncertainties about the estimation of several parameters, it is important to assess the effect of different parameter setups. The results for different out-of-control costs are plotted for both approaches. The results can be seen in the figure below. The code runs the optimisation repeatedly to collect the necessary data and then gathers the results into a ggplot-compatible data format.

results <- NULL
statds  <- NULL
for (i in 3:10)
{
for (j in c(0.9,1))
{
res_LDL_optim <- Markovchart(
statdist = stat_LDL_cost,
OPTIM = TRUE,
p = j,
cs = LDL$sampling_cost, cf = LDL$base_rep_cost,
coparams = c(0,i),
crparams= c(LDL$base_rep_cost,LDL$prop_rep_cost),
parallel_opt = parall)

results <- rbind(results,c(j,i,unlist(c(res_LDL_optim$Results[c(2,3)],res_LDL_optim$Parameters))))
statds  <- rbind(statds,res_LDL_optim$Stationary_distribution) } } results <- resultsbck <- as.data.frame(results) colnames(results) <- c("p","co","EC","SDC","h","k") results$p           <-  as.character(results$p) results <- cbind(results[,1:2],melt(results[,3:6])) results$variable    <-  as.character(results$variable) results$variable[results$variable=="h"] <- "Days between samplings" results$variable[results$variable=="k"] <- "Critical LDL increase" results$variable[results$variable=="EC"] <- "Expected cost" results$variable[results$variable=="SDC"] <- "Cost standard deviation" results$variable = factor(results$variable, levels=c("Critical LDL increase","Days between samplings","Expected cost","Cost standard deviation")) results$min <-  NA

results$min[results$variable=="Days between samplings"] <- min(results$value[results$variable=="Days between samplings"])
results$max[results$variable=="Days between samplings"] <- max(results$value[results$variable=="Days between samplings"])
results$min[results$variable=="Critical LDL increase"]  <- min(results$value[results$variable=="Critical LDL increase"])
results$max[results$variable=="Critical LDL increase"]  <- max(results$value[results$variable=="Critical LDL increase"])
results$min[results$variable=="Expected cost"]  <- min(results$value[results$variable=="Expected cost"])
results$max[results$variable=="Expected cost"]  <- max(results$value[results$variable=="Expected cost"])
results$min[results$variable=="Cost standard deviation"]    <- min(results$value[results$variable=="Cost standard deviation"])
results$max[results$variable=="Cost standard deviation"]    <- max(results$value[results$variable=="Cost standard deviation"])

app_LDL_params_plot <- ggplot(results, aes(x=co, y=value, group=p)) +
geom_line(aes(colour=p),size = 1.1) +
facet_wrap(~variable, labeller = label_bquote(rows=.(variable)), scales="free_y", nrow=2) +
scale_colour_manual(expression(italic(p)),values=c("black","#800000")) +
geom_blank(aes(y = min)) +
geom_blank(aes(y = max)) +
xlab("Out-of-control cost") +
ylab("Value") +
theme_bw() +
theme(strip.text.x = element_text(size = 11),
strip.text.y = element_text(size = 11,angle = 0), legend.position="top")

app_LDL_params_plot

The figure shows parameters, average total cost and cost standard deviation as function of the out-of-control cost. The critical value and time between samplings decrease and the average cost and cost standard deviation increase with higher out-of-control costs. Just as on the heat maps, one can observe here, that if the cost standard deviation is taken into account in the optimisation, then the critical value should be lowered and the time between samplings increased. It can be observed that a substantial decrease can be achieved in the cost standard deviation while the cost expectation barely changes.

# Glycated haemoglobin

In this section we will apply the Markovchart package to a randomised real-world data of diabetic patients. The package provides this patient-level pseudonymised and randomised data as the diabetes dataset. The help page of the dataset provides code for data extraction that is useful for control chart applications. This code is skipped here, as it is technical and can be vastly different between applications and datasets. Nonetheless, the reader is invited to study the help page if interested, but the analysis below is based on the already aggregated, useful information gathered from the patient-level data.

## Patient data

The analysis is based on a month-aggregated time series data of diabetic patients from Hungary which was gathered from the period of 2007 September - 2017 September. The data came from two sources: the National Health Insurance Fund of Hungary and the South-Pest Central Hospital. The first source provided information about diagnoses, treatments, health care event and related costs while the latter provided laboratory data regarding blood sugar level. Patients with International Classification of Diseases (ICD) codes (diagnosis) of E10, E11 and E14, and at least one blood sugar measurement were included initially. Only the data of patients with at least one E11 (type II diabetes) diagnosis in the study period was kept. An additional homogenising filter was the requirement of age above 40 at the time of the first diagnosis. Disease progression and therapy effectiveness estimation required at least two blood sugar (HbA1c) measurements with simultaneous therapy data. A total of 4434 patients satisfied all conditions out of which 2151 had at least two HbA1c measurements.

The example study focused on two therapy types: insulin analogues (artificial insulins) and glucagon-like peptides (GLP, promotes insulin secretion). Of course there are more treatment types, the database also lists oral antidiabetics (OAD) and human insulins, but we chose to make the analysis simpler by focusing on GLP and analogue therapies. For the sake of comparison the therapies are grouped in this way: the first group is insulin analogues with possible parallel OAD therapies but human insulin and GLP excluded. The second group is GLP therapies with possible parallel OAD and insulin analogue therapies but human insulin excluded. Essentially we are comparing the effect and cost of insulin analogues with the effect and cost of additional GLP therapies.

The monitored characteristic of the control chart is the blood sugar level, namely the HbA1c level. This is the glycated haemoglobin measured in percent and is slow-changing. The medical aspects of the disease, therapies and blood sugar measurement will not be discussed further as these fall outside of the scope of this example.

The data contains sensitive, patient level information, thus we would only able to show figures, statistics and aggregated data here. However, we pseudonymised and randomised the data in a way that creates fake data but still retains most of the important characteristics and the connections between variables. Namely, random numbers were added (or subtracted) to the date of sampling, the number of sampling per month, costs, HbA1c measurements (both average and standard deviation) and the age of the patient. Furthermore, a subsample was taken from the available patient sample to not even keep the order of patients. The otherwise often very complicated therapeutic regimens were simplified into GLP and analogue categories. This simplification sometimes meant the overwriting of the true therapy used. This is meant to further complicate the identification of patient by e.g., therapeutic pattern. The number of sample elements in the final, randomized dataset can be seen in the table below:

Patient group Total Analogue GLP
Total (all have E11 and are over 40) 800 630 170
E11 ICD 492 272 99

## Parameter estimation

This section provides succinct information about parameter estimation and thus the disease, therapies and the related costs. The aim of the application is to show as many aspects of the Markovchart package to an R user as possible. Due to this, some estimates come from simplified calculations. During industrial or academic applications a group of data scientists and healthcare professionals could provide more accurate estimates, but the work of such a research group is out of the scope of this paper.

Parameters related to disease progression, namely the expected number of shifts per unit time (days) (inverse of shift intensity argument s in the Markovchart function) and the expected shift size delta (delta argument in the Markovchart function) were estimated using the time series data of HbA1c measurements. The delta estimate was calculated from a filtered dataset where we required a HbA1c change larger than $$2\sigma$$ (process standard deviation) and more than 90 but less than 184 days between samplings. The former was necessary to find actual shifts and not just random fluctuations and the latter to try to estimate the size of one shift and not the size of multiple stacked shifts. The expected shift size was delta=1.16. The restriction during estimation of the expected time between samplings was less strict: it was required that the samplings should be at maximum one year (<367 days to account for leap years) from each other. The time between shifts were then gathered from this filtered database and averaged. The shift intensity was calculated by taking the reciprocal of this average and was s=0.0045.

Measurement error (the process standard deviation, sigma argument in the Markovchart function) was estimated using lucky anomalies: the HbA1c level should not be measured more frequently than three months, because the measured values should barely change, and, in relation to this, only 4 measurements are supported per year by the National Health Insurance Fund of Hungary. Nonetheless there were 221 cases in the original data set where HbA1c measurement occurred more than once per month. This number was boosted in the randomised dataset to 692. This data essentially provides information about the measurement error as the actual HbA1c level should change only mildly within such a short time frame. Our estimate with this simple method was sigma = 0.34. It is difficult to compare this result with existing literature due to different methodologies but our estimate is close to the one calculated by Phillipov and Phillips (2001) (even on the randomised data).

Therapy effectiveness was estimated using the definition of effectiveness in the Markov chain-based control charts: the proportion of distance from the target value after treatment compared to the distance before. The target value was set to be 4% HbA1c level, which is the lowest healthy level (Wang 2017). Setting it to a higher value would exclude data, as the target value is the lowest considered. When estimating therapy effectiveness, after the initial filtering of the data, we also restricted HbA1c data to cases where the initial value (after which improvement was seen) was 6%. This was because we wanted to see effect of the therapies only in cases where there is some notable deviation from the healthy state. Improvement was defined as $$>2\sigma$$ (sigma argument in the Markovchart function). Again, values lower than this threshold were considered stagnation (which could also be caused by an effective therapy followed by a degradation, thus making it unreliable) and were not included in the effectiveness estimation. To minimise autocorrelation and the effect of degradation parallel to therapies, only samplings where where more than 90, but less than 184 days elapsed between the two were considered. Parameters of the repair size beta distributions were estimated from the mean and the variance of the ratios of consecutive HbA1c values between samplings. The estimated distribution for analogue therapies was $$beta(2.63,2.05)$$ and for the GLP therapies $$beta(3.85, 3.11)$$ The estimated beta distributions for both therapies can be seen in

therap_plot <-  ggplot(data.frame(x = c(0,1)), aes(x)) +
stat_function(fun = dbeta, aes(colour = ' Analogue', linetype =' Analogue'), size = 1, args = list(shape1 = unlist(ANALOGUE)[1], shape2 = unlist(ANALOGUE)[2])) +
stat_function(fun = dbeta, aes(colour = 'GLP', linetype = 'GLP'), size = 1, args = list(shape1 = unlist(GLP)[1], shape2 = unlist(GLP)[2]), alpha = 0.9) +
scale_colour_manual('Therapy type', values = c('black', '#800000')) +
scale_linetype_manual('Therapy type', values = c('solid', 'dashed')) +
xlab('Therapy effectiveness \n(HbA1c level after therapy divided by the level before)') +
ylab('Beta distribution density') +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, size = 11), legend.justification = c(1, 0),
legend.position = c(0.20, 0.63), legend.title = element_blank(), legend.margin = margin(t = 0, r = 0.1, b = 0, l = 0, unit = 'cm'),
legend.key = element_rect(fill = 'white', colour = 'white'))
therap_plot

Therapy effectiveness comparison using the estimated beta distributions.

It can be seen that the GLP therapy is marginally better based on this dataset than the analogue therapy (judging by the mode). The effectiveness of analogue therapies have a somewhat greater variance. Of course, this method does not take into account the HbA1c level itself, only the change (as a proportion).

Random sampling (i.e., patient non-compliance) was not taken into account as our data did not provide information about this. Nonetheless a sensitivity analysis could be conducted to see the effect of different hypothesised compliance levels, but we have not pursued this idea here.

For daily therapy cost estimation (on top of the initial filtering) we also restricted the data to cases where HbA1c levels were less than or equal to 10%. We observed that cost data as a function of HbA1c level becomes unreliable above this value. This may be due to several factors, one being that patients in highly deteriorated states might have other illnesses (comorbidities) and complications. HbA1c values in between-sampling time periods were imputed using linear interpolation. This was necessary, as therapies were ongoing even between samplings. To be able to estimate cost variances by blood sugar level, the HbA1c level was discretised into 150 categories, as this number was proven to be a good compromise between adequate number of sample elements per category and fineness of the categorisation The relationship between the therapy costs, costs variances and the HbA1c level was then estimated using (non-)linear least squares. This was accomplished with the R function nls. GLP therapy and complication (i.e., OOC) costs in relation of the HbA1c level were estimated and calculated using the default linear cost and cost variance functions, (for more information, see the Markovchart function’s help page). However, the analogue therapy cost and cost variance estimation used a function of the form $c_{r}(v) = c_{rb} + \frac{c_{rs_1}}{v + c_{rs_2}},$ where $$v$$ is the distance from the target value, $$c_{rb}$$ is the base repair cost and $$c_{rs_1}$$, $$c_{rs_2}$$ are the parameters of the scaling repair cost. The resulting functions were the following (everything is given in euro): \begin{align*} c_{analogue}(v) & = 3.3 + \frac{-12.22}{v + 9.37}, \ c_{GLP}(v) = 1.69 + 0.07v, \\ \varsigma_{analogue}(v) & = 7.09 + \frac{-3.32}{v + 1.08}, \ \varsigma_{GLP}(v) = 7 - 0.2v, \\ \end{align*} where $$v$$ is the HbA1c level.

The same restrictions, discretisation and nls function was used when estimating the relationship between the OOC (health care event) cost, cost variance and HbA1c level. Additionally, as there may be a lag between a deteriorated patient state and the resulting health care event, a 6-month cumulative cost was calculated. The resulting functions were the following: \begin{align*} c_{o}(v) & = 0.62 + 0.0062{\cal A}_{h}^2(v), \\ \varsigma_{o}(v) & = 5.12 + 0.13{\cal A}_{h}^2(v), \\ \end{align*} where $$v$$ again is the HbA1c level and $${\cal A}_{h}^2(v)$$ is the expected total squared distance from this starting value on the condition that the sampling interval is $$h$$. Of course, the time between samplings, $$h$$, may not be 6 months in every case, but this was deemed to be a good compromise as there was not enough data to to take $$h$$ into account too during the estimation.

The resulting lines fitted to the data together with standard deviation bands can be seen in the figure below.

cost_plot <- ggplot(data = cost_plots, aes(x = HbA1c, y = Therapy cost)) +
geom_ribbon(aes(x = HbA1c, ymin = low, ymax = high, fill = '±SD'), alpha = 0.4) +
geom_line(aes(x = HbA1c, y = Therapy cost, col = 'Fitted line'), size = 1) +
facet_wrap(Data ~ . ,labeller = label_bquote(rows = .(Data)), ncol = 3) +
scale_color_manual(name = '', values = c('Fitted line' = '#800000')) +
scale_fill_manual(name = '', values = c('±SD' = 'grey70')) +
ylab('Daily cost in EUR') +
theme_bw() +
guides(color = guide_legend(order = 1),
fill = guide_legend(order = 2)) +
theme(legend.position = 'top')

cost_plot

Therapy and medical complication costs (€).

It can be seen that analogue and GLP therapies have similar costs, while compared to these, the OOC costs (especially for higher HbA1c levels) are considerably lower. This is expected as diabetes requires constant treatment while complications may or may not occur. Also, it is understandable to keep patients relatively healthy (thus complication-free) even at a high treatment cost.

There are some more parameters that are worth mentioning: for model fit time between samplings and the critical value was also estimated from the data. h was the mean time between samplings (206.22 days), calculated from the randomised patient data. k was determined according to the medical guidelines (7%, see Ton et al. (2013)). constantr=TRUE and ooc_rep=1 were used, as the therapy is constant and occurs even if there is no alarm (alarm states still play an important role as HbA1c reduction can only be initiated by an alarm).

## Stationary distribution and cost estimation with the Markovchart function

Before we can feed the estimated parameters to the Markovchart function, we have to define the above-mentioned custom cost and cost variance functions in R. Afterwards we can run the Markovchart function for both the analogue and the GLP therapy group.

crfun_ANALOGUE <- function(mudist, crparams){
mudist <- mudist
crb    <- crparams[1]
crs    <- crparams[2]
crs2   <- crparams[3]
cr     <- crb + crs / (mudist + crs2)
return(cr)}
vcrfun_ANALOGUE <- function(mudist, vcrparams){
mudist <- mudist
vcrb   <- vcrparams[1]
vcrs   <- vcrparams[2]
vcrs2  <- vcrparams[3]
vcr    <- vcrb + vcrs / (mudist + vcrs2)
return(vcr)}

stat_ANALOGUE <- Markovstat(
shiftfun = 'exp', h = 206.22, k = 3,
sigma = sigma_param, s = s_param,
delta = delta_param, RanRep = TRUE,
alpha = as.numeric(ANALOGUE[1]),
beta = as.numeric(ANALOGUE[2]),
Vd = 100, V = 18)
res_ANALOGUE <- Markovchart(
statdist = stat_ANALOGUE, p = 1,
constantr = TRUE, ooc_rep = 1,
cs = sampling_cost,
coparams = summary(mod.COST)$coef[ , 1], crfun = crfun_ANALOGUE, crparams = summary(mod.ANALOGUE)$coef[ , 1],
vcoparams = summary(mod_var.COST)$coef[ , 1], vcrfun = vcrfun_ANALOGUE, vcrparams = summary(mod_var.ANALOGUE)$coef[ , 1],
parallel_opt = parall)

stat_GLP <- Markovstat(
shiftfun = 'exp', h = 206.22, k = 3,
sigma = sigma_param, s = s_param,
delta = delta_param, RanRep = TRUE,
alpha = as.numeric(GLP[1]),
beta = as.numeric(GLP[2]),
Vd = 100, V = 18)
res_GLP <- Markovchart(
statdist = stat_GLP, p = 1,
constantr = TRUE, ooc_rep = 1,
cs = sampling_cost,
coparams = summary(mod.COST)$coef[ , 1], crparams = summary(mod.GLP)$coef[ , 1],
vcoparams = summary(mod_var.COST)$coef[ , 1], vcrparams = summary(mod_var.GLP)$coef[ , 1],
parallel_opt = parall)

Before inspecting the cost elements of the results, it is useful to check the stationary distribution of the Markov chains and compare them to the empirical HbA1c data. We can do this in a visually appealing fashion by creating a distance distribution from the stationary distribution (by not taking into account the type of the state, only the distance). The comparison can be seen in in the figure below.

int             <- seq(4, 22, by = (18 / (100 - 1))) - (18 / (100 - 1)) / 2
int[1]          <- 4

distance_dist_A <- res_ANALOGUE[[4]][c(1, (100 + 2):(100 * 2))] +
res_ANALOGUE[[4]][2:(100 + 1)]
theo.ANALOGUE   <- cbind('Analogue', as.data.frame(cbind(int,distance_dist_A)))
colnames(theo.ANALOGUE) <- c('Therapy', 'HbA1c', 'Probability')

distance_dist_G    <- res_GLP[[4]][c(1,(100 + 2):(100 * 2))] +
res_GLP[[4]][2:(100 + 1)]
theo.GLP                 <- cbind('GLP', as.data.frame(cbind(int,distance_dist_G)))
colnames(theo.GLP) <- c('Therapy', 'HbA1c', 'Probability')

hba1c_theo  <- rbind(theo.ANALOGUE, theo.GLP)
hba1c_empir$Therapy[hba1c_empir$Therapy=='ANALOGUE'] <- 'Analogue'

hba1c_empir$Density <- NA hba1c_empir$Density[hba1c_empir$Therapy=='Analogue'] <- hba1c_empir$Probability[hba1c_empir$Therapy=='Analogue']/ ((max(hba1c_empir$HbA1c[hba1c_empir$Therapy=='Analogue'])-min(hba1c_empir$HbA1c[hba1c_empir$Therapy=='Analogue']))/(100-1)) hba1c_empir$Density[hba1c_empir$Therapy=='GLP'] <- hba1c_empir$Probability[hba1c_empir$Therapy=='GLP']/ ((max(hba1c_empir$HbA1c[hba1c_empir$Therapy=='GLP'])-min(hba1c_empir$HbA1c[hba1c_empir$Therapy=='GLP']))/(100-1)) hba1c_theo$Density  <- hba1c_theo$Probability/(18/(100-1)) statd_therapies <- ggplot() + geom_bar(data = hba1c_empir, stat = 'identity', width = 0.01, aes(x = HbA1c, y = Density, fill = 'Empirical'), colour = 'black', alpha = .5) + geom_line(data = hba1c_theo, aes(x = HbA1c, y = Density, col = 'Markovchart'), size = 1.2) + facet_wrap(Therapy ~ . ,labeller = label_bquote(rows = .(Therapy)), ncol = 2) + scale_x_continuous(breaks = seq(4, 22, by = 2)) + xlab('HbA1c') + ylab('Density') + theme_bw() + scale_colour_manual(values = c('#800000')) + scale_fill_manual(values = c('black')) + theme(plot.title = element_text(hjust = 0.5, size = 11), legend.justification = c(1,0), legend.position = c(0.99,0.60), legend.title = element_blank(), legend.margin = margin(t = 0, r = 0.1, b = 0, l = 0, unit = 'cm'), legend.key = element_rect(fill = 'white', colour = 'white')) statd_therapies The figure shows empirical and theoretical stationary distributions. We can assess that the stationary distribution fits the data at an acceptable level. This is actually a very beneficial result, as the main focus of the function is cost estimation which is based on the stationary distribution. Now that we have assured that there are no anomalies in the stationary distribution estimation, let us check the expected costs and cost standard deviations. The costs statistics are shown below.  ### Empirical costs for comparison # ANALOGUE mean(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR)/30 +
mean(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) & !grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU) + sampling_cost/mean(deltats[,2]) #> [1] 3.15981 mean(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR)/30
#> [1] 2.323821
mean(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) & !grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU) #> [1] 0.7982125 sampling_cost/mean(deltats[,2]) #> [1] 0.03777634 sd(RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) &
!grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR/30 +
RANDOMISED_DATA[grepl("ANALOGUE", RANDOMISED_DATA$THERAPY) & !grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU + sampling_cost/mean(deltats[,2])) #> [1] 3.872857 res_ANALOGUE #>$Results
#>    G-value Expected cost Total cost sd Cost sd due to process var.
#> 1 3.167946      3.167946      3.752907                   0.2614257
#>   2nd process moment 3rd process moment 4th process moment
#> 1           10.10423           32.46643           105.1605
#>
#> $Subcosts #> Sampling cost Repair cost OOC cost #> 1 0.03777651 2.388222 0.7419474 #> #>$Parameters
#>   Time between samplings (h) Critical value (k)
#> 1                     206.22                  3

# GLP
mean(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR)/30 +
mean(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU) +
sampling_cost/mean(deltats[,2])
#> [1] 2.35354

mean(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR)/30
#> [1] 1.88242
mean(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU)
#> [1] 0.4333437
sampling_cost/mean(deltats[,2])
#> [1] 0.03777634

sd(RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$THERAPY_COST_EUR/30 +
RANDOMISED_DATA[grepl("GLP", RANDOMISED_DATA$THERAPY),]$COST_CUMU +
sampling_cost/mean(deltats[,2]))
#> [1] 3.493658

res_GLP
#> $Results #> G-value Expected cost Total cost sd Cost sd due to process var. #> 1 2.77835 2.77835 3.708439 0.2987063 #> 2nd process moment 3rd process moment 4th process moment #> 1 7.808456 22.24445 64.38612 #> #>$Subcosts
#>   Sampling cost Repair cost  OOC cost
#> 1    0.03777651    2.004543 0.7360308
#>
#> $Parameters #> Time between samplings (h) Critical value (k) #> 1 206.22 3 One can see that the total expected cost per day is 3.17€ for the analogue and 2.78€ for the GLP group, with considerable standard deviation (>3.7€) in both cases. (Note that p=1, thus the $$G$$-value is simply the expected cost.) The main source of the costs is the therapy cost, as we have seen during the parameter estimation. We can also compare the results to the empirical data for additional goodness of fit evaluation. It can be seen that the estimates of the Markovchart function are quite close to the ones estimated directly from the empirical data. Naturally there is room for improvement in some cases. For example the difference in the event costs is quite clear between therapies from the empirical data, while the Markovchart estimates do not reflect this. This is mainly an assumption problem: it is assumed that event (i.e., OOC) costs solely depend on the states of the chain (i.e., HbA1c level), thus due to the similar stationary distributions there will not be much difference between the event costs either. This of course can be solved by separate event cost estimates in the patient groups but since the absolute difference is small it is not important in this case. The sampling costs are exactly the same, as the empirical mean time between samplings (estimated from the whole patient population) was used in all cases. Overall, we can say that both the stationary distributions and the costs emulates the empirical data at an acceptable level. This is useful, because it means that we have a good baseline model which allows scenario exploration (i.e., changing various parameters to check their effects on the costs). We will now use the vectorisation feature of the Markovchart function and explore how the expected daily cost changes with different days between samplings and critical HbA1c levels. The code is almost the same as before, only the value of h and k is changed to numeric vectors. parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE) resmtx_ANALOGUE <- Markovchart( h = seq(90, 365, by = (365 - 90) / 19), k = seq(0.5, 6, by = (6 - 0.5) / 19), statdist = stat_ANALOGUE, p = 1, constantr = TRUE, ooc_rep = 1, cs = sampling_cost, coparams = summary(mod.COST)$coef[ , 1],
crfun = crfun_ANALOGUE,
crparams = summary(mod.ANALOGUE)$coef[ , 1], vcoparams = summary(mod_var.COST)$coef[ , 1],
vcrfun = vcrfun_ANALOGUE,
vcrparams = summary(mod_var.ANALOGUE)$coef[ , 1], parallel_opt = parall ) parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE) resmtx_GLP <- Markovchart( h = seq(90, 365, by = (365 - 90) / 19), k = seq(0.5, 6, by = (6 - 0.5) / 19), statdist = stat_GLP, p = 1, constantr = TRUE, ooc_rep = 1, cs = sampling_cost, coparams = summary(mod.COST)$coef[ , 1],
crparams = summary(mod.GLP)$coef[ , 1], vcoparams = summary(mod_var.COST)$coef[ , 1],
vcrparams = summary(mod_var.GLP)\$coef[ , 1],
parallel_opt = parall
)

We can use the plot method to quickly assess the results. The contour plot for the analogue therapy patient group is presented in the figure below.

resmtx_ANALOGUE[ , 2] <- resmtx_ANALOGUE[ , 2] + 4
plot(x = resmtx_ANALOGUE, y = 'Expected \ndaily cost',
mid = '#ff9494', high = '#800000',
xlab = 'Days between samplings',
ylab = 'Critical HbA1c level')

Contour plot of expected costs (€) related to analogue therapy.

It can be seen that more frequent samplings and lower critical values entail less daily costs on average. The optimal parameters fall outside of the plotted area. This is due to the fact that lower parameter values are not viable. Namely, HbA1c measurements within 90 days of each other become redundant due to the slow variation and HbA1c values lower than 4% are actually indicating too low blood sugar. Nonetheless the results provide useful information about the relationship between the parameters, for example we can see that the expected daily costs is less sensitive on the differences between greater parameter values.

Since we have two different therapies it may be beneficial to compare them. The figure below conveys the same information as the one above, but for the GLP therapy.

resmtx_GLP[ , 2] <- resmtx_GLP[ , 2] + 4
plot(x = resmtx_GLP, y = 'Expected \ndaily cost',
mid = '#ff9494', high = '#800000',
xlab = 'Days between samplings',
ylab = 'Critical HbA1c level')

It can be seen that the relationships between the variables are largely the same, however GLP is somewhat cheaper overall.

# References

Boekholdt, S Matthijs, Benoit J Arsenault, Samia Mora, Terje R Pedersen, John C LaRosa, Paul J Nestel, R John Simes, et al. 2012. “Association of LDL Cholesterol, Non-HDL Cholesterol, and Apolipoprotein b Levels with Risk of Cardiovascular Events Among Patients Treated with Statins: A Meta-Analysis.” Journal of the American Medical Association 307 (12): 1302–9.
Dobi, Balázs, and András Zempléni. 2019a. “Markov Chain-Based Cost-Optimal Control Charts for Health Care Data.” Quality and Reliability Engineering International 35 (5): 1379–95.
———. 2019b. “Markov Chain-Based Cost-Optimal Control Charts with Different Shift Size Distributions.” Annales Universitatis Scientiarum Budapestinensis de Rolando Eötvös Nominatae, Sectio Computatorica 49: 129–46.
Institute for Quality and Efficiency in Health Care. 2017. High Cholesterol: Overview. InformedHealth.org. https://www.ncbi.nlm.nih.gov/books/NBK279318/.
Phillipov, George, and Patrick J. Phillips. 2001. “Components of Total Measurement Error for Hemoglobin A1c Determination.” Clinical Chemistry 47 (105): 1851–53.
Thor, Johan, Jonas Lundberg, Jakob Ask, Jesper Olsson, Cheryl Carli, Karin Pukk Härenstam, and Mats Brommels. 2007. “Application of Statistical Process Control in Healthcare Improvement: Systematic Review.” BMJ Quality & Safety 16 (5): 387–99.
Ton, Van‐Khue, Seth S. Martin, Roger S. Blumenthal, and Michael J. Blaha. 2013. “Comparing the New European Cardiovascular Disease Prevention Guideline with Prior American Heart Association Guidelines: An Editorial Review.” Journal of Statistical Software 36 (5): E1–6.
Wang, Ping. 2017. “What Clinical Laboratorians Should Do in Response to Extremely Low Hemoglobin A1c Results.” Laboratory Medicine 48 (1): 89–92.
Zempléni, András, Miklós Véber, Belmiro Duarte, and Pedro Saraiva. 2004. “Control Charts: A Cost‐optimization Approach for Processes with Random Shifts.” Applied Stochastic Models in Business and Industry 20 (3): 185–200.

format: plain


text:
Microwave oven spot welder

Microwave oven transformer spot welder. Later I have replaced the welding pulse activation switch with a foot switch.

I made a spot welder out of a microwave oven transformer and a starter cable. Those are quite popular actually. They use a solid state relay to switch the transformer on and off. I use a Fotek SSR-40DA DC-controlled 40 A AC-zero-crossing-switching solid state relay. It was only 4 € and free shipping on eBay. First I needed to extract the transformer from a microwave oven and to replace the secondary coil with a starter cable.

Microwave oven (was about 1050 W consumed power) gutted. In addition to the transformer (left), some other goodies were found from the inside that I saved for later use.

The secondary coil sawed off. There was also another secondary coil (red) that I removed later. After sawing from one side, the secondary coil could be hammered out using a chisel.

Some initial trials using a 700 W rated 12 V automotive starter cable. Didn't produce satisfactory welds. The cable was the most expensive part in the project, about 20 €. There wasn't that much copper in the cable originally. Most of the thickness came from the rubbery insulation. Perhaps I should have been more patient to obtain the copper from somewhere else, as I ended up stripping the cable anyhow to double the conductive copper area. Minimization of resistive losses in the secondary is key.

After making my own cable for the secondary, with twice as much copper, I was still able to fit in a few more rounds. This increased the voltage and reduced resistive losses further.

At this point I got a nice 4.6 V AC out of the secondary.

One end of the secondary cable is routed "under the floor", where it will be grounded to the body of the welder, and then goes to the immobile jaw. The upper jaw is mobile. To reduce wear of the secondary cable here, it is shaped like a spiral. A transparent plastic block electrically isolates the upper jaw from the body.

In the bottom, the cable is crimped inside a piece of copper pipe, hammered flat.

Another view of the jaws. I wish I'd found perfectly fitting copper pipes for the sockets. Now they are a bit loose. The copper cable goes pretty much all the way to the end of the jaws to minimize resistance.

A sturdy 230 V AC fan will take care of cooling. It pushes air out through the front grid at great speed, so that sparks will have less chance of flying into the welder through the grid.

I wanted to automate the welding timing so I made this simple adjustable-duration pulse generator circuit out of some components that I had around:

Welding pulse generator circuit.

 Spot welder timer circuit CD4093BC pinout in the same orientation (C) Fairchild Semiconductor

Parts list:

• Resistor 1 kohm (R4)
• Resistor 9 kohm (R1)
• Resistor 40 kohm (R2)
• Multi-turn or log taper pot 10 kohm (VR1)
• 3x Capacitor 10 V 10 uF (you can put the extra capacitor in parallel with C2 to make the pulse longer if needed) (C1, C2)
• CD4093BC Quad two-input NAND Schmitt trigger logic gate

The variable resistor VR1 controls the timing. I used a helical potentiometer but a log pot would probably work just as well. The higher resolution is needed for long welding pulse times for which the whisker of VR1 is close to R2. R1 and R2 adjust the range of the voltage divider. Capacitor C2 is also charged through R1. The voltage divider range should be such that at one extreme, continuous welding is just barely enabled (by keeping the button pressed), and the other extreme just barely gives the shortest welding time that may be needed. If the obtainable shortest welding times are unnecessarily long, the capacitor C2 can be made larger. This will make it easier to get also very long pulse times. Doubling the capacitance should double the timing.

A rough approximation on how the voltage across C2 rises after the button is depressed at time $t = 0$ is given by:

$V_{\text{C2}}= \frac{5\text{V}\cdot t\cdot(R_{\text{R2}} + R_{\text{VR1}})}{C_{\text{C2}}\cdot R_{\text{R1}}\cdot(R_{\text{R2}} + R_{\text{VR1}}) + t\cdot(R_{\text{R1}} + R_{\text{R2}} + R_{\text{VR1}})}$

where $R$ is resistance and $C$ is capacitance across a component. $R_{\text{VR1}}$ is the nominal resistance of the variable resistor, not the current set value. The first half of the variable resistor VR1 and the second half of the variable resistor VR1 + resistor R2 act as a voltage divider that feeds a divided $V_{\text{C2}}$ to a Schmitt trigger which will trigger at nominally 3.3 V, ending the pulse.

Calculated pulse duration vs. potentiometer position. Drain by the IC, and filtering by C1 is not taken into account. In practice the longest duration achievable is about 2 s. Turning the pot further makes the pulse last indefinitely, or until the button is released. The shortest time is in reality also shorter than the 0.27 s shown. The function plotted is y = 2970/(x + 1060).

A logarithmic pot would give almost linear control of the pulse duration.

An initial version of the circuit took about a second to recover before it was able to produce another pulse after releasing the button. Resistor R4 was added to remedy this. It quickly drains the capacitors empty when the button disconnects the circuit from the power supply.

The circuit is fail-safe in the sense that it only gets power through the button. If the button is released, the output soon goes down, even if the IC or anything else in the circuit malfunctions for some reason. I use this circuit to control (the LED + series resistor of) an optically isolating solid state relay directly.

Noise can be a problem at longer pulse durations, where a sudden spike may end the welding pulse prematurely.

The timer is not synchronized to AC, so timing will also vary if the solid state relay is of the zero-crossing switching type.

Cabling added. As an afterthought, I might have used more grounding cable to ground all enclosure surfaces...

Solid state relay (right) and a 5 V power supply added

Enclosure almost fully assembled, only the top is missing. I put in rather hefty fuses into the 230 V AC input connector. The connector at the bottom right corner is for the foot switch.

Showing what a welding pulse does to a bicycle spoke.

Action! Just for the sparks.

Welding results are mediocre at best. But strong enough for some purposes. One problem is that when the workpiece starts to melt, the current increases dramatically, causing a small "explosion". I am not sure how this could be improved upon to give a steady current.

1. Could you somehow measure the current, and using that reading, feed it into a micro using the PID algorithm (http://en.wikipedia.org/wiki/PID_controller) to control the SSD. I’m using a thermocouple in a toaster oven with an atmega running a PID loop to (roughly) follow a reflow curve for surface mount soldering. If you want code, you can’t have it as it doesn’t currently exist; I’m in the process of writing it and developing the project! Plenty of others have done the same. I’d recommend Beauregard’s post on his implementation of the PID algorithm for Arduino: http://brettbeauregard.com/blog/2011/04/improving-the-beginners-pid-introduction/

Very cool, thanks for posting.

Comment by Will — 2013-06-17 @ 22:52

2. Will: I’ve been thinking the same. Perhaps the current could be measured from acoustic emissions, or inductively from the secondary coil. The SSD is zero-crossing switching so the control will be rather coarse. Perhaps automatic pulse stopping would be enough.

Comment by Olli Niemitalo — 2013-06-18 @ 11:05

3. “One problem is that when the workpiece starts to melt, the current increases dramatically, causing a small “explosion”. I am not sure how this could be improved upon to give a steady current.”

I’m no expert but I’d say that’s what it should do-
For welding a round rod as in the photo your timing would be critical so as not to melt it, welding two flat sheets would be a bit easier. But you still may need some means to get the right pressure on the tips. Otherwise looks good.
You do see portable commercial spot welders with no timing circuit at all and you just release the trigger at the right time to get a good weld, the timing circuit would be handy if you had a lot of similar welds to make. The commercial units are however set up so as to get some serious pressure on the tips.

Thanks for the posting, lots for us would/wanna bees to learn from it.

Comment by Old goat — 2013-12-13 @ 06:51

4. Thanks Oldie, that makes sense. For welding round pieces, if there is enough constant pressure, then the molten metal should flow out of the gap and not heat up too much, given that you will end the pulse in time. Probably more pressure is needed when welding plates together, as it will be more difficult to get a good contact with those. But maybe I could add a very small toroidal transformer around the fat secondary cable to sense the current that goes through the workpiece. Would be interesting to graph that at least, to see if it perhaps tells something about what’s happening.

Comment by Olli Niemitalo — 2013-12-14 @ 13:46

5. To measure current into a micro controller like arduino, take a torroid ring core that will fit over an output wire, cut a small section out of the torriod and place a Hall effect sensor in there, measure the magnetic force produced if you know what I mean

Comment by Anonymous — 2014-08-29 @ 17:33

6. Like this

http://www.machsupport.com/forum/index.php/topic,10339.60.html

Comment by Anonymous — 2014-08-29 @ 17:35

7. One thing that would prevent “explosions” is a purge gas like CO2. It would be a better weld too.

Comment by Mike — 2014-10-01 @ 23:19

8. Thanks Mike, that would be a nice addition.

Comment by Olli Niemitalo — 2014-10-02 @ 09:04

9. The sparks are caused by resistance between the copper tip and the metal rapidly heating that area causing expulsions.
Melted copper tips and flying molten metal is not due to too many amps.
Get heavy duty copper for tips and clamp down securely with 70 – 120 kg pressure before welding.
Commercial auto spot welders can put 50,000 amps through their tip and you don’t even see the area glowing.
They use mega amps, short time periods and clamping pressure. The spot of melting in a spot weld should be between the two plates or whatever you join, not on the surface between copper tip and plate.

Anyhoo my point is you don’t have too much current.

Regards doug..

Comment by Doug — 2015-01-07 @ 19:22

10. Just another comment for anyone looking at making a spot welder.
You want the most resistance in your circuit to be at the point between the two plates. That way the heat builds up their first and the molten nugget forms on in inside. Everything good but too little time will give an improperly formed small nugget which is weak. Just the right time will give an ideal sized strong nugget. More time after that does not increase the size of the nugget but starts to be detrimental weakening spot weld.

Size of spot weld tip diameter in contact should be plate thickness + other plate thickness + 3.4mm

So for 2 x 1mm plate a tip that gives a 5.4mm spot weld is ideal. A shallow domed curve is ideal as it assists with pressure initially and broadens quickly, this gives a larger diameter and less pressure towards the end of the spot weld.

Regards Doug.

Comment by Doug — 2015-01-09 @ 14:33

11. Nice project and quiet dangerously,……

Comment by adrian — 2015-03-07 @ 12:30

12. Hi folks, having recently acquired a commercial, “hand-held” (bloody heavy) spot welder WITHOUT any timers or any electronic additions I can say the following:
1) Electronic timers are NOT a requirement for a few incidental welds -just a bit of experience will do.
2) Clamping pressure is vital -you need a LOT of it. (that plus decently clean surfaces; rust / oxidation tends to be a problem due to the low voltage.)

I fully endorse Dough’s comment; he knows well what it’s about.

Comment by Antin — 2015-03-23 @ 18:49

13. You guys are so right, I need more pressure. I just tried a large-ish spot welder and had a hard time on applying enough pressure on the pedal to operate it. It made good welds without a timer.

Comment by Olli Niemitalo — 2015-03-24 @ 08:38

14. What output voltage on pin 4 its if input voltage its 5V. In my case its 3,2V i not enough for my SSR who need 5V input.

Comment by Boris — 2017-10-31 @ 02:54

15. Pulse width modulation of the input AC might control the current.

Comment by Paul Clark — 2018-07-11 @ 23:04

16. […] I’m not going to go into exhaustive detail on how I fabricated this spot welder. A lot of it is just frosting on top of other peoples’ tutorials. However, I am eager to share what I feel are some useful additions you might consider including in your own version. This particular machine was based off of designs by the King of Random, Albert van Dalen, and Olli Niemitalo. […]

Pingback by Scrap Microwave Spot Welder • HAR.MS — 2019-01-01 @ 21:34


text:
Share

# Find the Unit Vector in the Direction of Vector Bar(Pq) , Where P and Q Are the Points (1, 2, 3) and (4, 5, 6), Respectively. - CBSE (Science) Class 12 - Mathematics

#### Question

Find the unit vector in the direction of vector bar(PQ) , where P and Q are the points (1, 2, 3) and (4, 5, 6), respectively.

#### Solution

The given points are P (1, 2, 3) and Q (4, 5, 6).

Is there an error in this question or solution?

#### Video TutorialsVIEW ALL [1]

Solution Find the Unit Vector in the Direction of Vector Bar(Pq) , Where P and Q Are the Points (1, 2, 3) and (4, 5, 6), Respectively. Concept: Properties of Vector Addition.
S


text:
Primitive of Function of Arccosine

Theorem

$\ds \int \map F {\arccos \frac x a} \rd x = -a \int \map F u \sin u \rd u$

where $u = \arccos \dfrac x a$.

Proof

First note that:

 $\ds u$ $=$ $\ds \arccos \frac x a$ $\ds \leadsto \ \$ $\ds x$ $=$ $\ds a \cos u$ Definition of Arccosine

Then:

 $\ds u$ $=$ $\ds \arccos \frac x a$ $\ds \leadsto \ \$ $\ds \frac {\d u} {\d x}$ $=$ $\ds \frac {-1} {\sqrt {a^2 - x^2} }$ Derivative of Arccosine Function: Corollary $\ds \leadsto \ \$ $\ds \int \map F {\arccos \frac x a} \rd x$ $=$ $\ds \int \map F u \paren {-\sqrt {a^2 - x^2} } \rd u$ Primitive of Composite Function $\ds$ $=$ $\ds \int \map F u \paren {-\sqrt {a^2 - a^2 \cos^2 u} } \rd u$ Definition of $x$ $\ds$ $=$ $\ds \int \map F u a \paren {-\sqrt {1 - \cos^2 u} } \rd u$ $\ds$ $=$ $\ds \int \map F u a \paren {-\sin u} \rd u$ Sum of Squares of Sine and Cosine $\ds$ $=$ $\ds -a \int\map F u \sin u \rd u$ Primitive of Constant Multiple of Function

$\blacksquare$


text:
# Homework Help: Derivative of (x^2+5x-1)/(x^2)

1. Jun 15, 2010

### theBTMANIAC

1. The problem statement, all variables and given/known data Find the first derivative of [...] (below)

2. Relevant equations

((x^2)+5x-1)/(x^2)

3. The attempt at a solution

(x^2+5x-1)/(x^2)

$$y' = \frac{(x^2)(2x+5-0)-((x^2)+5x+-1)(2x)}{x^4} \\ y' = \frac{2x^3+5x^2-2x^3-10x^2+2x}{x^4} \\ y' = \frac{(-5x^2)+2x}{x^4} \\ y' = \frac{-5x+2}{x^3}$$

Is this correct? It feels "off."

Last edited: Jun 15, 2010
2. Jun 15, 2010

### matt_crouch

thats what i got.... =]

3. Jun 15, 2010

### theBTMANIAC

Cool. Thank you. :)

4. Jun 15, 2010

### HallsofIvy

However, the simplest way to do that is write
$$\frac{x^2+ 5x- 1}{x^2}= \frac{x^2}{x^2}+ \frac{5x}{x^2}- \frac{1}{x^2}= 1+ 5x^{-1}- x^{-2}$$

The derivative of that is, of course, $-5x^{-2}+ 2x^{-3}$.

Can you see that that is equal to your result?

5. Jun 15, 2010

### Anonymous217

You can also apply the product rule if HallsofIvy's method doesn't work out (for example, if the denominator is more complicated).


text:
LLVM  8.0.0svn
llvm::MCDwarfFile Struct Reference

Instances of this class represent the name of the dwarf .file directive and its associated dwarf file number in the MC file, and MCDwarfFile's are created and uniqued by the MCContext class where the file number for each is its index into the vector of DwarfFiles (note index 0 is not used and not a valid dwarf file number). More...

#include "llvm/MC/MCDwarf.h"

Collaboration diagram for llvm::MCDwarfFile:
[legend]

## Public Attributes

std::string Name

unsigned DirIndex

MD5::MD5ResultChecksum = nullptr
The MD5 checksum, if there is one. More...

Optional< StringRefSource
The source code of the file. More...

## Detailed Description

Instances of this class represent the name of the dwarf .file directive and its associated dwarf file number in the MC file, and MCDwarfFile's are created and uniqued by the MCContext class where the file number for each is its index into the vector of DwarfFiles (note index 0 is not used and not a valid dwarf file number).

Definition at line 50 of file MCDwarf.h.

## ◆ Checksum

 MD5::MD5Result* llvm::MCDwarfFile::Checksum = nullptr

The MD5 checksum, if there is one.

Non-owning pointer to data allocated in MCContext.

Definition at line 59 of file MCDwarf.h.

## ◆ DirIndex

 unsigned llvm::MCDwarfFile::DirIndex

Definition at line 55 of file MCDwarf.h.

## ◆ Name

 std::string llvm::MCDwarfFile::Name

## ◆ Source

 Optional llvm::MCDwarfFile::Source

The source code of the file.

Non-owning reference to data allocated in MCContext.

Definition at line 63 of file MCDwarf.h.

The documentation for this struct was generated from the following file:


text:
# Customizing the height of a figure

How do I to increase the height of a figure? and keeping its width like it's mention in the below code.

\begin{figure}[H]
\hfill{}\includegraphics[width=\linewidth]{images/cu_global_initial}\hfill{}
\caption{Diagramme de cas d'utilisation global initial}
\label{fg:diag_cu_glbl}
\end{figure}

• Additionally giving height and keepaspectration should do the trick, right? – Johannes_B Feb 14 '14 at 11:23
• If you want to center your image use \centering. The \hfill construction you show is a bit strange, the first \hfill will be discarded as it is at the start of the line, the second \hfill would make it flush left, except that the image is full width so there is no alignment possible. – David Carlisle Feb 14 '14 at 11:33
• Do you want to do "as if" the figure had a width of \linewidth and a given height, or do you want to zoom the figure up those two dimensions? – Bernard Feb 14 '14 at 12:10

Below is a code for both cases:

        \documentclass[12pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}

\usepackage{graphicx}

\begin{document}
Voici une figure :

\begin{figure}[!h]
\centering\fbox{\raisebox{0pt}[6cm]{\includegraphics{airepara}}}%[width=\linewidth]
\caption{Diagramme de cas d'utilisation global initial
\label{fg-aire}}
\end{figure}

Comme on le voit sur la figure \ref{fg-aire}, blablabla…

S’il s’agit d’effectuer un zoom :

\begin{figure}[!h]
\centering\fbox{\includegraphics[width=\linewidth, height = 8cm, keepaspectratio]{airepara}}%
\caption{Diagramme de cas d'utilisation global initial
\label{fg-aire}}
\end{figure}

\end{document}

(replace "airepara" by any graphics file that TeX can fing on your system).


text:
# trapezoid

## Primary tabs

\documentclass{article}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{pstricks}
\usepackage{psfrag}
\usepackage{graphicx}
\usepackage{amsthm}
%%\usepackage{xypic}

\begin{document}

A \emph{trapezoid} is a quadrilateral with at least one pair of \PMlinkescapetext{opposite} sides parallel.  Some \PMlinkescapetext{sources} insist that trapezoids have exactly one pair of \PMlinkescapetext{opposite} sides parallel, in which case parallelograms are not trapezoids.  Other sources do not restrict the definition in this manner, in which case parallelograms are trapezoids.  The convention in PlanetMath is to use the unrestricted definition.

In some dialects of English (\PMlinkname{e.g.}{Eg} British English), a trapezoid is referred to as a \emph{trapezium}.  Unfortunately, some confusion arises when this word is used, since in other dialects of English (e.g. American English), a \emph{trapezium} is a quadrilateral without any parallel sides.

Below is a picture of a trapezoid.

\begin{center}
\begin{pspicture}(0,0)(4,2)
\pspolygon(0,0)(1,2)(2.5,2)(4,0)
\end{pspicture}
\end{center}

The \emph{bases} of a trapezoid are its two parallel sides.  (If the trapezoid is a parallelogram, either pair of parallel sides can be declared to be its bases.)  The \emph{legs} of a trapezoid are the two sides that are not bases.  A \emph{height} of a trapezoid is a line segment that is perpendicular to the bases of the trapezoid and whose endpoints lie on the two lines formed by extending the two bases.  Typically, heights are drawn so that they intersect at least one base of the trapezoid.  (For some trapezoids, it is impossible to draw a height that intersects both bases.)  Below is a picture of a trapezoid with its bases labelled $b_1$ and $b_2$ and a height drawn in blue.

\begin{center}
\begin{pspicture}(0,-1)(6,4)
\psline[linecolor=blue](5,0)(5,3)
\rput[l](0,0){.}
\rput[r](6,3){.}
\pspolygon(0,0)(4,3)(6,3)(3,0)
\rput[a](1.5,-0.3){$b_1$}
\rput[b](4.5,3.2){$b_2$}
\psline[linestyle=dashed]{->}(3,0)(6,0)
\end{pspicture}
\end{center}

The \emph{median} of a trapezoid is the line segment whose endpoints are the midpoints of the legs of the trapezoid.  Below is a picture of a trapezoid with its median drawn in red.

\begin{center}
\begin{pspicture}(0,0)(4,2)
\psline[linecolor=red](0.5,1)(4,1)
\pspolygon(0,0)(1,2)(4,2)(4,0)
\end{pspicture}
\end{center}

In the \PMlinkescapetext{remainder} of this entry, only Euclidean geometry is considered.

If a trapezoid has bases of lengths $b_1$ and $b_2$ and a height of length $h$, then the area of the trapezoid is

$$A=\frac{1}{2}(b_1+b_2)h.$$

Note that the length $m$ of the median of a trapezoid is the arithmetic mean of the lengths of its bases; \PMlinkname{i.e.}{Ie},

$$m=\frac{1}{2}(b_1+b_2).$$

Thus, the area of a trapezoid can also be determined by

$$A=mh.$$
%%%%%
%%%%%
nd{document}


text:
Is there a standard way to read the Legendre symbol?

Is there any standard way to read (in English) the Legendre symbol (a|p) (http://en.wikipedia.org/wiki/Legendre_symbol), say, similar to "a choose b" which is used for the binomial coefficients?

• I wouldn't mind knowing how this is read in other languages, by the way :) Feb 16, 2010 at 15:02
• I am almost certain that the answer to your question is "no", and for everyone that posts an answer saying "well I say it this way" this is in some sense more evidence that the answer really is "no". Isn't that strange? Feb 16, 2010 at 15:23
• I think Kevin is correct. (As for myself, I say "the Legendre symbol of a with respect to p", because I don't know what else one is suppose to say.) Feb 16, 2010 at 15:34
• Mariano: In German,(a|b) is read "a nach b", according to a German mathematician who asked me the original question. Feb 16, 2010 at 16:02

• To use this in the statement of quadratic reciprocity, you'd have to talk about QR $p$ $q$ and QR $q$ $p$, which don't sound good to me. (An issue that came up in my thesis: If $\ell$ and $p$ are already taken for primes, and $q$ is already taken for a prime power, what is the notation for anotheer prime?) Jun 14, 2011 at 22:07


text:
# Proof of Euclidian Algorithim from Terrence Tao's Analysis 1

In Analysis 1 in chapter two section 3 on multiplication one of the exercises is the proof of the following.

(Euclidean algorithm). Let n be a natural number, and let q be a positive number. Then there exist natural numbers m, r such that 0≤r < q and n = mq +r.

With the hint to induct on n. I am having trouble navigating this proof with only without the idea of subtraction and only cancellation laws for addition and multiplication.

A proof would be helpful for me to reference.

• That's not what's usually called the Euclidean Algorithm, it's usually called the Division Theorem or the Division Algorithm. The Euclidean Algorithm is the one that calculates greatest common divisors (by repeated application of the Division Theorem). If you search this website for "Division Algorithm", you will probably find an answer to your question. – Gerry Myerson Jul 21 at 2:23
• @Gerry Try googling "Euclidean division" and "Euclidean division algorithm", e.g. Hendrik Lenstra wrote "Kummer noted that it would be sufficient to prove uniqueness of factorisation and for this he turned to the method that also Wantzel would apply: the euclidean division algorithm." – Bill Dubuque Jul 21 at 3:19
• We can extend $\Bbb R$ to a larger ordered field $\Bbb R^*$ in which all the basic rules for $+,-,\times, /$ and $<$ still apply, but $\Bbb R^*$ has (among other things) a member $q$ which is larger than every member of $\Bbb Q.$ The Q is not true if "number" means member of $\Bbb R^*$.... If "number" means member of $\Bbb R$ we must use some consequence(s) of the $definition$ of $\Bbb R$ beyond the basics of arithmetic to ensure we are not talking about $\Bbb R^*....$ Topics: Field, Ordered field, Complete ordered field, Archimedean property. – DanielWainfleet Jul 21 at 9:07
• @Daniel, I suspect $q$ is meant to be a positive integer. – Gerry Myerson Jul 21 at 9:13
• @Bill, thanks, but I'm specifically questioning "Euclidean algorithm", and nothing in your comment refers to that phrase. – Gerry Myerson Jul 21 at 9:15

Note that the base case is trivial, that is, for $$n=1$$ and any positive integer $$q$$ we have that $$n=0.q+1$$ if $$q>1$$ ($$n=1.q+0$$ otherwise).

Assume that the claim holds for every natural number $$n\leq N$$. (To be more explicit, it means that given a natural number $$n\le N$$ and a positive integer $$q$$ we have $$n=mq+r$$ such that $$0 \le r < q$$.)

We need to show that the claim also holds for $$n=N+1$$. To this end, recall that $$N=mq+r$$, $$0\le r < q$$ therefore $$N+1=mq+(r+1)$$. Now, analyse the following two cases separately: if $$r then we are already done. Otherwise, $$r=q-1$$ then $$r+1=q$$. Therefore, $$N+1=mq+q=(m+1)q+0$$. This proves that the claim holds for $$N+1$$.

P.S.: I do not have the book you have mentioned with me, so I am not sure how the set of natural numbers is defined. It may so happen that $$0$$ is included in the set or natural numbers. Most part of the proof still goes through, only the base case needs to be modified. In fact, it only makes it simpler.

Hint: the induction step is like adding $$1$$ in radix $$q$$ arithmetic. The analogy is clearer with $$2$$ digits $$r_i$$.

By induction $$\,n\!-\!1\, =\, m\, q^2 + r_1 q + r_0,\$$ all $$\,r_i < q$$

when $$\,r_0 < q\!-\!1\,$$ then adding $$\,1\,$$ yields $$\,r_0 < q\,$$ so we're done.

else $$\ \ \ r_0 = q\!-\!1\,$$ then adding $$\,1\,$$ we get a carry into $$\,r_1,\,$$ handled as above, i.e.

when $$\,r_1 < q\!-\!1\,$$ then adding $$\,1\,$$ yields $$\,r_1 < q\,$$ so we're done,

else $$\ \ \ r_1 = q\!-\!1\,$$ then adding $$\,1\,$$ we get a carry into $$m$$ and we're done (since $$m$$ can be any natural).

So it is essentially a radix $$q$$ increment operation restricted to the least $$2$$ digits (with all higher digits collected into $$m)$$. In the OP we use only the least digit (vs. least two).


text:
# iLand Tech blog Blog Actions

## let the seeds fly

Friday 30 of April, 2010

During the last weeks we spent some time discussing the whole topic of regeneration. We are quite excited about that, because this will be the first real 'landscape'-process and we are curious to find out how nice this plays with the concepts on which iLand is built upon... Besides that, Rupert did a very thorough testing and parameterization of PNW tree species. To that end we applied some modifications in the water cycle (as it seams, it is not easy to model the severity of the growth limitation due to drought). Anyway, so far the results look quite good.
But this post is about the seeds. The seed dispersion algorithm was the first part of the regeneration chain that made the jump from mere theory (e.g. regeneration, or dispersal) into silicio. Implementation details may be found here.

Here is an image showing the dispersal:
<img src='tiki-view_blog_post_image.php?imgId=3' border='0' alt='image' width='100%'/>

The test landscape is essentially empty (see the small screenshot in the lower right corner) and so the effect of the dispersal process is quite obvious.

The production of such nice looking debug images is very straightforward and can be done by one line of code:

// save a image in black/white mode (2nd parameter)
// store as png
gridToImage(seedMap(), true, 0.,1.).save("targetfile.png");

The (iLand)-Grid class offers a 'gridToImage()' function; and saving as a compressed .png is done - thanks to Qt - just by a call to 'save()'.

The next challenging topic will be to find way to model the establishment phase in a way that is both scientifically sound and computationally cheap... The answer (of course) is blowing in the wind...


text:
# "APIO2012" and "Luogu P1552" dispatch

### Description

In this gang, there is a ninja called Master. Except Master, every Ninja has only one superior. In order to keep secrets and enhance the leadership of ninjas, all instructions related to their work are always sent by their superiors to their direct subordinates, and are not allowed to be sent by other means.

Now you need to recruit a group of ninjas and send them to customers. You need to pay a certain salary for each dispatched ninja, and make the total salary not exceed your budget. In addition, in order to send instructions, you need to select a ninja as the manager. The manager is required to send instructions to all dispatched ninjas. When sending instructions, any Ninja (whether dispatched or not) can be used as the messenger. The manager may or may not be assigned. Of course, if the manager is not dismissed, you don't have to pay the manager's salary.

Your goal is to maximize customer satisfaction within your budget. Here, customer satisfaction is defined as the total number of dispatched ninjas multiplied by the leadership level of the manager, and the leadership level of each Ninja is also certain.

Write a program to give the superior \ (b\u i\), salary \ (c\u i\), leadership \ (l\u i\), and the total salary budget paid to the Ninjas \ (M\), and output the maximum value of customer satisfaction when the above requirements are met within the budget.

### Hint

• $$1 ≤ n ≤ 10^5$$ number of ninjas;

• $$1 ≤ m ≤ 10^9$$ total salary budget;

• $$0 ≤ b_i < i$$ the number of Ninja's superior;

• $$1 ≤ C_i ≤ M$$ Ninja's salary;

• $$1 ≤ l\u I ≤ 10^9$$ Ninja leadership level.

• For the data of the previous \ (30\%\), \ (n ≤ 3\times 10^3\).

### Solution

Obviously, the working relationship of Ninja forms a tree structure.

Then the problem can be formalized as follows: select as many nodes as possible in each subtree so that the total cost of these nodes does not exceed \ (m\). Finally, merge the answers of all subtrees.

It is not difficult to think of starting from the leaf node and maintaining the Ninja set with an appropriate data structure.

Whenever a point is reached, first merge the current set to the superior, and adjust the superior set so that the Ninja cost no more than \ (m\) and the set is as large as possible.

We greedily think that as long as we pop up the Ninjas that cost a lot each time, because the answer only involves the leadership of the roots and the number of ninjas.

After one point is completed, you can directly use your superior to update the answer.

For collection maintenance, we need an appropriate data structure that supports deletion of maximum values and efficient consolidation. Here, select the left leaning tree.

As for the complexity, since each Ninja will only pop up once, the time complexity is \ (O(n\log n) \).

### Code

#include <algorithm>
#include <iostream>

using namespace std;
const int N = 1e5 + 5;

struct ninja {
} nj[N];
int n, m;

struct ltNode {
int ch[2], val, dist, rt;
} t[N];
#define lc t[x].ch[0]
#define rc t[x].ch[1]

int merge(int x, int y) {
if (!x || !y) return x | y;
if (t[x].val < t[y].val) swap(x, y);
rc = merge(rc, y);
if (t[lc].dist < t[rc].dist) swap(lc, rc);
t[x].rt = t[lc].rt = t[rc].rt = x;
t[x].dist = t[rc].dist + 1;
return x;
}
int findrt(int x) {
return x == t[x].rt ? x : t[x].rt = findrt(t[x].rt);
}

struct LefTr {
int root;
int size;
long long sum;
inline void join(LefTr t) {
size += t.size, sum += t.sum;
root = merge(root, t.root);
}
inline void pop() {
sum -= t[root].val;
root = merge(t[root].ch[0], t[root].ch[1]);
--size;
}
} lt[N];

void init() {
for (register int i = 1; i <= n; i++) {
lt[i].root = t[i].rt = i;
lt[i].size = 1;
lt[i].sum = t[i].val = nj[i].cost;
}
}

signed main() {
ios::sync_with_stdio(false);
cin >> n >> m;
for (register int i = 1; i <= n; i++)
cin >> nj[i].fa >> nj[i].cost >> nj[i].lead;
init();

long long ans = 0ll;
for (register int i = 1; i <= n; i++)
ans = max(ans, nj[i].lead * 1LL);

for (register int i = n; i > 1; i--) {
int f = nj[i].fa;
lt[f].join(lt[i]);
while (lt[f].sum > m) lt[f].pop();
ans = max(ans, lt[f].size * 1LL * nj[f].lead);
}

cout << ans << endl;
return 0;
}

Tags: data structure

Posted by Penelope on Tue, 31 May 2022 12:43:21 +0530